Phosphorene under electron beam: from monolayer to one-dimensional chains

Powered by TCPDF (www.tcpdf.org) This material is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Vierimaa, Ville; Krasheninnikov, Arkady V.; Komsa, Hannu Pekka


Introduction
Exfoliation of graphene from a graphite crystal 1 indicated that other two-dimensional (2D) systems can be manufactured, provided that the "parent" bulk material has a layered structure with the layers being bound together by weak van der Waals forces. Indeed, a large number of 2D materials such as hexagonal boron-nitride (h-BN), 2 transition metal dichalcogenides (TMDs) 3,4 and phosphorene 5,6 have been produced from bulk counterparts by mechanical and liquid exfoliation techniques. The latter material derived from black phosphorus has recently attracted a lot of attention, as it is mechanically flexible, 7 has high electron mobility 8 and exhibits a gap of about 1-2 eV in the electronic spectrum, which is tunable by strain and the number of layers. [9][10][11] 2D materials have extensively been characterized by transmission electron microscopy (TEM), and their morphology ( just a few layers of atoms) together with rapid progress in TEM and related spectroscopy methods made it possible to not only visualize directly the atomic structure, but also follow the motion of atoms and atomic structure transformations in real time with atomic resolution. Besides, beams of energetic electrons, which can nowadays be focused into a sub-Ångström area, have been shown to work as cutting and welding tools on the nanoscale. 12 Specifically, nm-wide ribbons have been manufactured by the "top-down" approach, that is by electronirradiation-induced sputtering of atoms, from sheets of graphene, 13 h-BN, 14 and TMDs. [15][16][17] Moreover, as an ultimately narrow ribbon, free-standing atomic chains have been made starting from graphene 13,18,19 and h-BN, 14 which can serve as interconnects between 2D materials, 20 along with beamcreated nanowires in TMDs. 16,17 Few-layer phosphorene sheets have been characterized by TEM. 6,[21][22][23] Electron-beam engineering of nanostructures based on a phosphorene sheet can be of particular interest, as defects in this system have been predicted to be electronically inactive, 24 thus limiting the undesirable effects of unavoidable radiation damage of the system. Likewise, phosphorus chains, if stable, may exhibit intriguing electronic properties and unusual Peierls distortions due to fractional band filling. However, the efficient use of the electron beam for shaping and tuning phosphorene structures is virtually impossible without fundamental understanding of the response of this system to electron irradiation. At the same time, while the behavior of many other 2D materialsgraphene, [25][26][27][28] TMDs, [29][30][31] and h-BN 17,32-34 under electron beam has been studied at length, in phosphorene very little is known about the defect production mechanisms, defect dynamics, and irradiation damage development.
Here, using first-principles simulations, we study the production of defects in a single phosphorene sheet under electron beam. By employing the McKinley-Feshbach formalism and accounting for the thermal motion of atoms, we assess the cross section for atom displacement as a function of electron energy. We further investigate the energetics and dynamics of point defects and the stability of ribbons and edges under electron beam. Finally, we show that P atomic chains should be surprisingly stable, and their atomic structure is not linear giving rise to the absence of a gap in the electronic spectrum.

Methods
All simulations were based on density-functional theory (DFT) and carried out using the Vienna ab initio simulation package (VASP) 35,36 with the PBE functional. 37 The atomic structure of the pristine system was calculated using the plane-wave basis with a cut-off energy of 400 eV and a 12 × 12 k-point mesh.
The knock-on radiation damage was assessed by calculating the displacement threshold energy T d , which corresponds to the minimum initial kinetic energy needed to displace a phosphorus atom from its original lattice site without its immediate recombination with the associated vacancy. Assuming an instantaneous momentum transfer from the impinging relativistic electron to the atom during the impact (it happens at the attosecond time scale 38 ), an initial velocity was assigned to a selected P atom and the time evolution of the system was simulated using DFT molecular dynamics (MD). The initial kinetic energy was varied to trace the thresholds with an energy resolution of 0.5 eV. We used a time step of 0.5 fs, and a run time of 0.5-1 ps was usually sufficient to find out if the initial kinetic energy was high enough to produce a defect. The effect of spin-polarization (typically affecting the total energy by less than 0.1 eV) was found to have a negligible effect on the displacement thresholds at the adopted energy resolution, and thus calculations are carried out without spinpolarization. These calculations, as well as all calculations involving defects, were done using a 6 × 4 supercell. The Brillouin zone was sampled with a 2 × 3 Γ-centered k-point mesh.
Before knocking on another atom, the geometry of the structures with defects was fully optimized. At the same time, we calculated the formation energies of the produced defects: where E(defect) is the total energy of the system with defects composed of N P phosphorus atoms, and the chemical potential of phosphorus μ P is taken from the pristine phosphorene. Formation energies were calculated with spin-polarization. Barriers for defect migration were obtained using the nudged elastic band method. 39 The cross sections for atom displacements upon impacts of energetic relativistic electrons were assessed using the McKinley-Feshbach formalism which accounts for the thermal motion of ions. 26,40 3. Results

Pristine systems
The atomic structures of several phosphorus-based systems are presented in Fig. 1. Fig. 1(a) shows the top and side views of the pristine phosphorene lattice, highlighting its puckered character. For the pristine system, we obtained lattice constants of a = 3.30 Å in the x-direction and b = 4.62 Å in the y-direction, in agreement with the values reported in the literature. 7,41 For the edges in the zigzag direction, there are two inequivalent positions to cut the structure: in the middle of a terrace or in the middle of a cliff. Consequently, the system adopts two reconstructions, the zigzag (ZZ) and reconstructed cliff (RC), which, after structural relaxation become as shown in Fig. 1(c and d). The nanoribbon formation energy is defined as E NR f = E(ribbon) − N P μ P , where E(ribbon) is the total energy of the ribbon composed of N P phosphorus atoms and μ P is defined above. E NR f is plotted in Fig. 2 as a function of ribbon width for different types of the edges (two edges). Using the values for the widest ribbons considered, we obtain the edge energies of 0.70 eV and 0.52 eV per edge atom (or 0.21 eV Å −1 and 0.16 eV Å −1 ) for the ZZ and RC edges, respectively, in good agreement with those reported in ref. 42 about the electronic, optical, and mechanical properties of phosphorene edges and wider nanoribbons can be found in ref. 43-47. In the case of very thin ribbons, the system may relax strongly, which considerably lowers E NR f . The lowest energy structures are illustrated in Fig. 2. By running MD at 300 K, we ascertained that all the theoretical structures we found are at least metastable.

. Further information
Ultimately, the thin ribbons turn into chains. The structures of pristine atomic chains were studied in more detail in linear, zigzag, and armchair configurations. The optimized structures for zigzag and armchair chains are shown in Fig. 1(b) and the structural parameters are listed in Table 1, along with cohesive energies. Moreover, phosphorus has a tendency to form P 4 -molecules, and although phosphorene is in the stable phase, both the bond length and cohesive energy of the P 4 molecule and phosphorene are very similar.
Calculations for chains showed that the linear configuration is unstable, and the system forms either a zigzag or an armchair structure, as shown in Fig. 1(b). The optimal configuration turned out to be a zigzag pattern, where the bond angle is 95°and bond length is 2.12 Å. The cohesive energy for zigzag chains is 3.03 eV, a value much larger than that for the simple linear chain. For the armchair pattern, there are two inequivalent bonds: the bond length between atoms in the direction of the chain (x-direction) is 2.08 Å, while in the transverse direction the length is 2.17 Å. The bond angle is 107°a nd the cohesive energy is 2.93 eV, slightly smaller than that for the zigzag chains. As compared to carbon chains, 19 where the linear configuration is lower in energy because the sp 1 -hybridization naturally originates from the distribution of 4 electrons over the two sp 1 and two transverse π-type bonds, the extra electron in phosphorus results in non-linear configurations.

Displacement threshold calculations
Having studied the atomic structure of phosphorene sheets, ribbons and chains, we move on to assess the radiation tolerance of these systems under electron beam using DFT MD. We stress that here we study only the knock-on damage. We will discuss other channels for defect production such as electronic excitation and electron-beam-induced chemical etching later on. As T d can be correlated with the formation energies of vacancy-type defects, 29 we also studied the atomic structure and energetics of point defects.
3.2.1. Monolayer phosphorene. The atomic structures and displacement thresholds are presented in Fig. 3 for the pristine phosphorene and for the defects important during the first steps of the irradiation damage development. Fig. 3 also shows formation energies for defects and summarizes the damage pathways along with the corresponding displacement thresholds.
Starting from pristine phosphorene, no damage is observed for T < T d ∼6 eV. When T ≥ 6 eV, the P 1 atom has enough kinetic energy to break the bond to P 2 . The P 3 atom then moves in to fill the vacant site. Finally, as P 1 rebounds back to the original site, it forms a bond with P 3 instead of P 2 . The final structure is shown in Fig. 3(b), and it effectively consists of a local bond rotation leading to 4-fold and 2-fold coordinated P atoms in contrast to the typical 3-fold coordination. In accordance with the bond rotation of carbon atoms in graphene, this defect was also named Stone-Wales defect (SW) in ref. 24. Larger kinetic energies lead to displacement of P 1 to adatom sites or sputtering of P atoms.
From the SW defect the P 1 atom can be displaced to a more distant adatom site with T ≥ 4.5 eV, where it may easily diffuse (see section 4), and finally [ Fig. 3(d)] can be sputtered from the monolayer with T ≥ 4.5 eV. Thus, in order to start the process from pristine phosphorene, a minimum 6 eV kinetic energy is needed, but then irradiation would likely lead to sputtering of phosphorus atoms or possibly to formation of more complicated atomic structures.
Curiously, knock-on of the P 2 atom from pristine phosphorene with T ≥ 7.5 eV leads to the formation of a vacancyadatom complex (V P − P ad ), which is shown Fig. 3(c). It is similar to the SW defect, but P 1 is on a closeby adatom position and both P 3 atoms around P 1 have symmetrically moved in to bond with P 2 . Further knock-on events turn it into SW or lead to direct sputtering of the P 1 atom, and thus does not change the general picture. Depending on the magnitude and direction of the collision, both of these should occur, but as evident from the formation energies shown in Fig. 3, the SW defect is clearly more stable. On a related note, we occasionally  observed also other events. For instance, at T 1 = 9.5 eV, P 1 and P 3 atoms swapped places (full bond rotation) and thus resulted in the pristine lattice. The T 2 = 10 eV event first appeared to lead to the formation of the (V P − P ad ) complex, but the P 1 atom bounced back from the adatom site and the system returned to the pristine lattice. The sputtering of phosphorus or displacing them to distant adatom sites leads to the formation of single vacancies. The P atoms next to the vacancy are likely less tightly bound than those in the pristine lattice and thus easier to displace. In order to manage the computational workload, we only calculated T d for selected atoms around the vacancy. Nevertheless, this is sufficient to obtain a general picture of how easily the vacancy may expand.
The results are shown in Fig. 4. A single vacancy with formation energy E f = 1.67 eV is relatively unstable when compared to a double vacancy or larger vacancies, all having formation energies of about 0.7 eV per missing P atom. This is reflected also in the displacement thresholds. The P 1 atom in V P is very easily displaced with T 1 ≥ 3.5 eV, whereas the P 2 atom in V 2P requires T 2 ≥ 6 eV and P 3 atoms in V 3P require T 3 ≥ 5 eV.
Our formation energies are generally in agreement with those reported in previous studies, 24,48 which can be referred to for discussion on the electronic structure of these defects. However, it is worth commenting on a few differences in our work. Our formation energy for V 2P is in agreement with that in ref. 24, but is 0.5 eV lower than the rather similar-looking double vacancy in ref. 48. The slightly asymmetric structure given in ref. 48 can be obtained by displacing of one of the atoms next to the P 2 atom to the bottom layer in Fig. 4(b). Hu and Yang also found another SW configuration that is 0.3 eV lower in energy than the one shown in Fig. 3(b), where the bond rotation involves P 1 and P 2 instead of P 1 and P 3 . We did not observe this configuration in our molecular dynamics studies, but their formation is certainly possible. Nevertheless, this indicates that amorphization of phosphorene by bond rotations may be possible under electron irradiation similar to what is observed in graphene. 49 Finally, we found that the P ad configuration shown in Fig. 3(d) on top of the bond-center is favored over the position reported in ref. 24. All of the above indicate that due to its flexibility phosphorene may possess large number of defect structures corresponding to local energy minima.
We found that similar to graphene, 25 and unlike transition metal dichalcogenides, 29 calculating unrelaxed vacancy formation energies does not give a good estimate for the displacement threshold, but requires explicit MD calculations of T d . Physically, this is the corollary of equal masses of all atoms in the system, which gives rise to the redistribution of the initial kinetic energy of the recoil atom: its neighbors "follow" it, as opposed to TMDs, where chalcogen atoms are considerably lighter than the neighboring transition metal atoms.
3.2.2. Nanoribbons and edge stability. We next consider the irradiation stability of phosphorene edges that may be present either simply at the boundaries of the sample, or appearing during exposure to the beam due to agglomeration of vacancies or enlargement of voids. As evident from Fig. 4(d), the V 4P structure already looks like being formed out of two zigzag edges.
Here, we only determine the displacement thresholds, but do not attempt to trace further how the damage might progress. Although in practice the edge may be quite amorphous and expected to damage more quickly, our approach will still yield limits to the sputtering rate. In these calculations, to reduce the effects arising from the finite width of the ribbon, the non-damaged edge is passivated by hydrogen. 9,43,45 The supercell is 6 units wide along the edge of the ribbon.
The results of DFT MD calculations are shown in Fig. 5. Starting from the case of the reconstructed cliff edge, which has lower edge energy, the outermost P atoms can be sputtered with energy of 3.5 eV or 4 eV, depending on the direction. For the second outermost row of P atoms, we found that T d is between 5 and 6 eV, very close to the values in the pristine sheet. The collisions mostly led to sputtering of P 2 dimers from the edge, except for the T 2,in case where a whole P 5 ring was sputtered from the edge. Nevertheless, all cases resulted in formation of a zigzag edge. In the case of the zigzag edge, the outermost P atom can be displaced with an energy of 4 eV or 7 eV, depending on the direction, and usually instead of sputtering atoms this led to local reconstructions similar to those in the pristine systems. For three of the four cases involving outermost P atoms, T d is 3.5-4 eV and thus the edge should get easily destroyed. Only in the T 3,out case the edge is expected to withstand radiation better and thus might be present more often during imaging.
3.2.3. Chain stability. Formation energies listed in Fig. 2 indicate that the atomic chain is more stable ( per atom) than thin nanoribbons, suggesting that they might also be stable under electron beam.  When the collision between an electron and the zigzag chain is simulated, there are three relevant non-equivalent cases to be considered: displacement of the two inequivalent atoms in the y-direction and displacement of one of the atoms in the z-direction [cf. Fig. 1(b)]. The simulations were made with two different supercells, the smaller one containing 12 atoms and the larger one 24 atoms. The smaller system was first used to find an estimate for threshold energy for chain breaking, but a larger supercell was needed to obtain quantitatively accurate results.
It turned out that when the initial velocity vector points away from the chain (the y-direction), the threshold energy is the highest, being around T d = 8 eV. A slightly surprising result is that the second y-direction and the z-direction have the same threshold energy of T d = 7 eV. At the threshold, the kinetic energy received by an atom is not enough to make it escape, but the chain will break into two parts. Below the threshold, the chains were strongly bent and the zigzag pattern became distorted. Nevertheless, the chains were able to stay intact for the whole simulation time of 900 fs. Moreover, one cannot exclude reconnection of the chains on a microscopic time-scale.
To gain further insight into stability of the chains, we also carried out molecular dynamics simulations at high temperatures. The initial configuration was chosen to be zigzag, as it has the lowest cohesive energy. The temperature in the simulations was kept constant using the Nose-Hoover thermostat. The time-step was chosen to be 1 fs.
Starting from T = 300 K, the system was simulated for a few picoseconds at time. If the chain was still intact, the temperature was increased and the simulation was continued from the final configuration of the previous simulation. The second simulation was done at T = 600 K and after that the temperature was increased by 200 K between simulations. The chains turned out to be surprisingly stable. As the temperature increases, the chains start to rotate around their axes while keeping the zigzag-pattern intact at temperatures up to 1400 K. After this, the chain breaks up into P 2 and P 4 molecules. The ability to rotate seems to be the main reason for the chains' stability. It allows the atoms to have much more kinetic energy without escaping and breaking the chain. At temperatures of 600 K and above, the zigzag chains occasionally showed formation of armchair segments, but lasting for only a very short time. In simulations for the armchair chain, the system transformed back to the zigzag configuration within 1-2 ps, as expected on the basis of their cohesive energies ( Table 1).
The effects of strain, which may be present and considerably affecting the electronic properties, 50 were studied by elongating the unit cell in the x-direction. The percentage of elongation directly translates to strain and MD simulation allows us to find out the maximal amount the chain can sustain. The critical strain, after which the splitting happens, was estimated to be 13%. Also, a curious observation from the splitting is that the two atoms closest to the splitting point tend to rotate around the axis chains, terminating the chain ends with a P 4 -like configuration.

Cross sections
Using the values of displacement threshold energies calculated using DFT MD, we evaluated the cross sections for atom displacement. The knock-on damage cross sections were calculated at room temperature and are plotted in Fig. 6. The curve for T d = 6 eV represents the threshold for initiating damage in pristine phosphorene. An electron beam with 80 keV already has rather high cross section of 25 barns. Moreover, since the defects can grow rapidly, we conclude that 80 keV will likely quickly induce excessive damage. The energy of 60 keV should limit the damage of the pristine sheet, but will result in significant damage at the edges and defect sites. Note however that the divacancy (V 2P , T d = 6 eV) should be stable under 60 keV. According to our calculations, most of the knock-on damage can be avoided when using a 30 keV electron beam.
We note that our calculations only yield a lower estimate for the knock-on displacement, but other processes such as electronic excitations, charging, and chemical reactions 51,52 may also play a great role or even dominate. For instance, it was found that imaging at 15 kV or 30 kV leads to damaging of BN sheets, although their knock-on displacement thresholds are calculated to be clearly higher. 14 However, charging and electronic excitation effects should be less important in phosphorene as it has a much smaller band gap, while etching contribution to the damage may be quite substantial due to high reactivity of phosphorene as evidenced e.g. by its strong propensity for oxidation. 53,54

Defect dynamics
The molecular dynamics studies of knock-on events also hinted of efficient dynamics of defects, which include their migration and transformations between different types of defects. We first consider the barriers for the transformations Fig. 6 Calculated cross sections for knock-on damage induced by electron beam at selected displacement threshold energies and electron energies ranging from 0 to 300 kV. The shaded region indicates the variation in cross sections when moving from the static lattice to the Maxwell-Boltzmann velocity distribution at 300 K. between defect configurations that conserve the total number of atoms in the system, i.e., the Stone-Wales and the (V P − P ad ) defects. The energetics of the transformation with respect to pristine phosphorene is shown in Fig. 7. Formation of SW defects requires overcoming a barrier of 2.30 eV. This value is rather low, at least when compared to graphene, 25 indicating that the defects can be formed at high temperatures. On the other hand, healing of the SW defect only requires 0.81 eV and thus even if they are plentiful after e.g. irradiation experiments, they can also be quickly removed. The barriers between (V P − P ad ) and SW are small, but since SW is lower in energy, it is expected to be predominantly present.
Clustering of vacancies will be limited by the rate of vacancy diffusion. Diffusion of a single vacancy consists of two separate steps, which do not involve a clear jump of atoms from one site to another, but rather involve changes in the local bonding, similar to vacancy diffusion in carbon nanotubes. 55 As illustrated by configurations i and ii in Fig. 7(b), there can be a switch in the bonding between the two Jahn-Teller symmetry broken structures, which only involves movement of atoms P 1 and P 2 w.r.t. atom P 3 . The P vacancy remains next to atom P 3 . The energy barrier for this process is 0.25 eV. In addition, as illustrated by configurations ii and iii, atom P 2 can easily switch bonding to either atom P 1 or atom P 4 with only a very small change in its position. Consequently, the vacancy has moved to the neighboring lattice site. The small change in the geometry is reflected in a very small migration barrier, only 0.09 eV. Either of these barriers is low enough to warrant that rapid diffusion should occur even at room temperature.
In the case of adatoms, diffusion along the trench, from one "bond-center" site to another (from i to ii), is easy, with a barrier of only 0.23 eV. The path for migrating from one trench to another is more complicated as shown in Fig. 7(c) (configurations ii to v). The rate is limited by the barrier of 1.19 eVthis energy is required to move from the trench to the middle of the ridge.

Electronic structure of phosphorus chains
The band structure and the density of states for the zigzag chain are shown in Fig. 8. In addition to the total density of states, projected densities of states (PDOS) are also calculated to gain insight into the electronic structure. PDOS allows us to immediately recognize some of the bands to be formed by certain orbitals. The yellow-coloured bands are largely formed out of p z -orbitals forming π-bonds.
The blue-coloured bands originate mostly from s-orbitals with a small p x /p y contribution and the black bands from p x /p yorbitals with a small contribution from s. These bands are associated with the P-P bonds in the chain: the bonds lie in the xy-plane but they are not aligned with either of the major axis and therefore appear in both of the projections. Thus, they can both be considered as σ-bonds, the former arising from s-orbitals and the latter from p x /p y -orbitals. The minor mixing of the orbital characters indicates partial sp 2 -hybridization, which is also evidenced in the 95°angle found in the zigzag chain.
The system always remained (semi-)metallic with a halffilled π-band and the Fermi-level crossing it at the K-point. No indication for the Peierls-type transition was found apparently due to a non-linear structure of the chain.

Conclusions
Defect production in single phosphorene sheet under impacts of energetic electrons was studied using first principles calculations. Molecular dynamics was used to gain insights into time evolution of the system after impacts from electrons and the displacement cross sections were determined through the McKinley-Feshbach formalism. Although we focused on the stability under electron beam, the calculated displacement threshold energies, the most important quantities in radiation physics, are equally valid for ion or neutron irradiation, and the amount of damage in each case can be estimated using the Kinchin-Pease model. 56 In pristine phosphorene, the damage is expected to begin from the formation of Stone-Wales defects, from which P atoms can then be relatively easily displaced, finally leading to vacancy formation. The calculated displacement threshold is low enough that a 80 keV beam should rapidly lead to a considerable amount of damage. With 60 keV electrons, the initial SW production can be avoided. Next, we considered defect dynamics and stability of edges under electron beam. To avoid sputtering of P atoms from defect sites or from edges, the beam energy should preferably be even less than 60 keV. We stress that the knock-on damage mechanism should dominate under ultra-high vacuum conditions to prevent oxidation of phosphorus nanostructures and the associated etching. Phosphorene sheets may be encapsulated between graphene sheets, though, which should protect them from the interaction with the environment, similar to TMDs. 31,57,58 Finally, upon continuing thinning of nanoribbons by the electron beam, P atomic chains are expected to form. Calculations indicated them to be surprisingly stable, both thermally and under electron beam. They adopt zigzag configuration and their atomic structure is not linear giving rise to the absence of a gap in the electronic spectrum.