Superlubricity in phosphorene identified by means of ab initio calculations

Phosphorene possesses a great potential for tribological applications due to its layered structure and for the capability of phosphorus to reduce friction and adhesion in steel–steel contacts. Here we present a comprehensive analysis of the static tribological properties of phosphorene based on first principles calculations. The most suitable exchange-correlation functional for describing the structural and electronic properties of multilayer phosphorene is carefully selected. The interlayer binding energy and shear strength are then calculated for two relative orientations of the layers. Layers stacked with the same orientation (armchair–armchair and zigzag–zigzag) are slippery as common solid lubricants, as MoS2 and graphite. While the armchair–zigzag orientation shows a remarkable superlubricity, with a reduction of one order of magnitude for the shear stress. We uncover the microscopic origin of such superlubric phase by analyzing the electronic charge at the layer interface.


Introduction
Black phosphorous and its single layer constituent, phosphorene, are emerging as promising materials for technological applications in opto-and electronics devices. Thanks to its direct band gap of around 0.3 eV [1,2] phosphorene band gap lays between the almost zero band gap of graphene and the relatively large (1.5-2.5 eV) gap of transition metal dichalcogenides (TMDs). More importantly, the band gap of phosphorene films can be tuned by changing the number of layers [3] and by applying an in-plane strain to the film [4]. This intriguing property allows to modify the nature of phosphorene conductivity, from a direct to an indirect band gap semiconductor, semimetal or metal.
The mechanical properties of phosphorene are also extremely interesting. Thanks to its peculiar puckered structure, phosphorene shows highly anisotropic Young's moduli, the zigzag direction being about twice times more stiff than the armchair one [5], a property rather uncommon in 2D and, more in general, 3D materials. This intrinsic anisotropy can favor the development of new phosphorene-based devices useful not only for electronic, but also for photonic and thermoelectric applications [6].
Likewise other 2D layered materials, such as graphene and TMDs, black phosphorous layers are weakly bonded through van der Waals (vdW) interactions. Therefore, phosphorene layers have the intrinsic property to easy shear between each other, leading to a great potential for tribological applications. Moreover, P-based lubricant additives are commonly employed in engine oils to reduce friction and adhesion in high pressure applications [7,8].
Here we investigate the tribological properties of phosphorene by first principles calculations. This approach allows for an accurate description of the adhesive contribution to friction and has been successfully applied to describe the static friction properties of solid interfaces [29][30][31] and the interlayer slipperiness of layered materials, such as graphene [32,33] and TMDs [34][35][36]. In particular, it predicted the superior performances of graphene than MoS 2 [37] as recently verified by nanotribological experiments [38].
An ab initio study of the interlayer slipperiness of black phosphorous properties is lacking to the best of our knowledge. We first analyze the structural and energetic properties of phosphorene films of different thickness, adopting five different exchange-correlation functionals. We start from the most common schemes, then we introduce vdW corrections, which are implemented in both semi-empirical and self-consistent ways. Then, we proceed calculating the adhesion energy, potential energy surface, and shear strength for two different bilayer orientations. Extremely low corrugation and shear strength are detected when the upper layer is rotated by 90° around the z axis with respect to the lower layer. This orientation is stable and gives rise to a structural superlubricity, whose origin is explained by analysing the electronic charge distribution within the interlayer region.

Method
We perform density functional theory (DFT) calculations using the Quantum Espresso package [39,40], in which the electronic wave functions are expanded in plane waves and the ionic species are described by ultrasoft pseudopotentials [41]. Two different approximations are considered for the exchange-correlation (XC) functional: the local density approximation (LDA) [42] and the generalized gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE) parameterization [43]. Dispersion interactions are included according to different vdW schemes implemented on top of PBE, namely PBE + D2 [44], vdW-DF2 [45], and rVV10 [46]. According to the PBE + D2 scheme, a semi-empirical correction term is added to the total energy by evaluating atomic pair interactions, while, in the vdW-DF2 and rVV10 schemes the nonlocal contribution is calculated ab initio. A scaling factor s6 equal to 0.75 is chosen in PBE + D2 as it optimizes the interlayer distance of bulk phosphorene.
After test calculations on the bulk structure, a kinetic energy cutoff of 30 Ry is adopted to truncate the plane wave expansion in all the calculations apart from the vdW-DF2, where a value of 35 Ry is used. The energy cutoff for the charge density is chosen to be either 8 times (LDA, PBE, PBE + D2) or 10 times (vdW-DF2, rVV10) the one for the kinetic energy of the wave functions. The Brillouin zone of the bulk cell is sampled by means of a 12 × 9 × 4 Monkhorst-Pack grid [47].
We first model phosphorene films of different thickness and analyze the variation of the structural and electronic properties as a function of the number of layers, N, that is increased from one to five. After a full relaxation of all the structures, including the lateral cell dimensions, the interlayer distance, d int , is obtained as the average distance between the innermost planes of each couple of layers contained in the supercell. The interlayer binding energy is calculated as where E N is the energy of the film composed by N layers, E 1 is the energy of a single relaxed layer (calculated within the same supercell), and A is the in-plane area of the supercell. In the case of a bilayer, the interlayer binding energy is also referred to as adhesion energy, γ.
As a second step in our analysis, we study the interlayer slipperiness of phosphorene layers by applying a Workflow, which has been recently developed by our group to perform high-throughput calculations of interfacial properties, including tribological figure of merits [48,49]. In particular, the potential energy surface (PES) of the sliding interface is constructed automatically by calculating the adhesion energy γ, which is the opposite of the work of separation [50], as a function of the relative lateral positions of the two surfaces in contact, in the present case represented by two layers. For each layer stacking, the (x, y ) coordinates of all the atoms are kept fixed, while their z coordinate is allowed to relax. The whole bidimensional PES is then obtained from a finite number of data by interpolation with radial basis functions. The values reported on the PES correspond to the corrugation, ∆γ, obtained by subtracting to the interlayer energy its minimum value, γ. The maximum value of the PES corrugation, ∆γ max , is obtained for the most unfavorable stacking configuration.
The ideal shear strength, which corresponds to the maximum resistance to sliding for an interface, is calculated as the derivative of the PES profile along high symmetry directions and along the minimum energy path (MEP) [51]. The MEP is automatically identified by means of the improved simplified string method, which is implemented in our Workflow [52]. When an interface is formed by mating two isolated surfaces, the charge spilling out from each surface is typically removed from the near-surface region and accumulated at the middle of the interface. These charge displacement ∆ρ is calculated by subtracting from the total charge of the bilayer, the charge of the single layers as ∆ρ(r) = ρ bilayer (r) − ρ up (r) − ρ down (r). The total amount of charge that is redistributed ρ redist = 1 2zo +zo −zo |∆ρ|dz , where z o is half of the interlayer distance, has been shown to be directly related to the interfacial adhesion [53] and is analyzed according to the procedure described in [48].

Structural and electronic properties
Black phosphorus is characterized by non-planar layers, where each P atom is covalently bonded with other three P atoms ( figure 1(a)). Having five valence electrons, each P atom presents a 3s lone pair which may determine a small covalent character also in the interlayer bonds.
The nature of the layer-layer interaction is analyzed in figure 2(a), where the charge displacement is calculated for the phosphorene bilayer. We observe that the charge spilling out from each surface is removed from the near-surface region and accumulates at the middle of the interface for all the employed XC functionals. The only exception to this behavior is with vdW-DF2 functional, for which a weakly pronounced  negative peak is observed. Even though this behavior is in agreement with previous diffusion Monte Carlo calculations [64], this feature may be dictated by the different description of the nonlocal kernel term compared to rVV10 one, which is the other nonlocal functional considered in this study, and leads to bulk lattice parameters significantly different from the experimental ones.
The charge accumulation here observed is in agreement with previous calculations performed both on black and blue phosphorous and employing the same XC functionals [64,65], and suggests that the interlayer bonds possess a certain degree of covalency. Although diffusion Monte Carlo simulations predicted a huge charge depletion, an interlayer covalent character can not be excluded a priori in layered materials. Indeed, recent synchrotron x-ray diffraction experiments showed that in TiS 2 the interlayer bonds have a larger covalent character than the theoretical predictions [66]. TMDs materials present non-planar structures, as phosphorene, which may play a role in determining the nature of interlayer bonds.
The bilayer binding energy, calculated within different schemes, is represented as a function of the interlayer distance in figure 2(b). As expected, both LDA and PBE are not able to describe the long range energy contribution: the former shows a pronounced increase of energy with distance, until unphysical positive values are reached above 5 Å , whereas the latter predicts atomically weak interlayer binding. By adding the Grimme-D2 correction to PBE, the description of long range interactions improves, and the potential well becomes eight times deeper. The considered nonlocal functionals, i.e. rVV10 and vdW-DF2, predict stronger binding especially at long range. The energy minimum of rVV10 is the lowest one, whereas vdW-DF2 has a similar binding energy to PBE + D2, but shifted to higher distances. The relative positions between peaks and the depth of minima are in agreement with those reported by Shulenburger et al [64].
The interlayer binding changes with the number of layers as can be observed in figure 3(a). In particular, for every functional adopted in this study, the layers have a higher binding energy in absolute value, when their number is increased, progressively approaching the correspondent bulk energy. At the same time, as shown in figure 3(b), the vdW schemes show a decrease in the interlayer distance when the thickness of the slab is increased, whereas the opposite behavior is observed for LDA. In the former approaches, the values of the interlayer distances are the closest to their bulk counterpart with a coherent energy-distance   trend, which leads to closer layers when the adhesion is higher. In multilayer structures the two outermost pairs of layers are found to be closer to the neighboring layer than the bulky pairs due to the presence of a boundary. This effect is enhanced by the use of LDA, which tends to overbind structures, as confirmed by the higher charge accumulation that is predicted for bilayer phosphorene ( figure 2(a)). It is clear from this analysis that the use of vdW corrections is crucial to properly describe the structure and the energetic of phosphorene films. The geometrical parameters calculated for black phosphorus and phosphorene are reported in table 1 for all the analyzed methods. Despite its simplicity, PBE + D2 is found to provide the best agreement with experiments and we thus adopt it for our subsequent analysis. The same computational choice has been made by our group to investigate the tribological properties of graphene and TMDs [32][33][34]. This allows for a close comparison between these materials.

Tribological properties
To evaluate the potentiality of phosphorene as solid lubricant, we calculate its interlayer slipperiness and compare it with multilayer graphene and MoS 2 , which are commonly used to reduce friction of solid interfaces. We consider two different layer orientations for the analysis of phosphorene tribological properties: a 'parallel' orientation ( figure 4(a)), where the zigzag and armchair directions of the two layers are aligned, and a 'perpendicular' orientation ( figure 4(b)) where the upper layer is rotated by 90° with respect to the lower one. This latter structure is modeled with a supercell having 4 × 3 in-plane size, where an initial geometrical mismatch of 0.5% is present between the two planar lattice vectors of the supercell. The mismatch has been reduced to 0.04% after we performed a variable-cell relaxation, where both the atomic positions and in-plane lattice cell vectors are optimized.
The calculated PES for both the parallel and perpend icular layer orientations are reported in figure 4. The relative position of the interfacial P planes is also represented from a top view for the stacking configurations corresponding to the minima and maxima of the PES.
The most energetically favorable stacking for the parallel orientation, AB, corresponds to a configuration where the atoms of the two interfacial P planes are in a hollow position and in the typical stacking of black phosphorus layers. This bilayer stacking is characterized by a binding energy of −0.35 J m −2 and an interlayer distance of 3.1 Å . The less favorable staking, AB , corresponds to an on-top disposition of the interfacial atoms, where the Pauli repulsion caused by the overlap of the electronic clouds is maximum, leading to an increase of binding energy and distance to −0.17 J m −2 and 3.7 Å , respectively. The transition from AB to AB stacking can occur during interlayer sliding along the armchair direction, while it never occurs along the zigzag direction. This gives rise to a friction anisotropy, as can be seen from the energy profiles reported at the bottom of the PES in figure 4(a). From the derivative of the energy profiles it is possible to obtain the ideal shear strength, that is the resistance of the system to fail to shear. The calculated values of 1.5 GPa and 1.2 GPa along the armchair and zigzag directions indeed confirm a friction anisotropy.
This anisotropy has been also measured experimentally during the sliding of an atomic force microscopy tip on phosphorene films 30 nm thick. The friction measured along the armchair direction always resulting higher than that measured along the zigzag direction [9]. This behavior has been attributed to the different values of Young's modulus along the armchair (27.2 GPa) and zigzag (58.6 GPa) directions [5], determining lower lateral deformation and friction along zigzag direction than the armchair direction. Here we add a new piece of information concerning the anisotropy of the interlayer sliding. This effect may play a relevant role especially in the situation where the tip drags a phosphorene flake and the sliding interface is shifted from the tip-phosphorene tribopair to the phosphorene-phosphorene one.
The sliding path with the highest statistical weight is the MEP that connects the PES minima passing through the saddle points. The MEP is marked in white on the 2D plots of the PESes in figure 4, and the corresponding energy profile is described by a red solid line below, with an associated shear strength of 0.85 GPa.
The PES obtained for the perpendicular orientation appears flat when compared to the PES for the parallel one on the same energy scale, indicating that all the relative lateral configurations in the former case are almost energetically indistinguishable ( figure 4(b)). In particular, the minimum and maximum energy stackings, referred to as AZ and AZ in figure 4(b), differ only by 0.01 J m −2 being ten times smaller compared to the parallel orientation. The friction anisotropy almost disappears and the shear strength along the MEP is reduced by one order of magnitude, becoming equals to 0.08 GPa. This strongly points at the existence of a structural superlubricity that may result in ultralow friction coefficients, as evidenced in a classical simulation of a sliding tip on phosphorene [69]. As can be observed in figure 5, the AZ (AB ) configuration is less stable than the AB configuration by 20% (51%). This smaller decrease in the adhesion energy for the AZ stacking is nevertheless accompanied by a reduction of one order of magnitude in the corrugation and shear strength.
In table 2 the interlayer tribological properties of black phosphorus are compared with those of graphite and MoS 2 [37,67,68]. The adhesion energy and shear strength of the parallel orientation are of the same order of magnitude compared to other 2D materials, being only slightly larger than those of MoS 2 . Instead, the perpendicular orientation presents a resistance to sliding of about one half compared to graphene. These results confirm that phosphorene may be successfully employed as solid lubricant, suggesting possible applications for this material in nanotribological devices operating in vacuum or dry conditions, as well as in aerospace environments, due to its instability when exposed to air [70]. Phosphorene has been successfully employed as lubricant additive in water-based solution, giving rise to a superlubrication regime [12].
Since the parallel and perpendicular orientations present a similar stability, we expect the presence of both orientations between dispersed phosphorene flakes and this may play a significant role in achieving a superlubricant regime.
A similar electronic charge distribution between the PES maxima and minima results in a smaller PES corrugation, as explained in [53]. It is evident from the planar average of the total electronic charge reported in figure 6 that the charge rearrangement occurring   when moving from the minima to the maxima is more consistent for the parallel layers than for the perpendicular ones as the amount of charge at the middle of the interface is much higher in AB than AB , while it is almost the same in AZ and AZ . For this configuration, we observe also a superimposable charge displacement, reported in figure 1 of the supplementary information.
A possible explanation of the superlubricant behavior of the AZ stacking can be found looking at the distribution of the electronic charge at the interface, as can be seen in figures 7(a)-(d), where the total charge present on a plane located at the middle of the interface is represented. The AB stacking shows a non uniform distribution of electronic charge, characterized by horizontal strips connecting the mated atoms. On the contrary, in AZ stacking the charge is present only nearby the atoms of the two layers in contact, and is almost zero in the interstitial regions. This strong localization determines the formation of sparse charge bridges between the two layers. Due to the incommensurate disposition of the phosphorene atoms within the layers, this grid-like pattern is maintained for all the possible relative lateral shifts between the two layers, making all the lateral configurations energetically equivalent.

Conclusion
In this work we investigated, for the first time to our knowledge, the static tribological properties of phosphorene by means of ab initio calculations. We compared different exchange-correlation schemes in describing the structural and electronic properties of films with a different number of layers, and found that PBE + D2 provides a description of black phosphorus in closest agreement with experiments.
The analysis of the interlayer shear strength showed that black phosphorus has similar tribological Figure 7. The electronic charge density on the xy plane at the middle of the interlayer region for the AB (a), AZ (b), AB (c), and AZ (d) stackings. The charge for the AB and AB configurations is represented on a replicated 4 × 3 cell, in order to compare it with the supercell adopted for the AZ and AZ ones. Phosphorene unit cell is also reported on the plots.
properties to other layered materials commonly used as solid lubricants. Moreover, a huge reduction in the PES corrugation and shear strength has been detected for a perpendicular orientation of the layers, where the armchair direction is aligned to the zigzag one. The stability of this superlubric configuration of the layers is similar to that of the parallel orientation, thus flakes containing both phases are expected in phosphorene-based lubricant solutions. The microscopic origin of the observed structural superlubricity has been explained by analyzing the electronic charge at the interface. The perpendicular layers are bonded by several localized bonds where the charge is localized. In this situation all the relative lateral positions become equivalent from the energetic point of view and sliding can occur easily on the flat PES. Our study represents an important step towards the full understanding of the atomistic origin of phosphorene lubricant properties, which can be of paramount importance not only for the tribological community, but also for applications in several micro and nanodevices.