Mechanistic insights on the phosphorene degradation

TL;DR

This study provides detailed mechanistic insights into the structural and chemical degradation processes of phosphorene, using first-principles calculations to develop a comprehensive degradation phase diagram and identify key atomic-level mechanisms.

Contribution

It introduces a first-principles based phase diagram for phosphorene degradation, elucidating microscopic kinetic pathways and temperature effects on structural and chemical stability.

Findings

Degradation involves vacancy merger and annihilation at 400 K.

Structural degradation occurs via atom emission and sublimation above 650 K.

Chemical degradation is driven by oxygen dissociation at defect sites.

Abstract

The structural and chemical degradations of phosphorene severely limit its practical applications despite the enormous promise. In this regard, we investigate a plethora of microscopic kinetic mechanisms and develop a degradation phase diagram within the first-principles calculations. At 400 K, the degradation and the competing self-annealing proceeds through the merger and annihilation of vacancies, respectively, which are triggered via itinerant vacancy and adatom. A further increase in temperature beyond 650 K, the structural degradation results through the emission of the undercoordinated atoms from the defect and the concurrent pair-wise sublimation. The role of inter-layer vacancy diffusion is discarded in the context of structural degradation. The chemical degradation is routed through the dissociation of oxygen molecule that is either activated at the room-temperature on the pristine surface or spontaneous at the single-vacancy site. The present results are in agreement with the few available experimental conjectures and will motivate further efforts.

Tables

Table 1. Table 1: Formation energies E f superscript 𝐸 𝑓 E^{f} of the various defect complexes in SLP. The E f superscript 𝐸 𝑓 E^{f} is calculated as, E f = E tot − N P × E P superscript 𝐸 𝑓 subscript 𝐸 tot subscript 𝑁 P subscript 𝐸 P E^{f}=E_{\rm tot}-N_{\rm P}\times E_{\rm P} , where E tot subscript 𝐸 tot E_{\rm tot} is the total binding energy of the defect-containing supercell, and E P subscript 𝐸 P E_{\rm P} is the binding energy per P-atom in the pristine SLP. N P subscript 𝑁 P N_{\rm P} is the number of P-atoms in the defected lattice. The activation energies E a subscript 𝐸 𝑎 E_{a} , for the cohort of microscopic processes leading to the self-healing, degradation, and sublimation are calculated. Note that both vacancy and P-adatom P A are highly itinerant due to lower E a subscript 𝐸 𝑎 E_{a} and the corresponding diffusion are highly anisotropic along the zigzag and armchair directions.

Defect complex		Formation energy		Microscopic process		Activation barrier
		$E^f$ (eV)				$E_a$ (eV)
SV(5$|$9)		1.68		SV(5$|$9) $\to$ SV(5$|$9)		0.12 (zigzag), 0.38 (armchair)
SV(55$|$66)		2.03		$-$		$-$
DV(5$|$8$|$5)		1.56		DV(5$|$8$|$5) $\to$ SV(5$|$9) + SV(9$|$5)		1.64
SV(5$|$9) + SV(9$|$5)		2.52		SV(5$|$9) + SV(9$|$5) $\to$ DV(5$|$8$|$5)		0.68
Single PA		1.43		PA $\to$ PA		0.19 (zigzag), 1.05 (armchair)
				PA desorption		2.14
PA-pair		2.01		PA-pair desorption		barrierless
SV(5$|$9) + PI		1.56		SV(5$|$9) + PI $\to$ pristine SLP		0.69
				pristine SLP $\to$ SV(5$|$9) + PI		2.25
SV(55$|$66) + PI		1.93		SV(55$|$66) + PI $\to$ SV(5$|$9) + PI		0.10
DV(5$|$8$|$5) + PA		2.68		DV(5$|$8$|$5) + PA $\to$ SV(5$|$9)		0.46
				SV(5$|$9) $\to$ DV(5$|$8$|$5) + PA		1.46

Full text

I Introduction

Few-layered black phosphorus (BP) has garnered significant attention due to a layer-tuneable band gap between 0.3–2.0 eV and anisotropic mechanical, electronic and optical responses. (Liu et al., 2014a; Li et al., 2014; Xia et al., 2014; Zhang et al., 2014; Luo et al., 2015; Fei et al., 2014; Rodin et al., 2014; Tran et al., 2014) The staggered hexagonal layers of BP are stacked together by van der Waals interaction, similar to graphite. The isolated layers of BP, phosphorene, offers significant advantages in realizing novel semiconducting and optoelectronic devices due to its high current on-off ratio, extraordinary carrier mobility, ambipolar transport, and anisotropic electronic dispersion. (Liu et al., 2014a; Li et al., 2014; Das et al., 2014; Tran et al., 2014; Wang et al., 2015a; Liang et al., 2014) Moreover, the high mechanical flexibility and the retention of extraordinary electronic and optical properties under moderate strain allow to design functional heterostructures with other two-dimensional materials. (Rodin et al., 2014; Peng et al., 2014; Yuan et al., 2015; Muduli et al., 2018) Further, many exotic quantum many-body states have been theoretically predicted and experimentally demonstrated in few-layer phosphorene. (Liu et al., 2015a; Kim et al., 2015; Babar and Kabir, 2018a)

However, the major hurdle is the inferior stability of phosphorene at ambient conditions compared to graphene and transition-metal dichalcogenides. The degradation in two-dimension results from the dynamics of lattice defects and chemical interaction with the environment. (Queisser and Haller, 1998; Diebold, 2003; Janotti and Van de Walle, 2007; Song et al., 2017) Both surface and subsurface defects have been experimentally observed in the few-layer phosphorene due to their low formation energies, and are responsible for intrinsic p-type conductivity. (Liu et al., 2014a; Kiraly et al., 2017; Riffle et al., 2018; Wang et al., 2015b) Similar to graphene, the vacancy defects in phosphorene relaxes in multiple configurations, (Babar and Kabir, 2016; Cai et al., 2016) which act as anisotropic scatterers and have been further manipulated to form a many-body state with excited carriers. (Kiraly et al., 2017; Qiu et al., 2017) In addition to the disruption in the electronic properties, the vacancies in phosphorene induce mechanical degradation as well as produce structural changes under increasing strain, temperature, or vacancy concentration. (Liu et al., 2015b; Sha et al., 2016; Hao and Chen, 2016) Vacancy containing strained phosphorene reduces the fracture strength and structural integrity along the zigzag axis due to bond distortion and breakage. Further, regions with high defect density are predicted to initiate crack formation followed by the structural failure under the transverse force generated by an atomic force microscopy nanoindentation. (Sha et al., 2017)

Thermally activated diffusion of point-defects holds the key to structural changes. Vacancy transformation, diffusion, and aggregation leading to the formation of complex defects and grain boundaries have been extensively studied in graphene, (Skowron et al., 2015; Babar and Kabir, 2018b) where the weak bonding at the edges and vacancies initiate structural degradations. The few-layer phosphorene is observed to undergo anisotropic degradation, amorphization, and sublimation when heated above 650 K. (Liu et al., 2015b; Yoo et al., 2018; Fortin-Deschênes et al., 2016) However, the current description of structural degradation stands at the diffusion of isolated vacancy within a single layer, (Hu and Dong, 2015; Cai et al., 2016; Li et al., 2016) which severely limits our microscopic understanding. In this regard, here we investigate a cohort of mechanisms involving complex defect-defect interaction, vacancy and adatom migration, surficial desorption, which lead to larger defect formation. The competing self-healing mechanisms of defects are also investigated. Such a comprehensive understanding of degradation is necessary to identify the optimal experimental conditions to prepare and restore the defect-free phosphorene to improve its ambient stability.

Rapid oxidation of phosphorene in the ambient environment presents another significant challenge towards its stability. Although the exact mechanism continues to be investigated, it is understood that the strong chemical bonding between lone pair electrons of phosphorus and $\pi^{*}$ electrons of oxygen initiate the chemical degradation. (Castellanos-Gomez et al., 2014; Favron et al., 2015; Huang et al., 2016; Ziletti et al., 2015; Zhang et al., 2018; Abellán et al., 2017; Kuntz et al., 2017; Yang et al., 2015; Zehua et al., 2017; Qionghua et al., 2016) Further, degradation is accelerated as the hydrophilicity of phosphorene increases with oxidation. (Castellanos-Gomez et al., 2014; Ziletti et al., 2015; Qionghua et al., 2016) Moreover, the edges and steps undergo rapid degradation in few-layered samples, where the role of vacancy defects is understandably negligible. (Zhang et al., 2018; Kuntz et al., 2017) It is also observed that the oxidation rate increases with the reduction in the number of layers, (Favron et al., 2015; Abellán et al., 2017) which indicates a stronger interaction of oxygen at the surface.

However, the nature of oxygen-defect interaction in phosphorene is a subject of debate at present without any experimental understanding. While the oxidation is proposed to be an activated process, it is more favourable at a lattice vacancy than the pristine surface. (Kistanov et al., 2017) On the contrary, a weak interaction was suggested as the vacancy diffusion is found to be unaffected by oxygen. (Cai et al., 2016) In the present study, we thoroughly investigate the interaction of O2 on the pristine and defected phosphorene. Further, the microscopic dissociation mechanism is elucidated. We find the O2-dissociation to be a barrier-less process at the single vacancy defect, whereas oxidation on both pristine surface and at the divacancy defect remains an activated process. Once the dissociation is complete, the strong P$-$O binding essentially makes O-removal from the lattice impossible.

II Computational details

The spin-polarized density functional theory calculations are performed within the projector augmented wave formalism, (Blöchl, 1994) as implemented in the Vienna ab initio simulation package (VASP). (Kresse and Hafner, 1993; Kresse and Furthmüller, 1996) The wave functions are described in the plane wave basis with 500 eV cutoff for the kinetic energy, and the exchange-correlation energy is expressed with the Perdew-Burke-Ernzerhof (PBE) functional. (Perdew et al., 1996) Structural relaxations are performed until all the force components fall below 0.01 eV/Å  threshold. Calculations on the single-layer phosphorene (SLP) are performed using two supercells with different sizes, 6 (zigzag) $\times$ 4 (armchair) supercell containing 96 atoms and 9 (zigzag) $\times$ 5 (armchair) supercell containing 180 atoms. The reciprocal space integration for the two supercells was evaluated with 4$\times$4$\times$1 and 3$\times$4$\times$1 k-point grid respectively. A vacuum of 15 Å  was always maintained perpendicular to the surface to minimize the spurious interaction between the periodic images. Calculations for the few-layer phosphorene are performed using the 6 $\times$ 4 $\times$ $N$ supercells ($N$ $\leqslant$ 3). The nonlocal van der Waals functional (optB88-vdW) of Langreth and Lundqvist is used throughout the calculations to describe the oxygen-phosphorene and interlayer interactions. (Dion et al., 2004) The activation barriers are calculated within the climbing image nudged elastic band method formalism. (Henkelman et al., 2000) The true nature of the transition states is confirmed by obtaining the single imaginary vibrational mode.

III Results and Discussion

The two half-layers in the hexagonal staggered lattice of SLP are separated by 2.13 Å  (Figure 1). The in-plane atoms are arranged in a zigzag fashion with 2.23 Å  bonds (zigzag axis), whereas the atoms across the two half-layers have slightly larger bond lengths of 2.27 Å (armchair axis). Thus, every P-atom is bonded with two in-plane and one out-of-plane atoms. Given that each P-atom has five valence electrons, three electrons participate in bonding, whereas rest of the electrons produce a lone-pair charge cloud. The Coulomb repulsion among the lone-pairs leads to a staggered, anisotropic arrangement of P-atoms with comparatively soft bonding (Figure 1). Thus, Young’s modulus for phosphorene is substantially lower than graphene. (Liu et al., 2007; Wei and Peng, 2014)

Before we discuss the kinetics of point-defects, we investigate the formation of lattice vacancy in the light of recent experimental observations. The removal of a single P-atom introduces dangling electrons on the undercoordinated atoms, which concurrently undergo differential structural rearrangements to form two distinct defect configurations, SV(5$|$9) and SV(55$|$66) in Figure 1(a) and 1(b). The SV(5$|$9) configuration involves inward displacement of two atoms along the armchair direction to form a new P–P bond between the two half-layers, while the third P-atom remains undercoordinated [Figure 1(a)]. In comparison, the vacancy relaxation in graphene results in multiple Jahn-Teller distorted configurations, where the planar (5$|$9) vacancy is the ground state. (Babar and Kabir, 2018b) The second configuration SV(55$|$66) in phosphorene undergoes minor displacement of atoms to form pentagon and hexagon pairs [Figure 1(b)]. The P–P bond along the armchair direction (2.57 Å) is weaker than both pristine (2.27 Å) and SV(5$|$9) defect with 2.41 Å  bonds. Such four-fold coordination also appears in rippled graphene due to increased sp3 character. (Santos et al., 2012) The SV(5$|$9) vacancy is thermodynamically more stable than the SV(55$|$66) configuration [Figure 1(c) and Table LABEL:tablepvS1]. Further, owing to a much lower $E^{f}$, the vacancy formation is much easier in phosphorene than in graphene ($\sim$ 7.5 eV) (Babar and Kabir, 2018b) and in mono-elemental bulk semiconductors such as Si and Ge ($E^{f}$ = 3.0–4.5 eV). (Śpiewak and Kurzydłowski, 2013; Weber et al., 2013) This indicates phosphorene to be defect-prone, and thus the point-defect driven structural degradation becomes relevant. Further, the vacancy formation energy at the surface of a few-layer phosphorene decreases with the increasing layer thickness [Figure 1(c)]. While the relative stability of SV(5$|$9) and SV(55$|$66) differs by 350 meV in single-layer, it decreases to 48 meV in bulk BP. Similarly, the calculated $E^{f}$ for the subsurface SV(5$|$9) and SV(55$|$66) defects in 3L phosphorene is much lower, 1.54 and 1.60 eV, respectively, compared to the SLP (Table LABEL:tablepvS1).

III.1 Anisotropic vacancy migration

A complete understanding of migration mechanisms of a lattice vacancy in phosphorene becomes important due to their significant presence owing to their low formation energies, and as they drive mechanical degradation. (Riffle et al., 2018) Here, we investigate both intra- and inter-layer diffusion and consider jumps along both the armchair and zigzag directions to understand the effects of structural anisotropy. Relating the degradation mechanisms with the vacancy migration along with the adatom mediated self-healing will assist in constructing a better defect annealing strategy.

The SV(5$|$9) migration in the SLP is anisotropic along the armchair and zigzag directions with 0.38 and 0.12 eV, respectively [Table LABEL:tablepvS1 and Figure 2(a)-(b)]. Along the armchair direction the vacancy jump occurs via switching of bonds – the P-atom of the pentagon bond is displaced along the armchair direction and subsequently bonds with the undercoordinated atom [Figure 2(a)]. The undercoordinated atom then migrates to the bottom half-layer in the final configuration. The vacancy migration along the zigzag direction occurs analogously shown in Figure 2(b); however, the migration is facile due to a much lower activation barrier. These results are in good agreement with the previously reported values of 0.18$-$0.40 eV. (Li et al., 2016; Hu and Dong, 2015; Cai et al., 2016)

In addition to the anisotropic migration, the vacancy jump mechanism discussed so far is peculiar since the undercoordinated atom plays a minor role in the migration. Typically, the undercoordinated atom can detach easily and jump to the vacant lattice site. However, during SV(5$|$9) vacancy migration, phosphorene retains the staggered arrangement by switching of P–P bonds across the half-layers and lowers the activation barrier. The overall diffusion occurs through a combination of anisotropic vacancy jumps, and the diffusion along the armchair direction is the rate-limiting step with 0.38 eV barrier. In contrast, a relatively high barrier of 0.72 eV was calculated for the vacancy migration in graphene. (Babar and Kabir, 2018b) The aggregation of vacancy defects in an anisotropic environment will lead to the formation of line defects and grain boundaries, which in turn will act as a sink for diffusing vacancies. (Li et al., 2016)

Now we investigate the effect of layer thickness on the in-plane vacancy diffusion, and also discuss how the migration is affected while the vacancy is embedded in few-layer phosphorene (Figure 3). In 2L phosphorene, the activation barrier decreases to 0.28 and 0.07 eV, respectively, for diffusion along the armchair and zigzag directions (Figure 3). We do not anticipate further reduction in $E_{a}$ with increasing thickness. While the mechanism remains the same, the subsurface vacancy diffusion is easier as revealed by the lower activation barrier (Figure 3). For subsurface SV(5$|$9) migration in 3L phosphorene, the calculated $E_{a}$ are found to be 0.22 and 0.04 eV along the armchair and zigzag directions, respectively. Therefore, in 3L phosphorene, the $E_{a}$ for the subsurface diffusion converges to bulk BP, which we find to be 0.2 and 0.04 eV, respectively (Figure 3). These results are in contrast to the few-layer graphene and graphite, where the activation barrier is higher for the subsurface vacancy diffusion, (Babar and Kabir, 2018b; Asari et al., 1993; Paredes et al., 2009) and the microscopic origin has been discovered only recently. (Babar and Kabir, 2018b) In single-layer graphene, the strain generated in the lattice due to the vacancy diffusion is released by an out-of-plane buckling, whereas the corresponding strain in phosphorene is released by the lattice relaxation within the half-layers. As the out-of-plane buckling is restricted in multi-layered graphene, the $E_{a}$ for the surface and subsurface vacancy increases substantially. In contrast, the lower $E_{a}$ in few-layer phosphorene is due to the reduction in the corresponding formation energies as discussed earlier [Figure 1(c)]  (Wang et al., 2015b; Guo and Robertson, 2015) A similar Brønsted-Evans-Polanyi-type correlation between $E^{f}$ and $E_{a}$ has been reported in other bulk semiconductors and carbon nanotubes. (Carrasco et al., 2006; Cheng and Selloni, 2009; Kabir and Van Vliet, 2016)

While the $E^{f}$ for the subsurface SV(5$|$9) defect is lower than the same on the surface, it would be interesting to investigate the subsurface vacancy diffusion to the surface. Due to the AB-staking in few-layer phosphorene, the vacancies in adjacent layers do not align vertically and thus diffuses via inter-layer diagonal jumps [Figure 2(c)]. We considered the two possible pathways, where either of the nearest and next nearest neighbour atom from the adjacent undefective layer migrates to the SV(5$|$9) vacancy. The migrating P-atom, in both mechanisms, forms a Wigner-like defect, where the moving P-atom bonds with both active layers [Fig. 2(d)]. (Telling et al., 2003) While the inter-layer SV(5$|$9) migration via the nearest P-atom jump is calculated to be 2.17 eV, the second mechanism costs much higher activation energy of 2.80 eV. Such substantial inter-layer barrier restricts the vacancy migration across the layers to very high-temperature and does not play any significant role in the annealing process, which is usually performed at temperatures around 475 K. (Liu et al., 2015b) In comparison, the corresponding barrier is significantly higher in few-layer graphene, 5.5–7.0 eV, and becomes relevant only above 2000 K. (Trevethan et al., 2014; Liu et al., 2014b) However, the inter-layer vacancy migration may occur at a much lower temperature in the presence of complex inter-layer defect structures. The similar increase in inter-layer activity is reported in the few-layer phosphorene above 600 K, while such inter-layer defect complexes are present. (Yoo et al., 2018)

III.2 Vacancy-driven degradation, void formation and sublimation

The recent experiments indicate that heating of few-layer phosphorene creates oblate voids in the individual layers, which is followed by the inter-layer P–P bond formation. Further, the degradation proceeds via amorphization and eventual sublimation between 650–680 K. (Liu et al., 2015b; Yoo et al., 2018; Fortin-Deschênes et al., 2016) Thus, to achieve a better understanding of structural degradation, we investigate a cohort of mechanisms such as self-healing, [SV + SV] $\rightarrow$ DV, SV $\rightarrow$ [DV + PI] pair transformations, and sublimation (Table LABEL:tablepvS1).

The divacancy defect DV(5$|$8$|$5) is formed by the removal of two adjacent P-atoms from the same half-layer. The concurrent structural relaxation to pentagon-octagon-pentagon arrangement ensures three-fold coordination for all the atoms with a formation energy of 1.56 eV, which is lower than the SV(5$|$9) defect formation. Although the DV(5$|$8$|$5) has not been yet observed in the low-temperature experiments below 10 K, (Kiraly et al., 2017; Riffle et al., 2018) it should be formed during exfoliation or by the coalescence of two itinerant SV(5$|$9) defects. Further, the [SV(5$|$9) + SV(9$|$5)] $\rightarrow$ DV(5$|$8$|$5) merger (Supplemental Material (sup, )) is viewed as one of the relevant mechanisms for multi-vacancy defect formation, and the calculated barrier is found to be 0.68 eV [Table LABEL:tablepvS1 and Figure 4(a)]. The corresponding mechanism is shown in the Figure 4(b). The reverse mechanism of DV(5$|$8$|$5) splitting into two adjacent vacancies requires much higher energy of 1.64 eV. Thus, the itinerant SV defects will easily coalesce to form DV(5$|$8$|$5) and larger defects. In contrast to the present calculations, a much smaller barrier of 1.05 was estimated earlier for the DV(5$|$8$|$5) $\rightarrow$ [SV(5$|$9) + SV(9$|$5)] dissociation. (Cai et al., 2016) The discrepancy arises since the previous estimate is based on the energy differences between the DV and [SV + SV] structures only and without the kinetic consideration of bond reorientation. The split vacancy may also diffuse away, along the armchair and zigzag directions, with activation barriers that decrease with increasing the distance between them. The calculated $E_{a}$ converges to the isolated vacancy migration, while the vacancies are sufficiently separated. sup While compared with the single-layer graphene, the DV defects are known to be stable due to a much high barrier for dissociation above 5 eV. (Lee et al., 2005; Wu et al., 2013)

At present, the proposed mechanism for the anisotropic void formation, that is experimentally observed above 650 K, is the removal of edge atoms from the defects. (Liu et al., 2015b; Yoo et al., 2018; Fortin-Deschênes et al., 2016) However, the corresponding rate-determining atomistic process is not yet known. In this regard, we investigate the transformation of an SV(5$|$9) into the [DV(5$|$8$|$5) + PA] complex, where the under-coordinated P-atom of the SV(5$|$9) migrates to the nearest bridge site [Figure 4(c)]. An activation barrier of 1.46 eV is calculated for this mechanism [Table LABEL:tablepvS1 and Figure 4(c)], which is in excellent agreement with the experimental estimate of 1.64 $\pm$ 0.1 eV. (Fortin-Deschênes et al., 2016) Subsequent removal of the peripheral atoms will require a similar or lower barrier due to weak bonding at the edges and will result in void formation. Further, the PA-adatoms generated in the above process will diffuse along the zigzag direction with extremely low migration barrier of 0.19 eV, which is much higher 1.05 eV along the armchair direction (Table LABEL:tablepvS1).

As the temperature is further increased above 623 K, the degradation process is accelerated, and results in rapid sublimation between 650–680 K. (Liu et al., 2015b; Yoo et al., 2018; Fortin-Deschênes et al., 2016) We address the sublimation process by considering the desorption of single and paired P-adatoms that are produced during the SV(5$|$9) $\rightarrow$ [DV(5$|$8$|$5) + PA] process. The single PA on the pristine surface binds at the bridge site connecting the two armchair edges with 1.43 eV formation energy. The shorter adatom-surface separation of 2.19 Å  highlights the role of lone-pair electrons in the bonding. The desorption mechanism is modelled by vertically detaching the PA by 6 Å  from the surface, which requires overcoming a minimum barrier of 2.14 eV (Table LABEL:tablepvS1).

Highly itinerant PA on the surface forms PA-pairs, and we investigate the simultaneous desorption of both P-adatoms. Two adjacent PA (Supplemental Material (sup, )) along the zigzag directions is found to be thermodynamically most favourable with 2.01 eV formation energy. Thus, the pair formation is thermodynamically advantageous than the two non-interacting PAs, which favours aggregation of adatoms. Further, the desorption of PA-pair is found to be a spontaneous process without an activation barrier, and the free-energy is reduced by 0.82 eV/PA-pair. Therefore, we conclude that the itinerant PAs thermodynamically form PA-pairs, and its concurrent desorption is critical to the rapid sublimation of phosphorene. These results are in excellent agreement with the recent experimental hypothesis of pairwise sublimation. (Fortin-Deschênes et al., 2016)

We argue that the origin of anisotropic void formation is two-fold – directional anisotropy in both vacancy and adatom diffusion. Faster vacancy diffusion along the zigzag direction triggers multi-vacancy defects that are elongated along the zigzag direction. Further, the faster adatom diffusion along the zigzag path produces PA-pairs promoting sublimation. Thus, these processes will create anisotropic voids with the long-axis along the zigzag direction, which is in agreement with the recent experimental observations. (Fortin-Deschênes et al., 2016; Kumar et al., 2018)

III.3 Defect healing and the role of P-adatom

Point-defects can generate larger defects via various mechanisms that are discussed above. However, the point-defects can also be self-healed by interacting with the adatom or interstitial defects, which should be investigated in detail. The vacancy-interstitial (PI) Frenkel pairs with SV(5$|$9) and SV(55$|$66) vacancies [Figure 5(a) and 5(b)] are stable with 1.56 and 1.93 eV formation energies, respectively. In comparison, such Frenkel defects are absent in graphene due to the dominant nature of sp2 bonding, which prevents out-of-plane geometry of the C-adatom within the graphene vacancy. The self-healing of SV(5$|$9) by interstitial PI [Figure 5(b) $\rightarrow$ 5(c)] requires a moderate $E_{a}$ of 0.69 eV. In contrast, the relatively less stable SV(55$|$66) and interstitial pair [Figure 5(a)] first converts into an intermediate [SV(5$|$9) + PI] pair [Figure 5(b)] with 0.1 eV activation energy (Table LABEL:tablepvS1). In contrast, the reverse mechanism of [SV(5$|$9) + PI] Frenkel pair generation at the pristine lattice requires substantially higher energy of 2.25 eV (Table LABEL:tablepvS1).

Having studied the self-healing of single-vacancy, we now investigate the DV(5$|$8$|$5) healing. The partial healing of DV(5$|$8$|$5) to SV(5$|$9) by an itinerant PA requires an activation barrier of 0.46 eV [Table LABEL:tablepvS1 and Figure 4(c)], and this process is accessible at 300 K [Fig. 5(d–e)]. The complete healing of DV then proceeds via the absorption of another mobile PA into the SV with 0.69 eV barrier as discussed earlier, which is the rate-limiting process and occurs at a higher temperature above 450 K. (Liu et al., 2015b) In contrast to phosphorene, this process requires 0.9 eV of energy in graphene and thus, the adatom-DV(5$|$8$|$5) complex is experimentally observed at room-temperature. (Hashimoto et al., 2004; Lusk and Carr, 2008; Wu et al., 2013)

III.4 Oxidation of pristine and defected phosphorene

In addition to the defect-mediated structural degradation, phosphorene also undergoes chemical degradation in the ambient environment. In this regard, we investigate the O2-dissociation mechanism on the pristine and defected phosphorene (Supplemental Material). sup The binding energy is calculated as $E_{b}$ = $E_{\rm(SLP-O_{2})}-E_{\rm SLP}-E_{\rm O_{2}}$, where $E_{\rm SLP}$ is the energy of the pristine SLP, $E_{\rm O_{2}}$ is the energy of the O2-molecule and $E_{\rm(SLP-O_{2})}$ is the energy of the interacting composite system. The physisorbed O2-molecule is very weakly bound ($E_{b}$ = $-$0.09 eV) to the SLP at 3.17 Å  vertical height [Figure 6(a)]. (sup, ) The spin-state of the O2-molecule plays a crucial role in the activated chemisorption and concurrent dissociation. The O2-molecule undergoes a triplet to singlet spin conversion ($S=1$ $\rightarrow$ $S=0$) with increasing proximity to phosphorene. The chemisorption in the singlet state is an activated process with 0.48 eV barrier [Figure 6(a)], which is in good agreement with the previous calculation. (Ziletti et al., 2015) In this chemisorbed state, the O2 is still weakly bound to the SLP ($E_{b}$=$-$0.18 eV), while the O-atoms are weakly bonded to the P-atoms across the armchair direction with an activated O$-$O bond. The ultimate dissociation of the chemisorbed O2 requires 0.13 eV barrier, and the dissociated O-atoms form strong dangling P$-$O bonds with $-$4.22 eV binding energy [Figure 6(a)].

It would be interesting now to investigate the O2-dissociation in the presence of lattice defects. The physisorption at a 2.06 Å  height above the SV(5$|$9) vacancy is still weak with $-$0.21 eV binding energy [Fig. 6(b)]. Similar to the pristine SLP, the spin-state of the O2 is $S=1$. Remarkably, we find that the entire physisorption $\rightarrow$ chemisorption $\rightarrow$ dissociation process is spontaneous without any activation barrier. The present results are in contrast to the previous study, (Kistanov et al., 2017) where a 0.59 eV barrier was predicted due to a different physisorbed configuration, which we find to be a metastable state in the present calculations. The chemisorbed structure has a significantly higher binding energy of $-$1.17 eV compared to the same on the pristine phosphorene due to a strong P$-$O bond that is accompanied by an increased puckering at the SV site [Figure 6(b)]. Unlike the case of pristine SLP, the dissociated O-atoms embed within the lattice and results in a large gain in binding energy. (sup, ) A similar triplet to singlet spin-conversion takes place during the chemisorption at the SV(5$|$9) vacancy. While the O2-dissociation process is similar at the SV(55$|$66) vacancy, the only difference is that the physisorption $\rightarrow$ chemisorption process requires 0.53 eV activation energy. In comparison, we did not find the existence of a chemisorbed state at the DV(5$|$8$|$5), and the direct physisorption $\rightarrow$ dissociation process proceeds with a 0.56 eV barrier [Figure 6(c)]. In all cases, it is understood that once O2 is completely dissociated, it would be difficult to remove from the lattice due to strong P$-$O binding.

III.5 Degradation phase diagram

The defect formation energies and the kinetic pathways discussed so far (Table LABEL:tablepvS1) provide a comprehensive insight into the defect thermodynamics and their concurrent microscopic mechanisms leading to self-healing, oxidation, degradation and sublimation at relevant temperatures. To understand the temperature dependence of the various microscopic mechanisms, we calculate the kinetic rate and temperature phase diagram (Figure 7) using the Arrhenius equation $\Gamma=\nu_{0}\exp(-E_{a}/k_{B}T)$. The pre-factor $\nu_{0}$ is the Debye frequency, which is typically 1013 s*-1*.

At temperatures below 100 K, both the SV(5$|$9) and PA defects are highly itinerant along the zigzag direction, while all other kinetic processes are hindered (Figure 7). However, the anisotropic SV(5$|$9) diffusion in the armchair direction within measurable time-scale occurs only above 150 K. In contrast, the PA diffusion along the armchair direction of the phosphorene lattice is predicted to take place at a much higher temperature above 450 K.

Self-healing of the vacancy defects are rate-limited by the merger of SV(5$|$9) with PI, which appears above 400 K (Figure 7), which is in agreement with the experimental observation. (Liu et al., 2015b) It is important to note that at this temperature range, the itinerant SV(5$|$9) vacancies also merge to form bigger vacancy defects. Thus, we predict a complex interplay between the competing mechanisms leading to self-healing and extended defect formation. However, phosphorization above 400 K will anneal the vacancy defects.

Further, structural degradation will be accelerated beyond 650 K. In addition to much faster SV(5$|$9) mergers, the larger vacancy defects are also created from the SV(5$|$9) point-defects through simultaneous emissions of PA. The generated PAs concurrently form PA-pairs via diffusion and are detached spontaneously as PA-pair from the surface. Thus, above 650 K, the rapid pair-wise sublimation occurs, and the overall degradation mechanism is in agreement with the recent experimental observations. (Liu et al., 2015b; Yoo et al., 2018; Fortin-Deschênes et al., 2016) We also find that below 1000 K temperature, the inter-layer vacancy diffusion does not play any role in degradation in few-layer phosphorene. The chemical degradation via oxidation takes place above 250 K (Figure 7), and after that, it is impossible to remove the dissociated O-atoms from the lattice.

IV Summary

We investigate the structural and chemical degradation of phosphorene within the first-principles calculations. A cohort of microscopic mechanisms is studied to develop a degradation phase diagram. The vacancy diffusion is easily accessible below the room-temperature, which leads to their merger into larger vacancy defects above 400 K. The bond-rotation mechanism of two neighbouring SV(5$|$9) to form DV(5$|$8$|$5) is found to be the rate-limiting mechanism. The self-healing of vacancy via itinerant adatom absorption is also triggered at a similar temperature range. In addition to the merger of mobile vacancies, the emission of the undercoordinated P-atom from the point-defects also generates two-dimensional anisotropic voids above 650 K. Such P-atoms are highly itinerant and thermodynamically form PA$-$PA-dimer, and further degradation proceeds through the spontaneous pair-wise PA sublimation. While the inter-layer vacancy diffusion in few-layer phosphorene is mostly blocked owing to the high activation barrier, the merger of sub-surface vacancies is much faster at a given temperature compared to the same at the surface. However, the degradation through P-emission from the point-defect and the competing self-healing through the P-diffusion are understandably blocked at the sub-surface level. Therefore, the degradation activity including oxidation primarily occurs at the surface.

The chemical degradation that is observed in ambient condition proceeds via O2-dissociation, which we find to be accessible at the room-temperature on the pristine surface, and notably the process is spontaneous at the single-vacancy site. Further, it is impossible to remove the dissociated O-atoms from the lattice, and the most desirable electronic properties of phosphorene are permanently lost. Thus, the present work provides with the microscopic insights into the phosphorene degradation, that will have similarities with the other puckered layered two-dimensional materials.

V Acknowledgements

M.K. acknowledges funding from the Department of Science and Technology through Nano Mission project SR/NM/TP-13/2016, and the Science and Engineering Research Board through EMR/2016/006458 grant. The computational facilities at the Centre for Development of Advanced Computing, Pune; and at the Inter University Accelerator Center, Delhi; and at the Center for Computational Materials Science, Institute of Materials Research, Tohoku University are greatly acknowledged.