domain: phys.cond spi phys chim

text:
Baptiste Bienvenu

Hakim Amara

François Ducastelle

Lorenzo Sponza

Strain and electronic properties at the van der Waals interface of phosphorus/boron nitride heterobilayers

Past years saw a massive investigation of atomically thin materials, or 2D materials, demonstrating both theoretically and experimentally many of their interesting electronic and mechanical properties [1][START_REF] Xia | Ramasubramaniam[END_REF][3][4] . In 2013, A. Geim and I. Grigorieva opened new perspectives in this field by proposing the concept of van der Waals heterostructure 5,6 , i.e. a compound formed of different monolayers stacked on top of each other, as a way to combine and control the characteristics of several 2D materials on a single system. Since then, such structures have been investigated theoretically and experimentally, and lead to some breakthrough applications and insights [7][8][9] .

Among 2D materials, black phosphorus (BP) exhibits anisotropic mechanical, electronic and optical properties due to its relatively flexible "accordion" structure. Other peculiarities are a very high hole mobility along the armchair direction 3,4,10,11 and, quite uniquely within the 2D family, its fundamental gap remains direct all along the way from the bulk down to the monolayer. As in all 2D materials the layer properties are particularly sensitive to the surroundings, for example the gap width can be modulated on a large range by varying the dielectric properties of the environment 12 or the number of stacked layers [12][13][14] (from 0.3 in the bulk to more than 2.0 eV in the monolayer). All these characteristics make BP a promising candidate for ultra-flat and flexible opto-electronic technology. It has been synthesized using various techniques such as mechanical or chemical exfoliation 3,15,16 , or by Chemical Vapor Deposition 17,18 . However, at low thickness, its chemical stability is compromised by its high reactivity 16 , which implies encapsulation or passivation for any practical application.

Another 2D allotrope of phosphorus is the blue phosphorus (P blue ), which has been predicted to be stable by theoretical works 19 and has been produced by molecular-beam epitaxy 20,21 . It is homostructural to silicene, with a buckled honeycomb lattice. Its electronic properties differ from those of BP, with a larger indirect bandgap, similarly modified by the layer number 19 .

One way of protecting both the BP and the P blue is to cap them with passivation layers. But owing to their sensitivity to the external conditions, it is important to choose capping layers that either preserve target properties, or modify them in the desired way. To this purpose, hexagonal boron nitride (hBN) was proposed as a good candidate because of its chemical inertness and its insulating properties 22,23 which are predicted to preserve the electronic features of phosphorus layers upon stacking [24][25][26][27][28][29] . However, due to lattice mismatch and symmetry differences between the two layers, contact strain cannot be avoided. This is is likely to have an impact on both the structural stability and electronic properties 30 . Since previous studies on P/BN bilayers have taken into account structural relaxation only partially or lack of comparative studies, we think that it is worth carrying out a more systematic investigation. Complete studies of this kind have been carried out on different systems by Van Troeye and coworkers 31 on bulk heterostructures of BP on graphene, and by Peng, Wei and Copple 32 on monolayer BP.

In section I, we present the structures considered in this work. In section II we clarify our methodology and define useful quantities. The structural properties of the BP/hBN and the P blue /hBN heterostructures are presented and discussed in section III. Section IV is devoted to the discussion on how the electronic properties of the constituent layers are modified by the contact deformations. Conclusions will be presented in section V.

I. CONSIDERED STRUCTURES

The building-block monolayers are reported in Figure 1. The structural parameters reported have been obtained after full relaxation (atomic positions and cell parameters) of the unitary cells. We will refer to the two crystallographic directions of the BP cell using the common "armchair" and "zigzag" nomenclature. These directions are marked with blue arrows in Figures 1 and2. In constructing the BP/hBN heterobilayers, it is convenient to consider the non-primitive orthorhombic cell of hBN, in which b hBN = √ 3a hBN , which is drawn with a dashed red rectangle in Figure 1b).

We considered four different heterostructures, three of the BP/hBN type and one of the P blue /hBN type. The corresponding equilibrium structures are visualized in Figure 2 and their structural parameters are reported in Table I. These structures have been chosen to be good compromises between computational feasibility and size.

We start presenting the three BP/hBN systems presented in Fig. 2 found equilibrium distance d between the layers is 3.36 Å, which falls inside the interval of values reported by other studies done in the same cell 24,26,28 . In this cell, only tensile and compressive stress are exerted on each layer. Most of BP/hBN or BP/graphene simulations found in literature 24,26,28,29,31,34 are done in the C-cell because it is the smallest supercell accommodating the two layers without excessive deformations. Still, strain can be quite relevant in this configuration, so it is sensible to look for larger supercells where layers are less strained.

In the S-cell, the BP layer is sheared of 3 • , so that the angle γ between the armchair and the zigzag directions is equal to 93 • (see Fig. 2). This deformation allows us to accommodate 4 BP cells on top of 11 B-N pairs, for a total of 38 atoms, without changing significantly the cell parameters of either monolayer with respect to its free-standing structure. As a result, uniaxial strains are quite low in each sheet, at the price of applying a shearing stress and breaking the orthorhombic symmetry of the BP layer. The lattice parameters of the S-cell at equilibrium are a=6.62 Å, b=9.04 Å, c=2.13 Å , and d=3.43 Å. This type of cell, already considered in literature 35 , is compared here with other BP/hBN configurations. Moreover, shearing the BP layer without deforming the underlying hBN, seems quite a difficult condition to realize in actual experiments.

In the T-cell, strain is minimized through a twist angle of 19.14 • between the orthorhombic cells of the two layers. In this configuration, it is possible to fit 10 unitary cells of BP on top of 28 B-N pairs, for a total of 96 atoms. The structural parameters of the supercell at equilibrium are a=6.62 Å, b=22.95 Å, c =2.11 Å, and d =3.43 Å. As in the C-cell, no other deforming force is exerted on the layers of the T-cell apart from tensile and compressive stress. The same twist angle has been investigated by Van Troeye and collaborators 31 who scrutinized also other twisted configurations, however the systems they studied are BP/graphene bulk heterostructures, so quite different from ours in particular because of the absence of polar atomic bonding in graphene. Finally, Costantinescu and coworkers 27 have investigated similar systems (hBN/BP/hBN and BP/hBN/BP) with a twist angle of about 40 • between the layers, but the study lacks comparison with other structures.

The only cell type considered in the P blue /hBN bilayer is reported in Figure 2 and is called H-cell because of its in-plane hexagonal symmetry. It is constructed by fitting 9 unit cells of P blue on top of 16 unit cells of hBN, for a total of 50 atoms. In the resulting configuration, only biaxial tensile or compressive stresses are applied on the layers. Its equilibrium lattice parameters are: a=9.99 Å, c=1.23 Å , and d=3.37 Å.

II. COMPUTATIONAL DETAILS AND METHODS

All first-principle calculations in this paper are done within density functional theory (DFT) with the ABINIT simulation package 36 . The Kohn-Sham states are represented with a plane-wave basis with a cutoff energy of 60 Ha. Norm-conserving pseudopotentials 37 are used, including 3 electrons for B and 5 electrons for N and P. The exchange-correlation potential is approximated using the GGA-PBE functional 38 . In addition, van der Waals forces are taken into account with the vdW-DFT-D3 scheme introduced by Grimme et al. 39 . The in-plane k-point sampling used in the density calculations has one point along z and 5x3 in the plane of the C-cell, 4x4 in the S-cell, 3x5 in the T-cell and 5x5 in the H-cell. To avoid artifacts due to interacting periodic replicas along the z axis, a vacuum of 16 Å is included in all supercells.

The heterobilayers considered in this study can be divided into four groups: the C-cell, S-Cell, T-cell and Hcell. In order to analyse the energetic costs and gains to first deform and then stack the layers, we employed the structural optimization procedure sketched in Figure 3.

1. The first step consists of varying by hand the size of the cell while relaxing only the atomic positions. In doing this, we took care of preserving the hexagonal symmetry of the hBN layer, which fixes a b/a ratio in each cell. This procedure allowed us to draw stability curves of the bilayers, i.e. total energies as functions of the surface S. Successively we fitted these curves with the quadratic function

E(S) = E(S 0 ) + B 4S 0 (S -S 0 ) 2 (1)
where E(S) is the total energy at surface S, S 0 the surface which minimizes the energy fit, and B = S 0 d 2 E/dS 2 the uniaxial or biaxial 2D compressibility depending on the nature of the strain. We prefer expressing this quantity in meV/ Å2 instead of GPa, as commonly done because that would imply a normalisation on the thickness of the layer, which is a quite arbitrary concept we want to avoid. The minimum of the stability curve E(S 0 ) constitutes the starting point for the next step. In the case of the C-cell, this procedure has been initialized with eight inequivalent configurations, each corresponding to one way of stacking the four atoms of a unitary cell of BP on top of the two atomic species of hBN. For the other cells, the large size of the calculation hampered the exploration of several initial configurations.

2. Once S 0 is extracted from the fit, we allow for a full structural relaxation (atomic positions and cell parameters), with the intent to account for possible deviations from the fitting functions. At this step, the hexagonal symmetry of hBN may break due to internal stresses. The resulting structures, the parameters of which are reported in Table I, are our best estimate of the equilibrium configurations of the heterobilayers. Note that in the T-cell a full relaxation has not been possible owing to its large dimensions. Instead, cell parameters have been derived from the fitted S 0 and the hexagonal symmetry of hBN (b/a = 3.464) and then the relaxation has been done only for the atomic positions.

3. Successively, from the equilibrium heterostructures, we remove all the P atoms and we compute the total energy of the hBN which is left without performing any relaxation. So the hBN layer is isolated but its atoms occupy the same position as in the equilibrium heterostructure (deformed monolayer). Then we do the same by keeping only the phosphorus atoms and removing the hBN. In this way we can access the total energy of a hypothetical stage where each building-block is isolated but experiences the same deformations as in the bilayer.

4. Finally, fully relaxed calculations of the monolayers in their unitary cells gave us access to the total energy of the free-standing pristine materials (relaxed monolayers).

We use the FIRE 40 algorithm for the first optimization steps, and successively the BFGS 41 as they are implemented in ABINIT. The stopping condition for all relaxation calculations is that the change of the forces exerted on every atom must be lower than 10 -5 Ha/Bohr along all Cartesian directions. Within the heterobilayer, the strain percentage on the generic lattice parameter p of the material M (M=BP, P blue or hBN) is defined

η(p M ) = 100 p M -p REF M /p REF M , (2)
where p REF

M

is the corresponding cell parameter of the relaxed monolayer M.

III. STRUCTURAL PROPERTIES

In this section, we first investigate the influence of induced deformations on the total energy of the three iso- lated monolayers, then we pass to the discussion of the heterobilayers.

A. Isolated monolayers

The equilibrium structure of the monolayers, reported in Figure 1 have been obtained using a full-relaxation procedure (step 4). In addition, starting from each equilibrium structure, we changed the cell parameters and let the atomic positions relax in order to compute total energy variations as functions of the resulting strain. This allowed us to draw the stability curves reported in Fig- ures 4, from which we derived the 2D compressibility B of the monolayers by fitting equation 1 (see Table II and Figure 4). Our stability curves of BP and hBN are in good agreement with other published results 32,42 .

In the case of hBN (red diamonds) and P blue (blue downward-pointing triangles) we have applied a biaxial stress to preserve their hexagonal symmetry, while in the case of the orthorhombic BP layer, we have varied the a (zigzag) and the b (armchair) parameters separately, 2D compressibility and strain distribution on the constituent layers. In the H-cell biaxial strain is along both a directions.

* Because of the misalignment in the S-cell and the T-cell, the strain distribution in hBN is always given with respect to the armchair and zigzag directions of the BP layer (cfr Figure 2). Note that they are not perpendicular in the S-cell.

which corresponds to applying uniaxial stress along either direction. In this study, we are not interested in biaxially strained BP because, once inside the bilayer, its deformations are anisotropic. For a study on the effects of biaxial strain on BP, we refer to the work by C ¸akır and coworkers 43 . Finally, we have looked at the stability curve of the sheared BP monolayer (black upward-pointing triangles), where we decided to apply stress in such a way to keep a shear angle of 3 • . According to our analysis, the hBN is clearly the stiffest layer of all owing to its flat hexagonal structure which hinders any deformation. As expected 31,32 , the BP structure is found to be much softer along the armchair direction than along the zigzag one. The P blue layer is also quite soft because the in-plane stress can be partially transferred to the buckling parameter h, which passes from 1.29 Å when a is compressed of 6% to 1.20 Å upon a dilatation of the same amount.

As a result, in the P/BN heterostructures we expect strain to be larger on the phosphorus sheet (BP or P blue ) which will deform conveniently to adapt to the much more rigid hBN layer.

B. BP/hBN bilayer

In this section we discuss the effects of strain on the three BP/hBN cells. Results are summarised in Table III.

Cell type C

As explained already in the methodological section, in the case of the C-cell type, we have been able to carry out the relaxation procedure (step 1) starting from eight different configurations. Half of them correspond to the four possibilities to place the inequivalent P atoms of BP on top of a N atom placed at the origin (P/N configurations). In the other four configurations, the same P atoms are on top of a B atom placed at the origin (P/B configurations). Note that other configurations were possible, for example we could have put the P-P bonding on top of B or N, but configurations alike have been predicted to be less stable 24 . Moreover, they can be obtained from those we considered through a rigid translation in the armchair direction of the BP, but such a modification has been demonstrated to be energetically expensive 26,29,34 .

The resulting eight stability curves (not shown) have the same shape and all relax to the same cell-parameter with negligible total energy differences (order of 0.1 meV). However, it is possible to identify some structural selection mechanism since, at the end of the relaxation procedure, all the P/N cells have evolved into P/B ones. This tendency is rationalized observing that in the B-N ionic bonding, trivalent B atoms have a positive charge, whereas pentavalent N atoms are negatively charged. As P is also pentavalent, and hence attracted by B atoms. This repulsion is obviously not present in BP/graphene 31 and is likely related to the polar nature of the B-N bonding.

All further analysis on the characteristics of the C-cell have been carried out on the structure with the lowest total energy, chosen as representative in virtue of the very tiny energy and structural differences among them. In Figure 5, black squares mark the computed total energy The red bullet at the origin marks the equilibrium structure of each system.

as a function of the surface variation ∆S of this representative C-cell. In dashed black, the fit of expression 1 allowed us to extract the biaxial 2D compressibility B of the bilayer, as explained in Section II (step 1).

In order to form the heterostructure, both layers deform. In particular, the BP layer contracts along the armchair axis (-5.32%) and elongates along the zigzag axis (+1.44%) and the buckling (+0.77%). Similarly, the hBN layer undergoes an anisotropic deformation but with opposite sign, with a weak compression (-0.16%) along the zigzag axis aligned to the BP zigzag one, and a weak dilatation (0.29%) of the armchair axis perpendicular to the former. As expected, the preferential direction to release the stress is by far the relatively soft armchair axis of BP, while the hBN layer acts basically as a rigid support.

Cell type S

Data from the S-cell are reported in red triangles in Figure 5. The 2D biaxial compressibility of the S-cell is almost the same as that of the C-cell.

In this configuration, the strain is actually evaluated along directions that are not Cartesian, the shear angle between the BP (reference) armchair and zigzag direction being 3 • . Moreover in the hBN layer these two directions have no special interest (neither is the zigzag or armchair axis of hBN). So the analysis of the strain is somewhat more involved and results should be interpreted as qualitative indications. What can be said, however, is that also in this case most of the stress applies along the armchair axis of the BP layer which compresses to accommodate to the underlying hBN sheet, even though strain is much lower than in the C-cell, which is the reason it has been studied in literature 35 . Similarly to the previous case, the hBN layer plays the role essentially of a rigid substrate, undergoing only a very weak and anisotropic deformation with a different sign with respect to that of the BN layer.

Cell type T

Because of the large dimension of this cell-type, once we extracted the equilibrium surface S 0 from the fit, we could not make a full relaxation to permit anisotropic deformations of the hBN layer. So we constructed a cell of surface S 0 with fixed ratio b/a = 3.464, which keeps the hexagonal symmetry of the hBN layer (cfr. Table III), and we relaxed only the atomic positions. From this equilibrium position we changed the cell parameters to construct the energy curve, reported in green circles in Figure 5. Thanks to the twist angle between the BP and the hBN, this structure minimizes the deformations on the BP layer while keeping the deformation of hBN quite low (hBN: -0.14% biaxial; BP: -0.2% on armchair, 0.71% on zigzag and 0.16% on buckling). In addition it is found to be the most rigid of all bilayers.

In the work of Van Troeye and collaborators 31 on BP/graphene heterostructures, the BP layer always deforms by contracting along the armchair axis, so the authors suggest that this is a general behaviour of BP under stress. This is indeed what we observed in the S-cell and the C-cell. However, in the T-cell, most of the deformation undergone by the BP layer is a dilatation of the zigzag axis, despite its higher stiffness, whereas the armchair axis contracts only by 0.2%. This result can be ascribed to two differences with respect to the work of Van Troeye. First, our BN layer is less homogeneous than the graphene layer because of the polar nature of the B-N bonding, which acts on a short-range scale, probably through the same P/N repulsion we observed while relaxing the C-cells. We expect this characteristic to cause an angle-dependent deformation behavior of the deposited BP layer. Second, Van Troeye's systems are bulk heterostructures, not isolated bilayers. This can also influence the atomic arrangement through long-range interactions, even though we estimate these long-range contributions to be less important than the short-range ones. Unfortunately, our computational resources could not allow us to investigate further these hypotheses.

H-cell

The data relative to the H-cell are reported in Table III and as blue diamonds in Figure 5.

In relaxing this structure we applied biaxial stress on the bilayer supercellmaking the hypothesis that the hexagonal symmetry is kept in deformation because it is common to both layers. The results confirm once more that the hBN has the tendency to contract slightly (-0.4%) while the P layer undergoes much larger deformations to adapt to the underlying sheet. In this case, the biaxial in-plane dilatation of the P blue layer (+1.28%) is accompanied by a sizable reduction of the buckling of 2.18%.

Stability and energy costs

Our relaxation procedure allowed us to quantify the stability of the bilayers, while distinguishing how much energy is employed to deform the layers and how much is gained in the adhesion, as depicted schematically in Figure 6.

In order to compare the stability of the different bilayers, we introduce the formation enthalpy per unit surface

∆H = E P/BN -N P E P -N BN E BN 2S P/BN = ∆E 2S P/BN , (3)
where N M is the total number of atoms belonging to the P or the BN layer in the unitary cell of the bilayer, E M is the total energy per atom of the isolated monolayer of the material M, E P/BN is the total energy of the unitary cell of the bilayer and S P/BN its surface. The formation enthalpy is negative if the structures is stable and the lower it is, the more stable is the bilayer. It is simply related to the binding energy per P atom E b = ∆E/N P used in similar studies 24,26,28,29,34 and to the cohesion energy per atom E coh = ∆E/(N P + N BN ) employed by Van Troeye et al. 31 .

The number of atoms, the bilayer surface and the formation enthalpy are reported in the top part of Table IV. All bilayers are predicted to be stable. Among the BP/hBN cells, the S-cell exhibits clearly the weakest cohesion. On the other hand, the H-cell has the lowest formation enthalpy among all considered bilayers, comparable to that of the S-cell, despite the same symmetry of the two layers, although this result may result from the specific choice of the stacking. The difference between the enthalpy of the C-cell and the T-cell is very low so we are not in the condition of claiming that one is particularly more stable than the other. This result marks a difference from the conclusions of drawn for hBN/graphene bulk heterostructures 31 , where the C-cell was clearly favoured. As discussed earlier, this different conclusion can come from two factors: one connected to the weak polar nature of the B-N bonding and to the resulting P-N repulsion we observed, and the other to the different geometries considered in the two studies (isolated bilayers here, bulk heterostructures in 31 ).

For comparison, we checked that our binding energy of the C-cell (E b = -62.2 meV/P) falls in the interval of values reported in literature 24,26,28,29 . These range from a maximum of -84 meV/P 29 to a minimum of -49 meV/P 26,44 .

We analysed in more details the energy cost of the evolution from the monolayers to the bilayer by splitting the formation enthalpy into three contributions ∆H = ∆E P d + ∆E BN d + ∆E a . To this end, we introduce the deformation energy per unit surface of the sheet (M=P or BN)

∆E M d = E M d -N M E M /2S P/BN , (4)
defined as the energy difference between the isolated monolayer deformed as in the bilayer (E M d ) and the fully relaxed monolayer, divided by twice the surface of the bilayer. This quantity expresses the energy paid per unit surface to deform each monolayer to the same configuration it will have inside the bilayer. The adhesion energy

∆E a = E P/BN -E P d -E BN d /2S P/BN (5)
is defined as the energy cost per unit surface to stack the two monolayers once they have the appropriate atomic configuration.

The comparison between ∆E a and ∆E M d tells us the importance of taking into account deformations. Among the BP/hBN cells, the configuration where deformations cost the most is the widely studied C-cell, while that where they cost the least is the T-cell. In absolute terms however, it is in the H-cell where the deformation energy is minimized. More details can be understood if looking at the strain distributions reported in Table III. In all BP/hBN structures, the deformation energy concentrates essentially on the softer BP layer, but not with the same effects. Indeed in the C-cell, with the largest deformation energy, such a cost is explained by the amplitude of the compression of the armchair axis (-5.32%). In the S-cell the amount of axial strain is substantially lower, so most of the deformation energy is used to shear the structure. The result is that, despite the lower axial strain, the deformation energy of the S-cell is comparable to that of the C-cell. Finally, in the T-cell, the quite low deformation energy (less than half of that of the S-cell) is spent mostly to elongate weakly the relatively rigid zigzag axis. The low price of its deformations is the reason for the higher stability of the T-cell.

In fact, the energy hierarchy is actually inverted once the two layers are stacked. The lowest ∆E a is the one of the C-cell, while the T-cell has the highest one among the BP/hBN family. Despite this, the respective deformation energies are such that the T-cell turns out to be the most stable structure, with a ∆H comparable to that of the C-cell. The S-cell instead, having already the worst ∆E a results in being by far the less stable among the BP/hBN bilayers.

C-cell S-cell T-cell H-cell N P 12
Concerning the H-cell, it exhibits the worse ∆H despite the quite low deformation energies, comparable to those of the T-cell. Similarly to the S-cell scenario, this has to be ascribed to the bad ∆E a , more than to the deformation term since this is as low as in the T-cell. Quite interestingly, the higher rigidity of the P blue layer forces the hBN to deform more than in the BP/hBN heterostructures. However, in the case of the P blue /hBN structures, we stress that a single configuration is not enough to generalise this result and it could be that other P blue /hBN bilayers exhibit higher stability.

IV. ELECTRONIC PROPERTIES

In previous works on BP/hBN heterostructures 24,26,28,29 , the bilayers were always in the C-cell configuration, which by itself imposes a lot of stress onto the BP layer, even when relaxed, as we have shown in the previous section. Shrinking the hBN layer instead of the BP one has also been considered 29 , but we have shown that this is not a realistic account of the strain conditions in this system. These works all agree on the fact that hBN is a promising candidate for protecting BP while preserving its peculiar electronic properties, in particular its high hole mobility in the armchair direction 4 . However, owing to the extreme sensitivity of BP properties upon deformation 32 , we think it is worth reconsidering this statement on the light of the results we presented above on the distribution of strain.

To this end, the electronic properties of the bilayers are discussed on the basis of the strain distribution of their constituent sheets. In Table V we report the PBE energy gap of the four heterostructures. It is well known that DFT systematically underestimates the band gap of semiconductors and its value strongly depends on the choice of the exchange-correlation potential. However differences are quite reliable within the same approximation.

A. The BP/hBN bilayers

Previous studies carried out on the C-cell 26,28,29 , found that the top valence and the bottom conduction states were localised on the BP layer. In addition all these works report that the bilayer gap is between 0.2 and 0.07 eV larger than that of the isolated BP because of local electrostatic effects 28 . We confirm that this is the case also in the other two cells. As our PBE gap of the isolated BP monolayer is 0.81 eV (in quite good agreement with similar calculations 26,28,[45][46][47] ), we also confirm that there is a general trend of opening the gap in the bilayer configuration, even though not particularly in the C-cell.

We have shown that, in all BP/hBN bilayers, the BP sheet undergoes most of the deformations, preferentially along the armchair direction in the C-cell and the S-cell, and along the zigzag axis in the T-cell. In Figure 7 we report the band structure of the BP monolayer under uniaxial strain along the two crystallographic directions. Our results are consistent with the very complete study of Peng and coworkers 32 . The band plot shows that contractions along the armchair axis may spoil the remarkable transport properties of BP because of the rise of an alternative top valence state in the Γ -X (armchair) direction when strain is close to -6% (blue arrow in Figure 7a). This is the scenario we encountered in the C-cell and, in a lesser measure, in the S-cell. Note that, despite the shearing angle, the features of the band structures that we discussed, look similar in the C-cell and the Scell and evolve in a similar way. This consideration contradicts partially previous works on the subject, and sug- FIG. 8. Density of states of the relaxed monolayers (top panel) and the three BP/hBN bilayers. In the bilayers, the dashed BP and the dotted hBN curves refer to the monolayer contributions to the total DOS, the thin red curve, is the deformed BP monolayer. The top valence of the bilayers and the BP monolayer are aligned at 0, while the top valence of the hBN monolayer is aligned approximately to the hBN contribution of the bilayers. All data have been broadened with a Gaussian of σ = 0.05 eV.

gests that BP/hBN bilayers where the zigzag axis of BP and hBN are aligned may lead to a degradation of the transport properties of BP because of the contact strain which leads to a substantial contraction of the armchair axis of BP.

Instead, in the T-cell most of the deformation is a dilatation along the zigzag axis of about 1.5%. Such a deformation has a minor impact on the armchair hole mobility since bands are quite unchanged along Γ -X. Instead the zigzag hole mobility may result slightly enhanced because the curvature of the top valence in Γ -Y increases as the zigzag axis is expanded (blue arrow in Figure 7b). On the other hand, at low degrees of tensile strain, the conduction band is a bit flattened, leading potentially to a lower electron mobility. But, if the zigzag axis is highly extended (around 6% of tensile strain), a band comes down from high-energy regions and becomes a new bottom conduction, as highlighted by the yellowish arrow. Given the high dispersion of this new conduction state, electrons should acquire a lighter effective mass and one should observe the electron mobility suddenly increase as the band crossing takes place. To generalise this statement we claim that, if some twist angles allow one to force BP to stretch along the zigzag axis instead of contracting along the armchair axis, that would not only preserve the very good hole mobility along the armchair direction, but even lead to a substantial increase of the electron mobility in the zigzag direction on the proviso that the elongation is sufficiently important (around 6%).

To characterize further the electronic properties of the three considered BP/hBN bilayers, we report the layerprojected density of states (LP-DOS) in Figure 8. Thick black lines report the total DOS of the bilayers, while dotted blue and dashed-dotted red curves report the contribution of the hBN and the BP layer respectively. In thin solid red, we also report the DOS of the isolated BP layer deformed as in the heterobilayer.

Comparing the gap of the structures, we identified a competition between the opening of the gap due to the interaction with the hBN and the closing because of the deformation. This is true in all systems, but it is particularly evident in the C-cell where the two effects almost cancel (0.81 eV in the relaxed monolayer, 0.80 in the deformed one, 0.83 in the bilayer). This explains partially why we observed a much smaller opening of the gap with respect to the other works done on the C-cell configurations 26,28,[45][46][47] . A red arrow in Figure 8 highlights this effect.

All total DOSs exhibit specific structures in the range between -1 eV and 0 eV. Given the range of energy and the amplitude of these peaks, we expect similar features to be detectable by Scanning Tunneling Microscopy analysis, and possibly used as easy measures to characterise the stacking in future experiments. However, these features are just vaguely related to the DOS of the deformed BP sheet, indicating that they cannot be put straightforwardly in relation to the deformation of the BP sheet. Indeed, this indicates that the hBN layer does actually influence the electronic properties of the phosphorus layer beyond the simple application of some stress.

P blue band structure FIG. 9. Band structure of the isolated P blue monolayer biaxially strained.

B. The P blue /hBN bilayer

We now focus on the electronic properties and the LP-DOS of the P blue /hBN layer. As reported previously, the BN is slightly deformed (weak compression of 0.4%) and does not lead to any particular electronic modifications. Meanwhile, most of the deformation is sustained by the P blue layer by a stretching of 1.25% of the in-plane parameter and a reduction of the buckling height. If one refers to the band plot of the isolated P blue monolayer under deformation (Figure 9), one realises that stretching the P blue layer of an amount of 1.25% leads to a sizable flattening of the conduction band in the Γ -M direction (blue arrow), and much weaker effects along Γ -K. This may have some impact in enhancing directional transport in this material because of the drastic increase of the electron effective mass along the Γ -M axis. This may also have some notable impact onto the optical properties of the layer. In particular, upon formation of electron-hole pairs, the flattening of the conduction band would result in a strong localization of the electronic part of the excitonic wavefunction.

Considering the valence band, the LP-DOS reported in Figure 10 shows an unexpected and striking feature. While the bottom conduction is completely localized on the P blue layer, the top valence presents a non-negligible hybridisation between the two layers. From an excitedstate perspective, this suggests that electron transport would take place only through the P blue layer, while holes could hop from one layer to the other. Similarly, excitons would have a mixed character, being well localised in the P blue layer for the electronic part of their wavefunction (especially if the P blue is stretched, as discussed above), while presenting some charge-transfer characteristic in the hole part of it. This point is of particular interest because the localization of the electron and the hole into two different layers (interlayer exciton) suppresses most recombination channels, hence stabilizing the exciton even at room temperature leading to promising applications 7,8 .

V. CONCLUSIONS

We have studied the structural and electronic properties of BP/hBN and P blue /hBN bilayers by means of DFT calculations with the intent of evaluating the impact of hBN on the phosphorus electronic properties in realistically relaxed structures.

In the BP/hBN systems, we find that there is a systematic repulsion between P and N atoms. Because of this, structures where the zigzag directions of BP and hBN are aligned (we called them C-cell structures) are not preferential configurations for the bilayer, at variance with those predicted in BP/graphene heterostructures 31 . Shearing the BP to reduce strain, has eventually a very high energetic cost and the resulting structure has even worse stability. Instead, by misaligning the two layers (Tcell) through a twist angle it is possible to obtain configurations slightly more stable than the C-cell or, at worse, of comparable stability, with very low deformations of the BP layer.

The compressive strain along the armchair axis required to fit in the C-cell is expected to spoil the transport properties of BP because of a competition between the top valence in Γ and a secondary valence band raising in between Γ -X. Instead, upon a twist angle of 19 • , the BP layer mostly stretches along the zigzag axis leaving the band structure on the armchair axis almost unaltered. As a result, the very promising charge mobility of the BP layer is not spoiled, and may actually be enhanced upon stretches of about 6%.

We also computed the stability, the strain distribution and the electronic structure of a P blue /hBN bilayer. Also in this system, the phosphorus layer deforms more than the hBN (biaxial dilatation of 1.25%), even though the deformation energy is shared more equally between the two layers as a consequence of the higher rigidity of the P blue with respect of that of the BP. We observed that, at low compression, a flat band develops in the bottom con-duction and that there is a non-negligible hybridization between the valence states of the two layers. These observations imply that (i) the electron mobility of P blue can be suppressed and the electronic part of excitonic wavefunctions can be localised on the P blue layer by relatively low strain; and (ii) that, upon creation of electron-hole pairs, holes can hop from one layer to the other. This work received funding from the European Union's Horizon 2020 research and innovation programme under grant agreements No. 785219 (Graphene Core 2) and from the French national research agency (ANR) under the grant agreement No. ANR-17-CE24-0023-01 (EPOS-BP).

  FIG. 1. a) Top and side view of the relaxed BP monolayer. Blue arrows mark the armchair and the zigzag axis. b) Top view of the hBN monolayer. The non-unitary cell of hBN is drawn with a dashed line. c) Top and side view of P blue monolayer. All structures are visualized with the VESTA package 33 .

  . The C-cell is constructed in such a way that the b axes of the two orthorhombic cells of BP and hBN are equal and parallel (parallel zigzag directions). It contains 3 BP cells and 8 B-N pairs (4 orthorhombic cells), for a total of 28 atoms. At equilibrium, the inplane a and b parameters are 10.01 and 4.36 Å respectively, the buckling parameter c of BP is of 2.14 Å. The

FIG. 2 .

 2 FIG. 2. Representative cells of the heterobilayers (thin dashed contours). The cells of constituent layers are drawn with thick continuous (unitary) or dashed (non-unitary) lines, with the same color-code of Figure 1.

FIG. 3 .

 3 FIG. 3. Cartoon of the relaxation procedure adopted in this work

FIG. 4 .

 4 FIG. 4. Stability curve of the monolayers. Black Circles: BP strained along the armchair axis; Black squares: BP strained along the zigzag axis; Black upward-pointing triangles: sheared BP biaxially strained; Blue downward-pointing triangles: P blue biaxially strained; Red diamonds: hBN biaxially strained. Dashed colored lines: fit of 1.

FIG. 5 .

 5 FIG. 5. Phosphorus/hBN cells: Total energy as a function of the deviation from the equilibrium surface ∆S = S -S0. Black squares: C-cell; Red triangles: S-cell; Green circles: Tcell, Blue diamonds: H-cell. Dashed colored lines: fit of 1.The red bullet at the origin marks the equilibrium structure of each system.

FIG. 6 .

 6 FIG. 6. Schematic representation of the steps from isolated monolayers to the heterostructure. Colored arrows represent the deformation applied.

FIG. 7 .

 7 FIG. 7. Band structure of the isolated BP monolayer. a) Strain ranging from -6% to +6% along the armchair axis; b) same strain applied on the zigzag axis. Blue and yellowish arrows highlight notable changes (see text).

FIG. 10 .

 10 FIG.10. Density of states of the isolated hBN and P blue monolayers (top panel) and the H-cell. In the bilayer (bottom panel), the P blue and the hBN curves refer to the monolayer contributions to the total DOS. In red solid line, the deformed P blue monolayer contribution. The top valence of all materials are aligned at 0. All data have been broadened with a Gaussian of σ = 0.05 eV.

TABLE I .

 I Cell parameters of the equilibrium P/BN bilayers.

	† The equilibrium parameters of the T-cell are derived from
	a fit, not from a full relaxation (see text).

TABLE II

 II

		B	Strain hBN (%)	Strain phosphorus (%)
	Cell type	meV/ Å2	armch.	zigzag	armch.	zigzag	buckling
	C-cell	25.36	0.29	-0.16	-5.32	1.46	0.77
	S-cell	25.22	0.02 *	-0.15 *	-1.83	0.82	0.50
	T-cell	26.73	-0.14 *	-0.14 *	-0.20	0.71	0.16
	H-cell	27.94	-0.40	1.28		-2.18

. Equilibrium surface S0 (in Å2 ) and 2D compressibility B (in eV/ Å2 ) of the monolayers resulting from the fit of expression 1.

TABLE III .

 III

TABLE V .

 V PBE energy gap in the isolated monolayers and the four bilayers.

		relaxed monolayers		P/BN bilayers
		hBN	BP	P blue	C-cell	S-cell	T-cell	H-cell
	Eg (eV)	4.71	0.81	1.93	0.83	0.86	0.91	1.81