Jean Galy
email: galy@lcts.u-bordeaux.fr

Gérard L Vignoles

The role of P 3s 2 lone pair (E) in structure, properties and phase transitions of black phosphorus. Stereochemistry and ab initio topology analyses

Keywords:

An approach merging crystal chemistry and density functional theory (DFT) electron localization function (ELF) taking P 3s 2 lone pair (E) into account induces a full renewal of stereochemistry of black phosphorus, its crystal network evolutions and phase transitions under increasing pressures from atmospheric up to 32GPa. Orthorhombic (Cmce) black P at ambient pressure, shows a packing of puckered [P] n layers -orthogonal to [010] -separated by a large free interspace (3.071Å), which actually is partially filled by lone pairs (E) (P-E ~ 0.8Å). Each P exhibits its lone pair pointing outside the [P] n layer, sandwiching it between two [E] n layers into a new stacking sequence …[EP 2 E] n … denoted O-[PE] n . The free interspace between [EP 2 E] n layers is much smaller 1.858Å but allows sliding along [001].

The pressure evolving up to 2.66 GPa, all structural details have been followed and reported, including the layer thickness reduction along [010] and the sliding along [001] of consecutive layers. A mechanism for the phase transition occurring around 5.5GPa is proposed. Depicted in the trigonal system the new layered phase R-[PE] n involves a bond rearrangement through E-E layer in zigzag phosphorus layers and P-E rotation and alignment with the A -3 axis. Now, the phosphorene layers have P-E patterns oriented towards each other in their interspace. A very particular phenomenon occurs around ~11 GPa the lone pair centroid Ec (P-Ec = 0.73Å) splits into three partially occupied sites Ed around the A -3 axis which explains observed variations in properties at this critical pressure. So, we claim that there are two trigonal phases, R1-[PE] up to 11 GPa followed by a second form R2-[PE] directly caused by lone pair displacement from Ec to Ed and its influence on layer stacking. A further layer sliding brings the phosphorus atomic layers close enough to each other to establish new P-P bonds and then to cause an ultimate transition to cubic system, with a new structure, isostructural to Po. The mechanisms of the transitions are detailed.

Introduction

Black phosphorus is at the center of active research owing to its particular layer structure made of phosphorene layers, partly analogous to graphene ones. The stable solid form crystallizes in the orthorhombic system (space group Cmce) as reported by Hultgren et al. in 1935 [1], here designed by O-P (A17 in literature). Brown & Rundqvist published refined structural data in 1965 [2] confirming its original layered packing. A little bit earlier, in 1963, Jamieson reported that under ~5 GPa O-P develops a phase transition to the rhombohedral (trigonal) system exhibiting a distinct layered structure after a rearrangement of atoms, R-P (A7) [3,4]. Applying higher pressures, ~10 GPa [4] or ~11 GPa [5], a new transition occurs giving a simple cubic form C-P isostructural with polonium Po. Important to note are the large electrical conductivity changes with phase transitions from semiconductor O-P [1,3] to semi-metallic R-P and finally metallic C-P [3]. This has been confirmed in a study by Kikegawa et al who locate this transition at 5.5GPa; they also report the transformation of semi-metallic R-P into a metallic cubic phase C-P occuring at ~ 11 GPa [3,6]. More recently a systematic study has been made by Scelta et al. using synchrotron radiation (ESRF) at up to 30GPa [7]; it shows the evolution of R-P structure up to 27.6GPa indicating the presence of a pseudo cubic sub-lattice in its X-ray data; finally an ultimate transition to the cubic form is reported at pressures >~30 GPa. A study at low temperature and high pressure allowed Shirotani et al. [4] to show that the minimal pressures to produce the phase transitions were higher at 21K than those at 77K. Li et al. [START_REF] Li | Proc. Nat. Acad. Sci[END_REF] studied resistivity, magnetoresistivity and Hall resistivity at low temperatures and high pressures, accompanying the phase transitions and confirming the onset of superconductivity at very low temperatures (7K or below) in black phosphorus single crystal. Our deep interest for lone pair stereochemistry covering a large panel of atoms, a systematic work supported by Density Functional Theory (DFT) Electron Localization Function (ELF) [START_REF] Matar | [END_REF][10][11], pushed us to reinvestigate this P in these O-P to R-P and finally C-P allotropes, showing interesting behavior element at high pressure and room temperature It has been shown in previous papers [START_REF] Matar | [END_REF][10][11] that lone pairs, designed by E, possess some pertinent structural data attached to each ns 2 element as demonstrated by Electron Localization Function (ELF) maps from DFT computations, briefly recalled hereafter.

-Three dimensional electronic density data show an excess concentration of electronic density in a small volume, close to some atom M*, having an intensity maximum, the so called centroid Ec of a lone pair, the coordinates of which can be determined; -the atom-lone pair distance M*-Ec can be calculated; -Ec is accompanied by a larger electron cloud E appreciated by its volume, roughly an ellipsoid, which corresponds to a sphere of influence (radius rE) a volume often compared to the one of an oxygen or a fluorine atom [12]. Its centre is often slightly displaced from Ec, the electronic cloud being sensitive to its environment and showing a certain degree of plasticity; -in all circumstances, E also shows a resistance to compression. The purpose of the present article is to bring new elements of understanding about the structure of the black phosphorus allotropes and to the transitions between them under pressure increase at ambient temperature, based on the identification of lone pairs, their positions, their steric occupation, and their reactivity upon pressure underlying the phase transformations, for which we propose network rebuilding mechanisms. In the following part, we will briefly recall our investigation methods; the rest of this document features five more parts, devoted to O-P, R-P and C-P, respectively, and in between to the transitions from one structure to another.

Methodology. X-ray data

We start from well-known and established structural data for phosphorus allotropes. Black P structures in the first pressure domain up to 2,66 GPa were taken from the results of investigations by Cartz et al. [13], in the orthorhombic system and Cmce space group corresponding to the O-P phase. For higher pressures (6,05 up to 27,55 GPa) corresponding to rhombohedral R-P structures, we have used the data of Jamieson [3], Kikegawa et al. [6] and Scelta et al. [7]. Finally, cubic P crystallographic data from Jamieson [3] and Kikegawa et al. [6] have been used. Orthorhombic gallium [14] and cubic polonium [5] structures are also given for comparison. All these data are gathered in Table 1 of the Appendix.

Computation framework

Computations were carried out by quantum density functional theory (DFT) with the Vienna Ab initio Simulation Package (VASP), version 4.6 [15][16][17][18]. A detailed presentation and use of such ab initio treatment has been given by Matar et al. in preceding papers [START_REF] Matar | [END_REF][10][11].

Our main purpose is to analyze the electron localization function (ELF) around phosphorus. The ELF was an original idea from Becke et Edgecombe [19] then extended and refined in 1991 and 1992 by Savin et al [20][21]. All structures were studied by single-point calculations with the atomic and lattice parameters described in the aforementioned references, using the Generalized Gradient Approximation (GGA) density functional [22] with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation potential corrected for solids [23], and the Projector Augmented Wave (PAW) method of Blöchl [24]. The plane-wave cutoff was set to 270 eV. The k-point spacing was set to 0.1 A -1 . The integration scheme used was Methfessel-Paxton with a 0.2 eV smearing, using smearing functions of order 1. The SCF convergence criterion was set to 10 -5 eV. Pre-and post-processing were achieved with the Medea software, version 2.12 [17][18] and some home-made tools developed on purpose by G. Couégnat at LCTS [10]. ELF normalized between 0 (zero localization) and 1 (strong localization) -with the value of ½ corresponding to free electron gas behavior -enables analyzing contour plots following a color code: blue zones for zero localization, red zones for full localization and green zone for ELF = ½. ELF slice representations are useful in such cases for quantitative discussion (metrics) of the lone pair development and stereo-activity. Some pertinent cross-sections are reported hereafter. Calculated ELF data allows depicting phosphorus layer inter-space via Isodensity Curves Analyses. Finally, the lone pair centers were localized in the ELF slices with the help of contour plots and their coordinates were measured and reported in the following tables. All distance values other than experimental data are reported with an estimated precision of 0.005Å. It is known that using dispersion-corrected potentials and modifications of the GGA method allow better computations of the interlayer forces and to discuss the relative stability of the allotropes [25,26]. However, in the purpose of the present work -which is only to localize the lone pairs -such corrections would only bring minute variations to the ELF density extrema positions. It is also well known that Wannier functions are another excellent tool for the determination of electronic density localization analysis [27]. In this work we have chosen not to use them, based on the experience acquired in previous studies, that have shown accurate enough to support our conclusions. Since the purpose of this work is exclusively dedicated to the lone pair localization, no attempt to discuss the energetics and bond strengths is made here. Likewise, we have not attempted to recompute the pressure from the DFT computations; however, as will be seen later, our ELF analysis results give consistent structure changes at the experimentally determined transition pressures.

Orthorhombic black phosphorus and its evolution at moderate pressures. 2.1. Classical description of the black phosphorus crystal structure

A refinement of black phosphorus crystal structure was reported by Brown and Rundqvist in 1965 [2]. This O-P variety crystallizes in the orthorhombic system with a base centered space group Cmce. X-ray data are reported in Table 1 of Appendix. The P atom x,y,z coordinates occupy the peculiar 8f (0,y,z) Wyckoff position which corresponds to a mirror plane perpendicular to [100]. In Fig. 1 a perspective view of O-P structure is given according to X-ray investigations. The drawing respects exactly the O-P network like reported by preceding authors [1-3, 6, 11, 12]. Discussions on their structural evolution and electrical properties are directly based on this restricted structure representation. All geometrical parameters vary smoothly with increasing pressure. ] n puckered layers are separated. The shortest P-P2' interatomic distance between layers amounts to some 3.592Å (dotted red line) which is also given in Fig. 4 by P-P4,5,6. Each P atom (black circles) is bonded to three other ones (P1,P2,P3) making an isosceles triangular pyramid.

A double row of P atoms along [100] draws a jagged ribbon, see for example …P3-P-P1…, in Figure 1 -which, repeated along [001] and bonded by P-P2 bonds, form a puckered layer [P] n . These layers, generated with b/2 periodicity, appear largely separated (3.108Å), with the shortest distance P-P2' = 3.592Å (dotted red line) giving this open layer structure. The thickness of the [P 4 ] n layer is given by ∆ intraspace = 2.131Å and the layer separation by ∆ interspace = 3.108Å. Bonds, interatomic distances and pertinent angles are reported in Tables 1 and2 (Appendix).

Comparison with Gallium

The Ga metal network shows strong similarities: same space group Cmce, same type of jagged ribbons [14]. However an important difference occurs: three successive jagged ribbons along [001] are now connected making again puckered layer interconnected by Ga-Ga3 = 2.442Å bonds which are the shortest ones (Table 2). So the b cell parameter for this three dimensional network becomes drastically smaller than the black O-P one. The Ga a and c cell parameters are larger than the O-P ones, sounds consistent with the Ga element being clearly larger than P. On the other hand, the Ga b cell parameter is markedly smaller: 7.660 Å against 10.478Å. This much larger b parameter can be justified by the hypothesis that the inter-space between phosphorus layers could be occupied by P 3s 2 lone pairs. This will be done in the next section.

P 3s 2 lone pairs in the crystal network of black phosphorus O-P

DFT-ELF calculations allow taking phosphorus 3s

2 lone pairs into account via their precise localization. Selected slices were extracted from the three-dimensional ELF data so created (Fig. 2 2 lone pairs E (pale blue circles). E-E inter-distances are marked by red sticks drawing a puckered E layer analogous to the P one (Table 2). b) The phosphorus tetrahedron PP1P2P3 and its E's show that each P in the structure network is inside a tetrahedron EP1P2P3. (G is the barycenter of a P1P2P3 triangle, a help to make path G-P-Ec and P13-P-Ec). A section by a plane perpendicular to P-E also allows 3s

2 pair ellipsoid size determination with its a/b/c parameters). a short P…P distance of P rows limiting a zigzag ribbon 1.483Å wide building one face, the opposite face being open (a parameter wide), the P rows being distant from the P'' ones by 2.894Å. Such opening is due to the presence of the ribbon of lone pairs which maintains this aperture with E"…E rows 1.702Å distance and P-Ec = 0.85Å. Thicknesses of successive layers are reported in Table 2.

A drawing of the core of that atom sequence between two successive layers directly responsible of O-P evolution versus increasing pressure is given in Fig. 4. Lone pairs localized between two [P] n layers show similar layers [E] n forming a triangular antiprism (elongated octahedron along [010]) E2E4E5E8E6E7 in this ∆P-P inter-space. Important to note that if both opposite isosceles triangles E2E4E5 and E6E7E8 show equivalent bases E4-E5 = E6-E7 = a = 3.313Å the former shows edges 2.375Å while the second has clearly larger ones equal to 3.144Å.

Structural evolution under pressure from atmospheric up to 2660 MPa

Taking into account the important presence of lone pairs the new real picture of black phosphorus leads to the identification of a fully renewed mechanism explaining its structure evolution versus increasing pressures. In this paragraph we focus onto the pressure domain from atmospheric to 2660 MPa taking advantage of the study by Cartz et al [13]. In this limited domain O-P black phosphorus evolves while keeping its orthorhombic structure (space group Cmce) (Fig. 3a). DFT-ELF ab initio calculations have been performed on this O-P form and we report their analyses based on the experiments at 0.1 (RT&P), 1510, 1730 and 2660 MPa [13]. This suffices, owing to the regular evolution of cell parameters and ELF structural data linked to lone pairs. The lengths of the P-P bonds (P-P1,3 and P-P2) as well as the P1,3PP2 and P1,3PP2 angles, the P-E and E-E distances and the thicknesses of various layers show a smooth but continuous evolution in this range of pressures. All data are reported in Tables 1 and2 (see Appendix). In this experiment the thickness of phosphorene [EP 2 E] n layer diminishes from 3.378Å at atmospheric pressure to 3.344Å, i.e. -0.034Å, while ∆E-E is condensed by 0.132Å, i.e. almost four times more. So it is reasonable to predict a shear slip movement in the (010) plane right in the middle of ∆E-E, the gliding of [EP 2 E] n layers favoring a rebuilding of a new network.

It will occur at a higher pressure, around 5.5GPa as shown by Kikegawa et al. [6] via a sharp drop of reduced cell volume at this pressure.

Phase transition towards As form at ~ 5.5GPa: Generation of R-P from O-P.

The schematic representation of this O-P → R-P transition is given in Fig. 5 after the abrupt drop of reduced volume vs. pressure at 5.5GPa [6]. Highlighted in green we have drawn a precursor pattern of the zigzag phosphorus layers obtained after the transformation. In it, P2, P6 and P7 atoms are still separated by lone pairs E2, E6 and E7. When the transformation occurs, a local electronic rearrangement occurs, with an exchange between lone pairs and bonding pairs. Thus, with minute atomic displacements, the new R-P structure is obtained. 6a). The double layer of [E] n is similar to the [P] n layers although it is considerably flattened (see E-E distances marked by red sticks). The P 4 pyramid, a basic element of the [P] n polymer, shows a perfect P1P2P3 equilateral base (3.377Å) with shorter edges P-P1,2,3 = 2.127Å. The lone pair E lying on top of P forms around it a regular tetrahedron E-P1,2,3 ~ 2.7Å. However E is also close to P4P5P6 of the next [P] n layer (E...P = 2.329Å), this equilateral triangle constitutes with P1P2P3 an elongated trigonal antiprism (or distorted octahedron -see dotted blue lines P-P = 2.854Å) (Fig. 6b). Ec is also equidistant of E4, E5, E6 with short E-E distances of 2.004Å (see dotted red lines Fig. 6b). The lone pair network of a [P] n layer does not interpenetrate neighboring ones: a little free space, around 0.464Å high, also remains between E's. It will be called ∆E-E intra-space, whereas the layer containing P-E will be called ∆P-E intra-space, the [P] n layer ∆P-P intraspace, the double layer of E's corresponds to ∆P-(E-E)-P inter-space. These layers thicknesses are precise indicators of the pressure effect (Fig. 6c). From 6,05 up to 27,55 GPa the layer thicknesses were followed systematically (see Table 3 in Appendix). It is worth to note that P-P (P-P1,2,3) bond lengths show limited variations, evolving from 2.127Å up to 2.299Å (after the jump from atmospheric to 6,05 GPa pressure). To visualize the detail of the network situation of R-PE four ELF sections calculated for 9,72 GPa are given as example in Fig. 7a,b,c,d. 3; we have P-Ec = 0.75Å. Ec centroid is indicated by a little red circle. c) Full E layer (E6E4E5) parallel to (001) plane. In E6E4E5 triangle center (z = 0.676) appears at section of Ec which is slightly underneath (z = 0.657). Note the continuous electronic pathway between the lone pairs. d) Section EcE4E5 shows the small separation of the electronic clouds favoring conductivity between consecutive E layers.

The R2-[PE] phase:

-Lone pair splitting around 11 GPa At 11 GPa, an accident in the reduced volume vs. pressure curve shows up , with rapid slope changes. This event marks that Ec starts to be brushed away from its 6c Wyckoff position (Wy) on the A -3 threefold axis towards a three more "spacious" ones around a -3 inversion axis as revealed by ELF sections. They show by symmetry almost identical z coordinates but also between x and y a relation close y = 2x. We could retainfor the E localization the 18f Wy x,y,z of R -3 (SG 148) but the sequence x,2x,z rather suggest the 18h Wy x,2x,z of R -3 (SG 166) which looks more appropriate. The triplets of sites are named Ed to mark this delocalization. A reasonable P-Ed ~ 0.85Å is maintained and of course these new sites show a partial occupancy, the 3s 2 doublet being able to "jump" from one site to another. Ec occupancy which was 100% for 6c Wy in R -3m (SG 166) becomes 33.3% now distributed onto its 18h Wy sites (Fig. 8a,b,c).

The peculiar network evolution of the R-[PE] trigonal phase.

The evolution of the lone pair site (Ec to Ed) seems to have drastically disrupted the black P network. Accordingly, it is now chosen to follow the bonding system evolution through the joint behavior of the P-P1 strong bond of the phosphorene layer and the P-P distance between two consecutive layers, i.e. P-P6 vs. increasing pressure between 6.05 and 27.55 GPa. There is a covalent bond between P and P1 atoms with an interatomic distance 2.291Å at 6.05 GPa whereas P6 which belongs to the neighboring phosphorene layer is clearly not covalently bonded to P with P-P6 = 2.620Å. The values of these interatomic distances reported in Table 3 of the Appendix are plotted in Fig. 9. As could be expected their evolution is drastically modified at ~11GPa with smoother slopes after Ec site delocalization in three Ed. The transition is marked by the fact that the P-P1,2,3 distance recovers the value at 6.05 GPa; moreover, ∆P-P (intra) starts to increase from 1.212Å up to 1.287Å and P-Ed shows an expansion from 0.827Å up to 0.845Å. For pressures larger than 11 GPa, with a limit of 27.55 GPa, P-P1 shows a slow linear decrease (blue line)2 while P-P6 decreases, but more markedly and following a third order polynomial3 (magenta curve) down to 2.335Å at 27.55 GPa.

At that pressure, P-P6 is now close to the P-P1 bond length and a strong bond is being established. To illustrate this, six ELF sections calculated on the basis of the three atoms P, P1 and P6 enlighten this P-P bond generation as shown in Fig. A1a-f (Appendix) for the pressures indicated in Table 3. In Fig. A1a,b,c, before the transition, the P-P1 bond shows a marked electron density concentration in its middle, while P-P6 remains as a simple electronic pathway. Then, after the transition, Fig. A1d signals the beginning of electron density concentration right in the middle of P-P6, an evolution which continues up to 27.55 GPa. Through all these modifications the role of lone pairs in Ec or Ed is of paramount importance in the inter-space of P layers and they master also the condensation and stabilization of the phosphorene mono-layer [EP 2 E] n , the thickness of which diminishes from 2.416 to 2.323Å. Such condensation is due to Ed which drives opposite variations of ∆P-P, ∆Ed-Ed(intra) from 0.367Å (after 0.355Å) to 0.348Å and ∆(P-Ed) from 0.602Å to 0.518Å. Important also is to note that [P] n layer thickness increases (1.212Å up to 1.287Å), a fact which could be surprising when almost all lengths are diminishing under pressure increase. This is readily explained by considering the tetrahedron PP1P2P3 with its equilateral base P1P2P3. This base remains equilateral though the edge lengths diminish with pressure, whereas P-P1,2,3 bonds do not change appreciably; this results in a clear closing of the zigzag angle in the [P] n plane with a value decreasing from 99.1° at 6.05 GPa to 91.8° at 27.55 GPa.

This triggers a lengthening of the PP1P2P3 tetrahedron height and then of the ∆P-P layer thickness. A more precise insight into the R1 → R2 transition is available in Fig. 11, with ELF slices taken at 10.775 and 16.25 GPa. In Fig. 10a E is not anymore marked by a quasi spherical high electronic density but by a kind of crescent-shaped region with somewhat looser density.

After calculation of Ec geometric coordinates on A -3 section PP1E shows the particular deformation of the Ec electronic cloud (pale blue ellipse, Fig. 10b). In Fig. 11c ELF data reveal electronic concentration around A -3 but at a z coordinate closer to P than to Ec. Subsequently, recognizing that the Ec site has been abandoned, the site will be called Ed. The coordinates of such Ed sites were refined directly from ELF data; an equatorial section of this Ed delocalized site shows clearly the three crystallographic sites (18h) (Fig. 10c), of course partially occupied. They form a locally maximal electronic density region constituting a continuous electronic closed loop. The planar section EcEdbEdc shows that an electronic flux is maintained between the Ec site and the new delocalized ones Eda,b,c favoring a continuous pathway for electric current between E planes (Fig. 10d). When pressure further increases, at 16.52 GPa, Fig. 10e shows that Ec (6c) site marked by "E?" has disappeared and that the new structure takes place with the lone pair distributed only in Ed sites. Finally, refined Ed x,y,z coordinates around A -3 and the subsequent section plane clearly evidence the three crystallographic sites (18h), of course partially occupied, in a two-dimensional electronic density loop (Fig. 10f).

Taking advantage of these last ELF calculations for R2-[PE] at 27550 MPa a section into median plane of ∆Ed-Ed, right in the smallest interspace (0.348Å), is illustrated in Fig. 11 together with the Edabc plane. Going back to the graph in Fig. 9, we see that extrapolating the curves of P-P1 and P-P6 bond lengths (dotted blue line), there would be an intersection at ~34GPa with a common bond length value of 2.28Å, indicating that phase transition R2-[PE] → C-[PE] takes place there, as discussed in the next section.

C-[PE] cubic form (isostructural with Po) 5.1. R2-[PE] → C-[PE] phase transition mechanism

The phase transition from the layered R2-[PE] structure, to the three dimensional cubic C-[PE] is latent in the network architecture of the former. Upon increasing pressure there is a slip of the [P] layers and a flattening of the double E layers. The zigzag layer ... P4P7P6P8... slips toward the ...P2PP0P1P9...one formed by P-P bonds (black sticks), these two layers being still separated by the double layer of split lone pairs Ed with P-P distances marked by red sticks as illustrated in Fig. 12a,b. At 16.52GPa, looking at the pseudo cubic polyhedron P1P0P4PP6P7P8P9, P-P1 bonds (black sticks) equal 2.300Å, while the P-P6 ones (in red) across Ed bi-layer are longer 2.407Å. Similar to the two P1P0P and P7P6P8 black triangles of each phosphorus layer there are two red triangles P4P0P and P9P6P8; all triangles together constitute a deformed cubic-like polyhedron. At 27.55GPa, P-P1 and P-P6 distances are 2.279Å and 2.335Å, respectively: they are now very similar. In the distorted cube the dihedral angles P1(P0-P) P4 = P7(P6-P8)P9 = 178.1° due to little difference between edges P1-P9 = P4-P7 = 2.279Å and P-P6 = P0-P8 = 2.335Å as seen in Fig. 12c,d

The C-[PE] phase

According to the data of Scelta et al. [7] the R2-[PE] cell becomes rather close to the cubic cell obtained by Kikegawa et al. [6] at 32GPa. The phase transition R2-[PE] to C-[PE] occurs at pressures in between 28 and 32GPa and lone pair sites reorganize themselves jointly to become isostructural with the Po structure [5]. An ELF section by (001) plane is displayed in Fig. 13a and shows a channel of electron density running along cube edges. On edges, in between P atoms there is an island of maximum density with two peaks. These peaks occupy two Wyckoff positions 6e (x,0,0 and 1-x,0,0) with 16.7% occupancy. The x coordinate in these 6e positions is found to be 0.392. It comes that P-E = 0.89Å and the separation of two consecutive sites in x,0,0 and 1-x,0,0 is E-E = 0,493Å (Fig. 13c). These data show that the memory around P coordination derives well from R2-[PE] in spite of the phase transition. There is no electronic density in the cell centered on cube faces neither in its center as shown by the ELF section of Fig. 13b. Every P atom is surrounded by six E sites making an octahedron, only one E site being occupied by edges.

About pseudo-cubic phase presence in R2-[PE] according to X-Ray data

Scelta et al [7] report in their X-ray study under increased pressures the presence of a pseudocubic phase (pC-[PE]). Recalling our description of the phase transition R2-[PE] to C-[PE] we can make a proposal concerning this pC-[PE] phase. The calculated powder XRD patterns of R1-[PE] (blue) at 6.05 GPa, R2-[PE] (red) at 27.55 GPa and C-[PE] (black) one at 32 GPa [6] are drawn in Fig. A2 in Appendix. Comparing R1 and R2 peaks we mostly observe a shift towards higher Bragg angles from R1 to R2 indicating a lower cell volume with increasing pressure. Both diagrams show an intense sublattice of reticular planes (012), ( 104)-( 110), ( 006)-( 202), ( 116), (214), … markedly more intense than any other peak. Very interestingly, these marked peaks are lying close to the C-[PE] peaks. For instance, (012), ( 110) and (006) of R2-[PE] almost match (100), ( 110) and (111) of C-[PE], respectively, as suggested by Fig. 15 (Appendix). This confirms that the phase change from R2 to C does not require any large atomic rearrangement.

Recapitulation of black P allotropes and phase transitions.

To get a global view of black phosphorus structure evolution at room temperature and various pressures in the range 0 up to 32 GPa a diagram of the reduced volume vs. pressure V(Å Data points are available for some pressure ranges: from 0.1 to 2.66 GPa [5] and to 5.5 GPa [13] for O-[PE], from 5.5 GPa to 9.7 GPa [13] and from 6.05 and 27.55 GPa [14] for R-[PE] and at 11 GPa [5] and 32 GPa [13] for C-[PE].

All rV vs. p data sets have been fitted with 1 st , 2 nd or 3 rd order polynomial equations 4 . The cell volume diminution from atmospheric to 9.7 GPa, according to Kikegawa et al. [13] (cyan color -curve IIa,b), allows locating the transition of O-P to R-P at 5.5 GPa. The rV of O-[PE] form has an initially steady linear decrease, whereas R-[PE] shows distinct evolutions in the interval of pressures between 6.05 and 10.775 GPa as compared to the rest of the graph, up to 27.55 GPa. Two different curves (III red and IV blue) were fitted. They intersect at 10775 MPa which is close to the admitted pressure of the R1 R2 transition. Maybe by coincidence, curve (I) also intersects the other two at the same pressure. In the R1 allotrope, as soon as pressure increases the rV decreases quickly, principally because the phosphorene layers approximate to each other: ∆(Ec-Ec) drops from 0.218Å to 0.091Å. After having aligned P and E along the threefold inversion axis, the Ec site splitting in favor of the three partial Ed sites leads to a very distinct pressure evolution, of the reduced volume. Concerning the C-[PE] cubic cell, there are only two rV [PE] reduced cell data values available, plotted in Fig. 16 as magenta triangles: one is 13.43Å 3 at 11.1 GPa [3], close to the transition pressure 10.775 GPa of R1-P to R2-P with rV = 13.60Å 3 , the second with rV = 11.7Å 3 at 32 GPa [6]. The two data points have been connected by a dotted line which evolves close to the blue one of R2-[PE] rV. This illustrates again the fact that over the whole range of R2-PE existence confusion can be made between its cubic sub lattice and the true cubic phase. Summarizing, we have described the structural evolution of black phosphorus as going through three phase changes, accompanied by changes of electronic conduction type: from semiconducting O-P, to semi-metallic R1-P, then to metallic R2-P and finally to (possibly) superconducting C-P.

Concluding notes

Orthorhombic black phosphorus exhibits in its crystalline network wide empty space between its puckered [P] n layers. Such peculiar puckered network architecture is also the case of Ga but in the latter case the layers are firmly associated by a Ga-Ga bond resulting in a much shorter parameter perpendicular of the layers. The difference resides in the existence of lone pairs in the case of phosphorus. These lone pairs play a major role in black phosphorus structure and transformations under high pressures. So new hypotheses implying phosphorus 3s 2 lone pairs were developed and supported by crystal chemical and DFT-ELF analyses in order to explain structure evolution versus increasing pressure and to propose a precise description of phase transition mechanisms.

• After DFT and ELF calculations the omnipresence of lone pairs has been established. In the orthorhombic phase O-[PE] the puckered [P] n layers appear sandwiched between E layers: these electronic doublets E make a double layer fully occupying the ∆P-P inter space while the ∆P-P intra space remains empty.

• As in our preceding papers E distinguishes itself by an important electronic density concentrated in a small volume, close to the P atom, accompanied by an electronic cloud exhibiting some plasticity. Ec coordinates in the network have been systematically determined from ELF analysis giving P-E ≈ 0.85Å. The E size appreciated as a sphere of influence shows a radius rE ≈ 0.90Å.

• Following O-[PE] network evolution when submitted to increasing pressures from atmospheric to 2660 MPa, c and b parameters decrease while a remains almost constant as well as P-P bonds (around 2.222Å). Actually E plays a major role by its ability to move, despite being always associated to P with P-E remaining quasi constant. The double E layer allows an easy sliding movement between them and also a compaction which has been described in detail.

• Based on the study of the lone pairs, a mechanism is proposed for the phase transition around 5 GPa conducing to a rhombohedral -trigonal structure characterized by the packing of zigzag [PE] n layers, the R-[PE] form. E continues its major influence associated with P, with P-E patterns aligned along A -3 axis. When pressure increases the P-E pattern undergoes a striking rearrangement around 11 GPa: E is expelled from A -3 and split on three new sites related by symmetry, keeping the link P-E but with a partial occupation of the new sites. We propose to distinguish two rhombohedral structures: R1-[PE] n up to 11 GPa and R2-[PE] n for pressures up to 27550 MPa.

• The compression of R2 structure leads to a transition towards a cubic structure by a very simple mechanism involving minute atom and electron density displacements. This transition can be anticipated by the fact that a cubic sublattice already exists in the R2 phase.

• Thanks to lone pair sliding, rearrangements and splitting into triplets, the shrinkage of the phosphorene layers submitted to increasing pressure leads to a reduced volume rV [PE] diminution from 19.0Å

3 at atmospheric pressure to 11.7Å 3 at 32 GPa (some 38% loss).

APPENDIX

Figure 1 -

 1 Figure 1 -Perspective view of black phosphorus orthorhombic network (O-P) at room temperature and pressure according to ref. 2. The zigzag sequence P1PP3 aligned along [100] defines ribbons alternatively linked by P-P2 bonds making the puckered layer [P 4 ] n parallel to the (010) plane. [P 4] n puckered layers are separated. The shortest P-P2' interatomic distance between layers amounts to some 3.592Å (dotted red line) which is also given in Fig.4by P-P4,5,6. Each P atom (black circles) is bonded to three other ones (P1,P2,P3) making an isosceles triangular pyramid.

  ). The first slices a), b) and c) correspond to cross-sections containing three P atoms, P-P1-P3, P-P13-Ec (P13 middle of P1-P3) and P2-P-Ec planes. The d) cross-section passes through three consecutive lone pairs Ec, E1 and E3. In slices a) and b) it is worth to notice between P-P1, P-P3 and P-P2 a peak of electronic density right in their middle corresponding to strong covalent P-P bonding. Besides, another important high-density zone on the top of P atoms appears corresponding to a section of the electronic cloud around P. A search for a maximum E density through isodensity curves was realized and the lone pair centroid Ec coordinates were extracted. Calculations based on sections via P, P2, P13 and Ec allowed delineating E as an ellipsoidal form covering P like a helmet (Figs.2b,c). Section d) passing through three consecutive E lone pairs shows their intense density and the electronic cloud zigzag along[100]. The values of P-Ec ~ 0.85Å and rE ~ 0.91Å are in rather good agreement with the ones obtained in investigations on phosphorus oxides and fluorides. Structural data are reported in Table2of the Appendix.

Figure 2 -

 2 Figure 2 -a) ELF section PP1P3 revealing electronic density in this ribbon face. P-P1 and P-P3 bonds are clearly marked by an electronic density in their middle. Also notice a section of lone pairs density covering P, P1 and P3 and the continuous electronic pathway (red color) connecting …E3, PP3 bond, E, PP1 bond, E1. b) P2PEc section which, combined with c) E planar section perpendicular to P2-Ec, allows measuring E ellipsoid parameters a/b/c. d) Section determined by three E3EcE2 lone pairs showing that lone pair electronic clouds may have an active participation to electric conduction together with the [P] n network. As a consequence, the black phosphorus O-P network can be revisited, as illustrated in Fig. 3a & b. The inter-space between puckered [P] n layers is now filled by a double layer of lone pairs making along [010] a periodic sequence ...∆P-E/∆E-E/∆E-P/∆P-P/∆P-E. Out of the [P] n

Figure 3 -

 3 Figure 3 -a) Black P network (O-P) with the "free" interlayer space of [P] n puckered layers packing filled by a double layer of P 3s2 lone pairs E (pale blue circles). E-E inter-distances are marked by red sticks drawing a puckered E layer analogous to the P one (Table2). b) The phosphorus tetrahedron PP1P2P3 and its E's show that each P in the structure network is inside a tetrahedron EP1P2P3. (G is the barycenter of a P1P2P3 triangle, a help to make path G-P-Ec and P13-P-Ec). A section by a plane perpendicular to P-E also allows 3s 2 pair ellipsoid size determination with its a/b/c parameters).

  c) "Exfoliated [EP 2 E] n layer" of black phosphorus with its [P] n puckered layer as framework (highlighted by grey color), sandwiched between two planes of lone pairs E.The green line in the middle of [E] n layers inter-space marks a cleavage plane perpendicular to [010] (green plane S) possibility to exfoliate the structure and to obtain phosphorene [EP 2 E] n single layers showing analogy with graphene (Fig.3c). It is worth to note the various thicknesses of the successive layers depicted in Fig.4. If no attention to E is paid, the intra-layer and inter-layer ∆P-P distances are 2.131Å and 2.973Å, respectively. Actually, the inter-layer value corresponds to the sum of three layer thicknesses ∆P-E + ∆E-E + ∆E-P = 0.587 + 1.799 + 0. 587 = 2.973Å (see Fig.3a,c). After cleaving the crystal into isolated phosphorene [PE] n layers, their thickness amounts to ∆E-P + ∆P-P + ∆P-E = 0.587 + 2.097 + 0.587 = 3.271Å. Before assessing the effect of pressure on the structure some other features must be underlined. There are two slightly different P-P bonds; in the ribbons along [100] P-P1,3 = 2.223Å while their connection along [010] is assumed by P-P2 = 2.245Å with an angle P1PP2 = 102.1°. The puckered [P] n layer shows, on both crenelated faces orthogonal to [010],

Figure 4 -

 4 Figure 4 -A scheme to follow O-[PE] black phosphorus at room temperature and various pressures. Black circles represent P phosphorus and blue circles lone pairs E with centroid marked by deep blue ones accompanied by its electronic cloud (pale blue circle).Thicknesses of successive layers are reported in Table2.

Figure 5 - 4 .

 54 Figure 5 -a) In the O-P network the green shaded area (normal to [021]) shows a succession of zigzag for …P-P-P-P files parallel to [100] separated by E -E sequence (dotted red lines); b) Such layer is extracted as shown in 5b; c) Then applying pressure (~5 GPa) a network mechanism takes place affecting lone pairs of P2P6P7 which are repulsed along blue arrows conducing to the O-P → R-P phase transition; 5d) Resulting new [PE] n layer with all P-E links parallel to [010] and showing that the zigzag P layer is sandwiched by a pair of E layers.

Figure 6 -

 6 Figure 6 -a) Zigzag P layers formed by P atoms (black circles) set up onto 6c Wyckoff sites of inversion -3 axes (SG 166 R -3 2/m) together with their Ec; note the face developed by …PP1P3… atoms in zigzag bond sequence P1-P and P-P3, while the other face shows parallel P-P2 bonds. The zigzag angle of the layer is P13PP2 where P13 is the middle of P1-P3. b) The basic unit: a PP1P2P3 tetrahedron with P at apex and the base P4P5P6 shared with the next unit along threefold axis create together a triangular antiprism (elongated octahedron) around P. Note also that P2P1P3 and P6P4P5 form perfect equilateral triangles. c) The R-[PE] crystal network with its stacking of successive P and E double layers for pressures <10GPa generating phosphorene layers separated by a double layer of lone pairs ∆E-E.

Figure 7 -

 7 Figure 7 -Four sections of R-[PE] ELF data at 9720 MPa showing details of lone pair 3s 2 Ec localization and topology in the network. a) ELF section by P3PP2 plane showing the presence of the lone pair E on top of P atoms. b) A section plane based on the P1PEc triangle, allows determining the P-Ec direction c parameter of lone pair ellipsoid formed by Ec centroid accompanied by an electron cloud. Ec coordinates are reported in Table3; we have P-Ec = 0.75Å. Ec centroid is indicated by a little red circle. c) Full E layer (E6E4E5) parallel to (001) plane. In E6E4E5 triangle center (z = 0.676) appears at section of Ec which is slightly underneath (z = 0.657). Note the continuous electronic pathway between the lone pairs. d) Section EcE4E5 shows the small separation of the electronic clouds favoring conductivity between consecutive E layers.

Figure 8 -

 8 Figure 8 -a) Schematic representation of R-[PE] black phosphorus right in the moment (domain around 10.8GPa) where lone pair starts to partly delocalize its presence on the previous one Ec onto A -3 axis and three Wy Ed sites around it. The little red torus links up these new sites marking a potential electronic pathway. Ec is shown by a little blue circle and Ed by a pale blue one. b) Section plane passing through Ed(a,b,c) showing the two dimensional electronic linkage between these Ed sites triplets. c) Section plane EcEdbEdc shows that electronic density is marked in between the Ec and Ed levels indicating the easiest electron delocalization between these crystallographic sites, explaining metallic character of this new phase called R2-[PE].

Figure 9 -

 9 Figure 9 -P-P1 bond (black & blue curves) and P-P6 interatomic distance (red and magenta curves)vs pressure plot for pressures between 6.05 and 32.0 GPa.In the [6.50-11.0] GPa range (the domain of R1-[PE] phase) second order polynomials fit the evolution of P-P1 and P-P6. Then, up to 27.55 GPa, the P-P1 evolves linearly while P-P6 follows a third order polynomial. Extrapolation of these two curves give an intersection close to the green diamond, which defines a C-[PE] phase at ~32GPa.

Figure 10 -

 10 Figure 10 -ELF sections evolving from 10775 MPa to 16520 MPa. a) PP1P3 section plane shows electronic density changes around Ec calculated position. b) The section PP1E, E position on A -3 calculated. c) ELF data reveal electronic at z very close to P level at Ed. They form a locally maximal electronic density region constituting a continuous electronic closed loop. d) Planar section EcEdbEdc shows an electronic flux between Ec site and Eda,b,c. e) At higher pressure 16.52 GPa the Ec site marked by "E?" has disappeared. f) Ed around A -3 is surrounded by a well delineated electronic closed loop.

Figure 11 -

 11 Figure 11 -The Edabc plane at 27550 MPa which shows Ed spots around threefold axis (left) at z = 0.3551 and the median plane (right) which shows a rich electronic pathway in between consecutive Ed layers.

  . The angles PP1P0 = P6P9P8 = 91.2° and P1PP4 = P9P6P7 = 90.1° are close to 90°.

Figure 12 -

 12 Figure 12 -a) and b) Under 27550 MPa black and red thick sticks delineate the polyhedron which generates the cubic cell of C-[PE]; c) the slightly distorted cubic cell with P and E; d) The plasticity of the lone pair cloud will favor their crystallographic site reorganization.

Figure 13 -

 13 Figure 13 -a) ELF section of (001) plane showing the electronic concentration along cube edges and particularly the density island with two peaks (6e) on P-P bond at P-E = 0.889Å; b) ELF section at z = ½ showing that cell cube faces and center are free of electronic density; trace of edges are clearly materialized; c) cell of crystal network of C-[PE] n . Potential sites for E going past are indicated by pale blue circles.

  3 ) = f (p) is given in Fig.14. rV is the reduced volume of [PE] obtained by the cell volume divided by the number Z of [PE] in the cell, rV = V/Z.

Figure 14 -

 14 Figure 14 -Evolution of black phosphorus reduced volumes (rV = V(cell)/Z([PE])) versus pressures up to 27550 MPa and the three transitions conducing from O-[PE] to R1-[PE] and R2-[PE] [6, 7, pw] then to C-[PE]. Two magenta triangles at 11GPa [13] and 32GPa [7] show rV(Å 3 ) of C-[PE] cubic form; the magenta triangle for p = 27550 MPa has been calculated using linear formula proposed for C-[PE] cubic evolution. It is seen that 11GPa is a critical point, opening on both R2-[PE] and C-[PE] final structure perhaps coexisting on largest 11-32GPa pressure scale. However the R2-[PE] sublattice appears as a good alternative to explain the presence of a pseudo-cubic structure in this domain.

Figure A1 Figure A1 -

 A1A1 Figure A1

Figure A2 Figure A2 -

 A2A2 Figure A2

Table 1 -

 1 X-ray cell crystal data at room temperature and high pressures of P varieties, orthorhombic (O-P), rhombohedral (R-P) and cubic (C-P).

	P (Po type) Cubic, Pm -3m (221) [3,5] noted C-P.
		a(Å)	b(Å)	c(Å)	V(Å 3 )	Z	rV[PE] (Å 3 ) Ref.
	O-P RTP	3.32	10.52	4.39	153.3	8	19.2	[1]
	O-P RTP	3.3136	10.478	4.3763	151.94	8	19.0	[2]
	O-P RT P = 0.1,	3.3133	10.473	4.374	151.8	8	18.98	[13]
	up to 2660 MPa	3.3120	10.140	4.229	142.0	8	17.75
	5500 MPa	3.31	9.97	3.95	130.4	8	16.3	[6]
	R-P RTP	3.377	3.377	8.806	86.97	6	14.50	[3]
	P = 9700 MPa	3.39	3.39	8.63	85.89	6	14.32	[6]
	R-P RT..P = 6050,	3.3997	3.3997	8.7514	87.60	6	14.60	[7]
	up to 27550 MPa	3.2576	3.2576	8.0103	73.62	6	12.27
	C-P RT P 11.1GPa	2.377	2.377	2.377	13.43	1	13.43	[3]
	C-P RT P 32 GPa	2.27	2.27	2.27	11.70	1	11.70	[6]
	Ga RTP	4.5198	7.6602	4.5258	156.69	8	19.6	[14]
	Po RTP	3.359	3.359	3.359	37.90	1	37.90	[5]

Black P: Orthorhombic, Space group Cmce (64) [1-3] noted O-P. P (As form): Rhombohedral, Space group R -3m H (166) [3,13-7] noted R-P.

P-P1 (Å) = 2.306 -9.097.10 -7 .p 2 P-P6 (Å) = 2.935 -6.548.10 -5 .p + 2.691.10 -9 .p 2 -4.002.10 -14 .p

(p expressed in MPa).

I rV = 18.958 -4.666xE -4 xp (p expressed in MPa) IIa rV = 18.962 -7.861xE -4 xp + 1.537xE -7 xp 2 -1.477xE -11 xp 3 IIb rV = 19.505 -0.00109xE -4 xp +

5.549xE -8 xp 2 III rV = 13.578 + 3.826xE -4 xp -3.532xE -8 xp 2 IV rV = 14.992 -1.5xE -4 xp + 1.871xE -9 xp 2 V rV = 14.308 -8.151xE -5 xp

Acknowledgements

J. Galy wishes to thank LCTS -CNRS for providing laboratory facilities. Thanks also to G. Couégnat (LCTS) for his efficient support in computing.

Table 2-Interatomic distances (Å) and angles (°) from XRD [13,2]

and ELF data introducing E [present work (pw)].

O-P [13,pw] Ga [14] Ga-Ga1 = 2.801 Ga-Ga2 = 2.742 Ga1-Ga2 = 2.712 Ga-Ga3 = 2.442 Ga1-Ga1a = 4.520 E coordinates have been determined and refined from three-dimensional ELF data. Inter layer distances between [P] n puckered layers are noted ∆P-P (Å) inter, ∆P-P or ∆Ε-E (Å) (intra) indicate the thicknesses of P or E layers, ∆P-E (Å) between P and E planes perpendicular to [010] (see Fig. 2).

Table 3-Interatomic distances (Å) and angles (°) from X-rays [7] and ELF data ([pw]).