Theoretical description of the physical/chemical contributions determining the stability of transition-metal doped phosphorene nanosheets

. Phosphorene is a two-dimensional semiconductor material with promising uses in nanotechnology, and its properties have been proposed to be tailored by transition metal doping; however, there is an incomplete understanding of the physical/chemical effects driving the stability of different dopants and doping schemes (adatom vs. substitutional doping). In this report, we explore the stability of doped-phosphorene nanosheets with first-row transition metals in the framework of density functional theory and by using a bonding characterization and energy decomposition analyses. The adsorption energies of Sc-Zn atoms onto phosphorene shows that the substitutional doping reaches higher stability in comparison to the adatom doping, where coordinate covalent bonding is the primary binding mechanism. The substitutional doping reaches the highest stability as a result of the strong interplay of charge-transfer, polarization, and Coulombic electrostatic effects. In contrast, adatom doping is mainly stabilized by charge transfer effects. It is also found that the magnitude of the polarization stabilizing effect is related to the polar character of the single closed-shell interactions between dopants and phosphorous atoms, while the fulfill of 3 d orbitals controls the magnitude of stabilizing charge-transfer energies. Moreover, the magnitude of the physical stabilizing contributions is associated with the electronic/structural properties of the doped phosphorene nanosheets, e.g., the topology of electron density, coordination environment, and ligand-field exerted on the dopant atom. These results allow us to expand the understanding of the physics and chemistry of low-dimensional phosphorene-based materials.


Introduction
Phosphorene is a layered low-dimensional nanomaterial composed of phosphorous atoms via strong P-P covalent bonds and weak interlayer Van der Waals forces that resembles a puckered hexagonal structure that is different of the flattened structure of graphene [1]. This structure has been successfully isolated employing liquid-phase exfoliation and mechanical cleavage of black phosphorus [2][3][4], and it shows interesting properties such as experimental stability, sizable direct bandgap (1.5 eV for monolayer phosphorene), and high electron mobility at room temperature (up to 1000 cm 2 V -1 s -1 ) [5]. The study of phosphorene-based materials is an increasing field of research, focusing on their applications in nanoelectronics devices, energy storage and solar cells, fieldeffect transistors, photocatalysis, optoelectronics, and photovoltaic devices [5][6][7][8][9][10]. Doped phosphorene have also been proposed as potential next-generation materials for spintronics, photoelectronic and optical devices, catalysis, storage, pollutant, removal, and gas sensing .
The introduction of impurities or dopants allows the tailoring of the phosphorene properties, where dopants can be experimentally introduced as in related low-dimensional nanomaterials, this is by chemical modification, intercalation, low-energy ion implantation, and defect-assisted doping by electron beam into irradiation [35,36]. In this sense, transition metals (TM) display distinctive properties across the periodic table (e.g., atomic radius, electronegativity), which are strongly influenced by the increasing number d electrons. These distinctive properties are used as a strategy to modify and modulate semiconductor, magnetic, and chemical reactivity properties [11,14,16,19,22,23,37].
The bonding between TM and phosphorene is reported to take place by covalent bonding in two different doping schemes: adatom and substitutional doping (see Fig. 1). In adatom doped phosphorene, the dopant atom is adsorbed at the surface of phosphorene in atom top, in-plane bridge, and hollow conformations, where the hollow position is the most stable [14][15][16]18]. In substitutionally doped phosphorene, the dopant atom replaces one phosphorous atom, forming stable chemical bonds with unsaturated P atoms [19,22,23,28,30,38]. In terms of stability, the adatom doping of phosphorene with TM has proven to achieve higher stable interactions (high magnitude of binding energies) compared to the doping of graphene with 3d TM atoms (which reach binding energies in the range of 0.3-3.4 eV) [39][40][41]. TM adatoms also show higher stability on phosphorene compared with non-transition metal atoms such as C, N, O, Li, Na, and Mg, which are bonded to phosphorene with a binding energy of 0.6-5.3 eV [12,14,15,18]. On the other hand, the substitutional doping of phosphorene with stable binding energies, and comparable to those on substitutionally doped graphene with metal and semimetallic atoms, where binding energies of 0.6-9.6 eV are reported [11,[42][43][44][45][46]. In this way, a wide series of atoms (including metals, alkaline metals, transition metals, metalloids, non-metals, and halogens) can be substitutionally introduced into phosphorene with binding energies of 2.1-9.2 eV [18,22,38].
From the structural point of view, TM and P atoms are strongly bonded, straining the phosphorene surface due to that the TM−P bond lengths are slightly shorter than the sum of their covalent radii. Consequently, TM dopants display different values of elevation, mostly accommodated inward into phosphorene and forming up to five bonds [23]. In the case of the TM adatom doping, the TM atoms are outwards from the phosphorene surface, where the dopant elevation decreases monotonically with respect to the dopant covalent radius [19]. Regarding the stability, the substitutionally doped phosphorene has higher binding energies (up to 8.3 eV) than the adatom doped phosphorene(up to 4.6 eV) [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], but both doping schemes are expected to be stabilized by chemical bonds (excepting Zn) as a result of the hybridization between d orbitals of TM dopants and 3p orbitals of phosphorene [14,16,19,22]. The wide range of binding stability reached by TM dopants in phosphorene has been explained due to the filling of the d orbitals through the TM series, causing several changes in the properties of phosphorene such as bandgap decreasing due to midgap states, electron charge-transfer from TM to the phosphorene nanosheet due to electronegativity differences, and induced magnetic moments by spin-polarized dopants [11, 12, 14-16, 19, 20, 22, 23]. Despite the latter, there is no major information about the specific role of physical/chemical phenomena driven the doping stability among different dopants (such as orbital interactions, electrostatic, polarization, dispersion, charge transfer, and/or Pauli repulsion).
It is clear that the number of 3d electrons cannot be the only key factor determining the dopant stability in phosphorene nanomaterials, but the study of physical/chemical phenomena determining the dopant stability remains some unknown.
Considering the interest in the understanding of the properties of phosphorene-based materials, this work is aimed to clarify the origin of the quantitative and relative stability of TM dopants on phosphorene by density functional theory. We focus the attention on the physical/chemical phenomena dominating the doping stability in adatom and substitutional schemes. From the physical point of view, energy decomposition analysis based on absolutely localized molecular orbitals (ALMO-EDA) will allow us to characterize the binding strength in terms of physical terms easy to interpret (electrostatic, polarization, dispersion, charge transfer, and Pauli repulsion). From the chemical point of view, the atoms-in-molecules (AIM) and natural bond orbital (NBO) analyses will be used to obtain the chemical fingerprint of the bonding (electron density at the bond critical points and strength of donor-acceptor orbital interactions).

Computational methodology
A finite phosphorene nanosheet (P125H30) was adopted (surface size of ⁓440 Å 2 and volume of ⁓1012 Å 3 ), where the dangling bonds at the edges were saturated with hydrogen atoms to avoid vacant orbitals. The selection of this model is based on the local nature of the doping and well-converged dopant binding energies; also, a molecular model allows to use an extended allelectron basis set, and it is required for the implemented binding and energetic analysis in this work (see below). Two doping schemes were modeled: substitutionally (S) doped phosphorene (TM-P125H30, abbreviated as TMSPhos), and adatom-doped phosphorene (TM-P126H30, abbreviated as TMAPhos) (Fig. 1). Third-row transition metals (Sc-Zn) were adopted as dopants.
The doped phosphorene nanosheets (TMPhos) have a dopant concentration of up to 1.6% to agree with related theoretical studies [16,23]; note that increase in dopant concentration is reported to slightly affect the binding energies in doped-phosphorenes [18]. DFT calculations were performed in the ORCA4. code [47,48] with the PBE [49] functional. The all-electron triple-split valence and polarized def2-TZVP basis sets were used for core and valence electrons of all the atoms, this is for the best quality of our results [50]; then, pseudopotentials were not implemented. The DFT-D3 method with the Becke-Johnson damping function was used for dispersion corrections on the SCF energies and gradients [51,52]. Basis set superposition errors (BSSE) were corrected with the Boys and Bernardi method [53]. respectively [54]. Atoms-in-molecules (AIM) and related wavefunction analyses were performed in the Multiwfn 3.6 program [55]. In addition, the orbital picture of the bonding was analyzed employing the natural bond orbital (NBO) module in Gaussian16 [56].
Dopant atoms form a characteristic coordination arrangement in the TMPhos systems. The geometrical distortion (gd) of these coordination arrangements with respect to ideal solids was characterized using the Shape 2.1 program [57], where the gd function is obtained as (Eq. 1): Then, the terms gd can be approximately given by an ideal or distorted structure/polyhedron composed of N vertexes. Q0, Qk, and Pk correspond to the center of mass, distorted, and perfect vector coordinates, respectively. The gd function is obtained by the minimization between the difference of distorted coordination and perfect coordination vectors relative to the center of mass vectors of a molecular fragment. In this sense, the distortion of the coordination sphere is analyzed from a general point of view without dependence on the bond length but rather from the deviation from ideal geometries. On the basis of such definitions, the boundaries for any symmetry measure are 100 ≥ gd ≥ 0, where lower limits stand for structures that match the target symmetry, and higher values stand for distorted structures [58].
The binding energy (Eb) is the energy required to pull-out the dopant atom from phosphorene, characterizing the stability of the dopant-phosphorene systems (Eq. 2): where ETM, EPhos, and ETM-Phos are the total energies of the isolated TM atom, isolated Phos nanosheet, and the TMPhos nanosheet, respectively. Then, the more positive the Eb value, the more stable the TMPhos system is. Binding energies were further decomposed into physical contributions by the second-generation energy decomposition analysis based on absolutely localized molecular orbitals (ALMO-EDA) in the Q-Chem 5.2 program [59][60][61][62][63]. In this scheme, the E b value of an A-B complex is decomposed into the sum of six physical terms as (Eq. 3): where EELEC represent the permanent Coulombic electrostatic interactions between fragmental charge distributions. EPOL is polarization energy term, resulting from the intramolecular density relaxation due to the perturbating presence of the other fragment in the complex. ECT is the charge transfer energy term, accounting for the stabilization due to the inter-and intra-molecular charge flow in the fragments after the complex formation [61]; note that ECT is not the number of electrons transferred between fragments. EDISP stands for stabilizing dispersion forces. EPREP is the preparation energy, tracking the energy penalty due to the geometric/electronic preparation of the fragments to reach the complex geometry. Finally, EPAULI stands for the steric repulsion of electrons with the same spin when two fragments are compressed (Pauli principle).  (Table 1).    3a-c show that in the substitutionally doped phosphorene nanosheets (TMSPhos), the metal atom form up to four TM-P bonds (d1, d2, d3, and d4 in Fig. 3b) with bond distances ranging from 2.2 to 2.9 Å (Fig. 3c). Moreover, a fifth metal-P bond (dh) is perpendicular to the plane of the d1-4 bonds in a downwards direction not longer than 2.6 Å; the bonding can be associated to a metal-centered square pyramid geometry based on a TM-P5 coordination unit as noted from the polyhedral representation (Fig 3a). The exceptions are the CuSPhos and ZnSPhos systems, where dopants form up to two bonds with bond length values of 2.2-2.3 Å. It is also noticed that d3 and d4 TM-P bond distances are longer than 2.9 Å (Fig. 3c), resembling the bond distance that the P atom has on pristine phosphorene. In the case of the adatom doped phosphorene (TMAPhos), the adatom prefers to stay at the hollow site of one hexagonal ring and bonded with three P atoms [12,16,18,19].

Binding energies
The TMAPhos systems have a trigonal pyramidal shape in the dopant region ( Fig. 3d; TM-P3), displaying three TM-P bond lengths (d1, d2, and d3) in the range of 2.2-2.9 Å (Fig. 3e). d1 is equivalent to d 2 , and d 3 differ up to 0.13 Å compared to d 1 -d 2 (Fig 3f). In addition, the dopants are protruded over the phosphorene plane from ~1.05 to ~2.08 Å (dh) according to the increase in the atomic radius of the isolated dopant atom [64]. Finally, the highest TM-P bond length (~2.8 Å) and elevation (~2.1 Å) in the ZnPhos system agree with its lowest binding energy.

Energy decomposition analysis
The physical contributions to the stability of the TMPhos systems were grouped as stabilizing (E ELEC, E DISP, E POL , E CT ), and destabilizing terms (E PREP , E PAULI ) in Table 2; while Figures 4a-d display the relative single percentage contributions. In general, it is noted that the stability of TMSPhos systems is dominated by the interplay between electrostatics, polarization, and charge transfer effects. In contrast, the stability of the TMAPhos systems is dominated by charge transfer effects. In all the cases, the destabilizing effects are mainly coming from Pauli repulsion (EPAULI). A detailed analysis of each contribution is given below.

Charge transfer (ECT)
The charge transfer is the dominant stabilizing physical effect, this is inter-fragment orbital interactions play a key role in the stability of doped phosphorene materials. For instance, the ECT values in the substitutionally-doped TMSPhos systems are in the range of -8.89 to -4.65 eV ( Table   2). Since ECT is related to the magnitude of the donor-acceptor orbital interactions responsible for a net charge flow between molecular fragments, the NBO analysis was used to locate the interactions between occupied and unoccupied orbitals. In this way, the magnitude of the ECT term is related to a mixture of electron delocalization from the occupied (bonding) orbitals of phosphorene to the unoccupied (antibonding) acceptor orbitals of metal atoms [such as dLP*, pLP* and d Rydberg (RY*)]. Occupied donor orbitals of phosphorene are bonding P-P orbitals (the so-called BD orbitals), sp lone pairs (the so-called LP orbitals), and s-core orbitals (the socalled CR or non-valence orbitals).
The role of the donor-acceptor (bonding-antibonding) interactions is exemplified by two limiting cases: ScPhos and ZnPhos (Fig. 5). In the case of the substitutionally-doped ZnSPhos system, the interaction shows phosphorene→metal electron delocalization of the BD→pLP* and LP→pLP* type (Fig. 5a, right); this is the P-P bond, and 3p lone pair orbitals of phosphorene and the low-occupied 4p orbitals of Zn (in n=4) (Fig. 5b, right). The energy mismatching between donor and acceptor orbitals decreases the stability gained by charge-transfer effects. Conversely, the donor-acceptor interactions in ScAPhos system are due to LP→dLP* orbital interactions (Fig.   5b, left), this is between the occupied 3p lone pair orbitals of P atoms of phosphorene (in n=3) and the low-occupied 3d orbitals of Sc (in n=3) (Fig. 5b, right). Then, the energy matching between donor and acceptor orbitals increases the gained stability by charge-transfer effects.  Additionally, the ECT term was compared against the electron density difference [Δρ(r)], which shows the sites with accumulation (red surfaces) or depletion (blue surfaces) of electron density after the doping (electron density redistribution). Fig. 6 shows the Δρ(r) function for all the systems (at the same isovalue) and the magnitude of the ECT term. After doping, the accumulation of electron density is located at the TM-P bonds, while the electron density depletion is located at the TM atom as a result of the hybridization of the metal orbitals during the bonding.
The higher magnitude of the ECT term is qualitatively related to the larger Δρ(r) surface. For instance, the ScSPhos and MnSPhos systems (relative higher |ECT| values of 8.9 and ⁓8.4 eV, respectively) show a relatively high electron density redistribution. Conversely, the CuSPhos and respectively, which agree to the low magnitude of the ECT values and smaller Δρ(r) surface. The same interpretation is obtained from the analysis of the adatom-doped phosphorene nanosheets.

Permanent electrostatic interactions (EELEC)
The EELEC term is related to the Coulombic electrostatic interactions between dopants and phosphorene due to its permanent moments of partial charges, showing an R -6 dependence (where R is the distance between fragments) [65]. The electrostatic interactions are attractive in all the cases, this is only negative EELEC values are obtained (Table 2). In the case of substitutionallydoped TMSPhos nanosheets, the EELEC energy is strongly determined by the displacement of the dopant atom around the hollow site of the phosphorene lattice (Fig. 3); the latter changes all the bond distances and perturbs its coordination environment. This dependence is due to that the electrostatic interactions are inversely proportional to the atom-atom bond distances r12 [EELEC1/r1-2]. Otherwise, the crystal field theory states that the five d-orbitals of free TM atoms have the same energy. In the presence of a ligand (in this case, the ligand is phosphorene), dorbitals are perturbed by the coordination environment, inducing the energy splitting of d-orbitals due to the electrostatic environment of the ligand (attractions and repulsions) [66]. Therefore, the distortion on the coordination geometry will affect the crystalline field on the TM atom, determining the EELEC values due to its dependence on the interatomic distances r1-2. Considering the above, Fig. 7a shows that the geometrical distortion (gd) of the square pyramid geometry (solid orange line) affects the permanent electrostatic moments in the TMPhos nanosheets, and subsequently, the EELEC energy (solid blue line). In fact, the magnitude of the stability gained by electrostatic effects (EELEC) is decreased as the distortion increases (Fig. 7b).
For instance, the relatively high distorted geometry of the NiSPhos, CuSPhos, and ZnSPhos system (resembling a quasi-trigonal pyramid geometry) is associated with a decrease in the magnitude of the EELEC term, this is in the stability gained by electrostatic effects. Otherwise, the adatom-doped TMAPhos systems hold a trigonal pyramid geometry, where the geometric distortion can be simply attributed to the change in the dopant elevation (dh) with respect to the upper phosphorene sublattice. Fig. 7c shows the dh parameter (solid green line) and the magnitude of the EELEC term  Fig. 7d, which is consistent with the 1/r 1-2 dependence of electrostatic interactions.

Polarization (EPOL)
The induced electrostatic or polarization term EPOL accounts for the on-fragment NiAPhos systems show <ρBCP> values of ~0.10 e/Bohr 3 , which are in the range of the highly polarized interactions (coordinate covalent bonds). In this regard, the magnitude of the EPOL energy is increased as the <ρBCP> value is increased (Fig. 8b). Conversely, the lowest stable system ZnAPhos shows a low <ρBCP> value of ~0.03 e/Bohr 3 , which is consistent with a non-covalent interaction of electrostatic character; as a consequence of the non-covalent state, the magnitude of the EPOL value is decreased (EPOL=-0.62 eV). It is important to point out that the EPOL terms also can be associated with the dh parameters since the ability of one fragment to polarize another fragment is depending on the distance between them. In addition, the relation of the EPOL term to chemical and structural properties of the substitutionally-doped systems is less intuitive due to the dopants and phosphorene are in close contact; in this scenario, the tagging of electrons to fragments is ambiguous and similar difficulties arise in the definition of intra-fragments [65].

Dispersion interactions (EDISP)
Since the dispersion energy (EDISP) shows an asymptotic distance dependency (R -6 ) [65], the contribution of dispersion forces is significantly reduced (0.30-0.59 eV) by the shorter TM-P distances as seen in Fig. 3. Then, dispersion energy does not play an important role in the stability of the doped systems, with a single percentage contribution of 2-8%.

Pauli repulsion and preparation energy (EPAULI and EPREP)
The destabilizing energy in TMSPhos and TMAPhos systems mainly comes from the repulsive interaction due to volume exclusion effects EPAULI (up to 99%, Figures 4b and 4d). The destabilizing EPAULI energy is higher in the TMSPhos systems compared to those in the TMAPhos nanosheets (Table 2), where the minimum EPAULI values are ~9.8 and ~3.7 eV, respectively.
However, this difference in the magnitude of the repulsive interactions cannot be associated with the differences in stability between substitutional and adatom doping, because the relative stability among doping schemes is mainly determined by stabilizing effects as noted above.
Finally, the energy penalty due to the geometric/electronic rearrangement of the fragments to reach the bonded state (EPREP) must increase when a higher order of hybridization takes place in the chemical bonding [67,68]. This statement applies in the case of the MnSPhos system, where the contribution of EPREP is up to 41% as a result of the change in the number of unpaired electrons of the isolated metal for its spin state in the MnSPhos system [16]. In other words, the resulting changes in the number of unpaired electrons can be associated with the difference of EPREP term values between doped systems (TMSPhos having one more unpaired electron than the TMAPhos complexes).

Conclusions
In summary, our results show that first-row transition metal dopants (Sc-Zn) bind with up to five phosphorous atoms in phosphorene, where the substitutional doping reaches higher stability with respect to the adatom doping, where the bonding is mainly of coordinate covalent character.
The highest stability of the substitutional doping emerges from the interplay between chargetransfer, polarization, and electrostatic stabilizing physical effects. Otherwise, the adatom doping