Metastable Interlayer Frenkel Pair Defects by Dipole-like Strain Fields for Dimensional Distortion in Black Phosphorus

Herein, our focus lies on close-proximity I-V configurations, in which vacancy-induced two-coordinated atoms strongly interact with vicinal self-interstitials. Particularly, we can classify I-V pair defects of this nature into two broad categories: Type-I and Type-II, depending on the relative location of the vacancy with respect to the self-interstitial. The configurational details (top view) of bilayer phosphorene pertaining to each of these categories are schematically illustrated in Fig. 1(b), with only key topographic features depicted for visual clarity. As shown, Type-I I-V pairs consist of V_in and I_in point defects, and present two symmetry-inequivalent armchair-bridge sites, B₁ and B₂, for self-interstitials. Conversely, Type-II I-V pairs entail V_out and I_in point defects, and host one unique bonding site, B₃, for self-interstitials. All three defect configurations B₁, B₂ and B₃ are determined to be metastable at 0 K, showing no signs of spontaneous decomposition or recombination into its pristine state. Their energetics in both the monolayer and bilayer are presented in Fig. 2(a).

The formation energy of I-V pairs is found to be minimized in the B₂ configuration and can be as low as 1.54 eV on the free surface of the monolayer. While a recent theoretical study by Gaberle and Shluger has contemplated the prospect of Frenkel defects in phosphorene, substantially higher E_f of $\sim$3 eV had been predicted,Gaberle and Shluger (2018) possibly without fully accounting for interdefect coupling. Furthermore, our results indicate that the coalescence between self-interstitials and vacancies is energetically favoured in the monolayer to form B₂ I-V pairs. However, in the bilayer and beyond, intimate I-V pairs tend to possess much higher E_f than their fundamental point defects. This implies a higher concentration of isolated defects than that of Frenkel pairs under equilibrium conditions, which explains the general lack of evidence of Frenkel defects in practical samples of black phosphorus following standard growth processes.Kiraly et al. (2017) Instead, Frenkel defects could be generated under high-temperature non-equilibrium conditions, such as during electron irradiation whereby the lattice is exposed to significant knock-on bombardment events,Vierimaa et al. (2016) or via high-temperature annealing in vacuum.Liu et al. (2015) Appropriate thermal annealing would also help promote the stable formation of Frenkel pairs due to the exothermic process of coupling between isolated defects, as reflected by their positive interaction energies E_i (see Fig. 2(a)).

The interaction energy is at its highest for the B₂ structure, with values of up to 1.68 eV and 1.48 eV in the monolayer and bilayer, respectively. For a deeper insight into the differences between B₁, B₂ and B₃ I-V pairs in the bilayer, we examine the local interlayer distance at the interstitial site d_loc, as shown in Fig. 2(b). This geometric parameter is defined by d_loc = ($(\textit{z}_{1}+\textit{z}_{2})-(\textit{z}_{1^{\prime}}+\textit{z}_{2^{\prime}})$)/2, where z_i is the z-coordinate of the neighbouring atom i to the self-interstitial (see given schematic in Fig. 2(b) for the reference case I_in). The values of d_loc for B₁ (4.10 Å) and B₃ (4.23 Å) are observed to be much larger than that of I_in (3.91 Å), while that for B₂ (3.82 Å) is relatively small in comparison. This suggests that the stabilising effects of interlayer binding are the strongest in the B₂ structure such that it displays the highest resistance against dissociative and recombinative mechanisms, especially under elevated temperatures. Hence, in agreement with earlier analysis, the B₂ I-V pair is the most energetically preferred adaptation of the Frenkel defect. Moreover, by performing finite temperature (canonical ensemble) ab initio molecular dynamics simulations at 300 K and for a period of 10 ps, we establish that the B₂ I-V pair can remain thermodynamically stable at room temperature (see Section II of the SM,Sup and referencesTheobald (2005); Welford (1962) therein).

In Fig. 2(c), we consider the effect of net charge on the formation energy of the B₂ I-V pair in monolayer, bilayer and bulk phosphorene. For all cases, both negatively ($-$1) and positively ($+$1) charged states are energetically unstable as compared to the defect in its neutral (0) condition. This implies that, in the near proximity for forming a Frenkel pair, deviation of the phosphorus atom from its equilibrium site tends not to alter the electrostatic environment of the host lattice. Such neutral tendency is reasonable for a non-polar homoelemental sheet like phosphorene as the Frenkel pair defect can remain pseudo self-compensated with weak charge accumulation/donation.

In this section, the effect of interdefect coupling on the electronic band structure of bilayer phosphorene is investigated. Spin-polarized band structures are calculated in Fig. 3(a) along the high symmetry path $\Gamma$-X-S-$\Gamma$-Y-S defined for a 2D orthorhombic lattice. Note that while the optB88-vdW functional tends to underestimate the band gap (see Table S2 in the SM,Sup and referencesCai et al. (2014); Dai and Zeng (2014); Perdew et al. (1996); Heyd et al. (2003, 2006); Qiao et al. (2014); Liu et al. (2014); Keyes (1953); Warschauer (1963); Maruyama et al. (1981); Narita et al. (1983) therein), it remains sufficiently reliable for making qualitative predictions for the purposes of our study.Cai et al. (2014); Pan et al. (2017) In its pristine state, bilayer phosphorene is a non-magnetic semiconductor with a direct band gap (0.46 eV) at the $\Gamma$ point. The presence of an inner vacancy (V_in) is shown to introduce defect states that intersect the Fermi level, reflective of a metallic system. On the other hand, for the case of an inner self-interstitial (I_in), the band structure contains two defect bands which straddle the Fermi level. Both V_in and I_in are spin-polarized with net magnetic moments of 0.39 $\mu$_B and 0.97 $\mu$_B, respectively. Interestingly, interdefect coupling in the form of the B₂ I-V pair is found to effectively quench the net magnetic moments originating from the fundamental point defects, while widening the direct band gap at the $\Gamma$ point to 0.60 eV (blue-shift).

The underlying physical mechanism behind the widening of the band gap is ascribed to the breaking of lattice periodicity in the vicinity of the Frenkel pair, through (i) a (5$\mid$9) reconstruction at the vacancy site, and (ii) a localized increase in the interlayer distance of the bilayer structure such that it essentially behaves as two decoupled monolayers in the defect neighbourhood. As discussed previously (see Fig. 2(b)), the local interlayer distance associated with the B₂ I-V pair can be as high as 3.82 Å at the self-interstitial site, significantly exceeding the equilibrium separation distance of 3.20 Å in the pristine bilayer. This mechanism is further corroborated by an analysis of the projected density of states, as shown in Fig. 3(b). Notably, we see that the regions surrounding the valence band maximum and conduction band minimum are primarily dominated by p_z orbital states (red), which illustrate the strong role of interlayer interactions on the electronic properties of black phosphorus. In this regard, the widening of the band gap is attributed to a net shift of the valence band edge towards lower energy levels (black arrow), while that of the conduction band remains relatively fixed with respect to the Fermi level.

The local atomistic structure of the B₂ I-V pair is analysed using discrete geometry within the context of 2D material nets.Lee (1997); Bobenko and Suris (2008) As opposed to conventional continuum-based methods for parameterizing the strain field, discrete geometry goes beyond first order elasticity to provide an exact description of shape under arbitrary deformations. The discrete geometry of a surface is characterized by four invariants of its metric tensor g and curvature tensor k, denoted by the trace and determinant of g (i.e. Tr(g) and Det(g)), mean curvature H and Gaussian curvature K.

The discrete geometry is formulated using the method of triangulations over a finite mesh of atomic positions,Bobenko and Suris (2008); Weischedel et al. (2012); Mehboudi et al. (2015) as illustrated in Fig. 4(a). Each triangulated element comprises of three adjacent atoms (denoted by atom index p = 1, 2, 3), and is defined by the directed edges s₁, s₂ and s₃ such that $\textbf{s}_{1}+\textbf{s}_{2}+\textbf{s}_{3}=0$. A discrete analogue of the infinitesimal length ds² is given by ${\textit{Z}_{p}}^{\text{I}}=\textbf{s}_{p}\cdot\textbf{s}_{p}$, which yields the square of the shortest distance between atoms. The variation in orientation between the normal vectors $\textbf{m}_{q}$ and $\textbf{m}_{r}$ is projected onto their common edge s_p following the relation ${\textit{Z}_{p}}^{\text{II}}=(\textbf{m}_{r}-\textbf{m}_{q})\cdot\textbf{s}_{p}$. Here, $\textbf{m}_{p}$ is the mean of the normal vectors of triangulated elements which immediately surround atom p (see dashed outlines). Another relevant quantity, known as the dual edge, is given by $\textbf{s}_{p}^{*}=\textbf{s}_{p}\times\textbf{u}$, where u is the normal vector of the triangulation described by atoms 1, 2 and 3. The area of the triangulated element in its pristine (reference) condition is denoted by A₀, while that under arbitrary deformation is indicated by A₁. The metric tensor g and curvature tensor k can then be expressed as follows:

The resulting tensors g and k are in the form of 3 $\times$ 3 matrices, with the parentheses (p, q, r) representing the set of three contributions (1, 2, 3), (2, 3, 1) and (3, 1, 2) included in the summation. The eigenvalues {k₁, k₂, 0} of the curvature tensor k provide the principal curvatures k₁ and k₂ for each triangulation such that $\textit{H}=(\textit{k}_{1}+\textit{k}_{2})/2$ and $\textit{K}=\textit{k}_{1}\textit{k}_{2}$ can be evaluated. The four invariants Tr(g), Det(g), H and K corresponding to each atomic position are averaged over their respective values at all triangulated elements that share the same vertex. As shown in Figs. 4(b)-(c), the reference values of these invariants are Tr(g) = 1, Det(g) = 1, H = 0 and K = 0, which correspond to the case of pristine phosphorene in its strain-free, planar configuration. Tr(g) and Det(g) take on values greater (less) than 1 for varying degrees of in-plane tension (compression), while H and K, in the usual mathematical sense, measure the out-of-plane deviation from planarity and character of the surface profile (i.e. elliptical (K $>$ 0) or hyperbolic (K $<$ 0)). Due to the puckered structure of phosphorene, the normalized layer thickness $\overline{t}$ provides another quantity of interest for our analysis.

The discrete geometry of the B₂ I-V pair is evaluated for the inner and outer sub-lattice of each layer individually, as shown in Fig. 5. Note that these results correspond to the case of the defect in the dilute limit and are independent of supercell size (see Fig. S3 in the SM).Sup From Tr(g) and Det(g), we find that the vacancy-induced (5$\mid$9) reconstruction generates significant tensile strain in S${}^{1}{\textsubscript{in}}$, while the encapsulation of the self-interstitial yields low-to-moderate compression (tension) in the inner (outer) sub-lattice of both layers. The resulting strain distribution consisting of local tensile ($\epsilon+$) and compressive ($\epsilon-$) zones creates a dipole-like mechanism which drives the coupling and binding of the Frenkel pair. Since the strain fields imposed by the vacancy and self-interstitial are highly anisotropic, the coupling strength of Frenkel defects will strongly depend on the relative arrangement of interstitials and vacancy sites in the lattice.

Another plausible factor that accounts for the binding of Frenkel defects is related to the electronic interaction. It is noted that both isolated vacancies and self-interstitials are accompanied with spin-polarized dangling states (see Fig. 3(a)). However, upon formation of the Frenkel pair, such shallow empty defect states are saturated and a new hybridized state is formed and buried within the valence band; accordingly, the I-V pair is stabilized. In addition, the coupling between the vacancy and self-interstitial gives rise to regions of strong, I-V oriented rippling in both the inner and outer sub-lattices of layer 1, as reflected by H and K. The structural signature of the B₂ I-V pair is also shaped by changes to the layer thickness, which decreases by up to 20% in the neighbourhood of the vacancy in layer 1. The distinct out-of-plane dimensional distortion introduced by such defects can enhance the chemical activity of locally rippled regions around vacancy sites.Utt et al. (2015); Kistanov et al. (2016) Interestingly, this may serve to facilitate the substitutional doping process during ion implantation treatments of phosphorene-based semiconductors.

By virtue of the close proximity between vacancies and self-interstitials in the B₂ I-V pair, short-range recombination can take place. In phosphorene, recombination events may also be driven by the inherent presence of highly mobile vacancies, which are particularly faster along the zigzag direction.Cai et al. (2016) Here, two paths to recombination - direct (red) and indirect (blue) - are preliminarily proposed, as indicated in Fig. 6(a) (top panel). In the direct route, the migrating self-interstitial simply recombines at the inner vacancy site of layer 1. On the other hand, in the indirect route, both the vacancy and self-interstitial undergo simultaneous migration towards one another such that recombination takes place at the original site of the two-coordinated host atom.

The variation in bond length d_m between the self-interstitial and two-coordinated host atom during migration is monitored in Fig. 6(a) (bottom panel). Region I (shaded grey) denotes the range of d_m which is less than the atomic radius of P (0.95 Å), and is therefore highly unphysical. Clearly, the direct route appears to be energetically unfavourable owing to the excessive overlap of atomic cores for a large part of the recombination process. In contrast, the indirect route presents a more viable mechanism as the bond length d_m is well-preserved throughout the course of migration. The energy profiles corresponding to the indirect path to recombination are shown in Fig. 6(b), while the transition state (TS) is given in the inset. We determine the activation energy E_a to be 1.46 eV in the recovery of the pristine bilayer, with a net energy release of $\sim$2.4 eV. By thinning down phosphorene even further into its monolayer, the recombination barrier increases to 1.77 eV due to quantum size effects that enhance the strength of P-P bonding. The recombination barrier of phosphorene-based I-V pairs is found to be higher than that reported for atomically-flat graphite (E_a = 1.30 eV) which is responsible for Wigner energy release at 470 K.Mitchell and Taylor (1965); Ewels et al. (2003) This is surprising considering that graphitic C-C bonds are generally much stronger than that of P-P in phosphorene, and goes to show the pivotal role of lattice curvature in the formation of highly metastable Frenkel defects.

Our results are consistent with recent reports on the formation of skeletal structures of black phosphorus annealed at around 670 K.Liu et al. (2015) The nature of such dimensional distortion has long been open to speculation as an amorphous red phosphorus-like skeleton. Here, we suggest that this structure could be closely linked with the onset and development of intimate interlayer I-V pair defects. While not been studied, these complex defects could even lead to higher-order defect rearrangements and spontaneous release of energy with either electron irradiation or thermal treatment.

On closer inspection of the intermediate states during recombination, we are able to obtain a clearer understanding of how the band gap might evolve with structure (see Fig. 6(c)). Lattice instabilities and distortion, which tend to introduce defect states near the Fermi level, ultimately lead to minimum band gap values at the TS for both the monolayer (0.42 eV) and bilayer (0.16 eV). The paths connecting the TS to either the I-V pair reconstruction or pristine lattice seem to provide alternative modes by which atomic orbitals can interact and hybridize, thereby giving rise to entirely different band gap widths at equilibrium in each case. Through verification of the band structure, we can also attest that in low concentrations, the presence of Frenkel defects should not alter the direct-gap character of phosphorene at the $\Gamma$ point.