Dynamics and Mechanisms of Exfoliated Black Phosphorus

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Letter

Dynamics and Mechanisms of Exfoliated Black Phosphorus Sublimation Matthieu Fortin-Deschenes, Pierre L. Levesque, Richard Martel, and Oussama Moutanabbir J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.6b00584 • Publication Date (Web): 11 Apr 2016 Downloaded from http://pubs.acs.org on April 11, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Dynamics and Mechanisms of Exfoliated Black Phosphorus Sublimation Matthieu Fortin-Deschênes,1 Pierre Levesque,2 Richard Martel,2 Oussama Moutanabbir1,* 1

Department of Engineering Physics, École Polytechnique de Montréal, C. P. 6079, Succ. Centre-Ville, Montréal, Québec H3C 3A7, Canada 2 Département de Chimie, Université de Montréal, 2900 boulevard Edouard Montpetit, Montréal, Canada *Corresponding author. Email: [email protected]; Phone: +1-514-340-4711.

Abstract. We report on real time observations of the sublimation of exfoliated black phosphorus layers throughout annealing using in situ low energy electron microscopy. We found that sublimation manifests itself above 375±20°C through the nucleation and expansion of asymmetric, faceted holes with the long axis aligned along the [100] direction and sharp tips defined by edges consisting of alternating (10) and (11) steps. This thermally activated process repeats itself via successive sublimation of individual layers.

Calculations and simulations using density

functional theory and kinetic Monte Carlo allowed to determine the involved atomic pathways. Sublimation is found to occur via detachments of phosphorus dimers rather than single atoms. This behavior and the role of defects is described using an analytical model that captures all essential features. This work establishes an atomistic-level understanding of the thermal stability of exfoliated black phosphorus and defines the temperature window available for material and device processing.

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of contents Image

Keywords: Black phosphorus; sublimation; in situ electron microscopy; 2D materials; Nanoribbons; Atomistic Simulations; Density Functional Theory.


Page 2 of 24

Page 3 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Black Phosphorus (bP), a van der Waals (vdW) solid that can be exfoliated to form a few atomic layer-thick structures, is now being considered as a novel family of two-dimensional (2D) semiconducting materials1. 2D-bP exhibits some interesting properties including a tunable direct bandgap in the 0.30-2.05 eV range, depending on the number of layers2-3, and a relatively high hole mobility4 in addition to very interesting anisotropic properties2, 5-6. However, to this day, most of the predictions and properties as 2D materials remain difficult to test experimentally compared to graphene and other mainstream 2D materials7. The main difficulty has to do with the thickness dependent instability in ambient conditions8, which hinders the development of a reliable production process to generate high quality 2D-bP having few layers only. Overcoming these challenges is of paramount importance to develop a deep understanding of the basic properties of bP and assess its technological relevance. With this perspective, recent investigations have shed new lights on bP on the environmental degradation and various strategies have been proposed to enhance its stability upon exposure to ambient atmosphere9-10. As most synthesis methods and conventional device processing require relatively high temperature treatments, it is crucial as well to understand the thermal stability and degradation processes of 2D-bP layers. Atomic- and nano-scale materials often exhibit different thermal properties compared to bulk, mostly due to the increasing importance of surface effects as dimensions shrink. In this regard, vdW and 2D materials are self-passivated and hence expected to behave differently: their thermal stability should depend principally on the presence of edges and on the atomic structure of these edges. For instance, multilayer graphene exhibits upon annealing a layer-by-layer sublimation with preferential formation of zigzag edges11-14. For bulk bP, the first studies were carried out in the early 20th century. It was demonstrated that bP is the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

most stable form of phosphorus and can sustain thermal annealing up to 560°C15. Nevertheless, a significant vapor pressure was measured by Bridgman at temperatures as low as 357.1°C16, which raises questions about the accuracy of the data reported in early studies15. It was also shown that bP partially converts to violet phosphorus through evaporation and condensation starting from 445°C17. For exfoliated bP, based on in situ transmission electron microscopy (TEM), it was recently observed that the sublimation of bP occurs at 400°C18. In this letter, we elucidate the exact thermal behavior of exfoliated bP using detailed real time observations of structural and morphological changes of bP upon thermal annealing and define the temperature window for its sublimation. Moreover, we describe the energetics of different pathways of P atoms detachment and provide the detailed atomistic-level mechanism governing the sublimation of bP. Herein, the dynamics of sublimation are studied using in situ Low Energy Electron Microscopy (LEEM). The experimental observations are discussed on the basis of detailed complementary theoretical calculations and simulations using density functional theory (DFT) and kinetic Monte Carlo (KMC). Finally, we present an analytical model that describes well the general behavior and the role of defects in the observed sublimation process. LEEM analyses demonstrate that exfoliated bP is thermally stable up to 375±20°C, the temperature at which the first signs of sublimation are observed. Sublimation occurs at the edges of layers as well as at the surface. The surface sublimation takes place in the form of expanding asymmetric holes in agreement with early observations18.

Interestingly, LEEM movies

demonstrate that the thermal decomposition of bP takes place through successive sublimation of individual layers as illustrated in Figs. 1a-1h. This figure exhibits snapshots of a representative LEEM movie recorded at an increasing temperature and at time intervals of two seconds. We notice that an asymmetric hole appears in Fig. 1b and rapidly expands as shown in Fig. 1c and 4 ACS Paragon Plus Environment

Page 4 of 24

Page 5 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1d. In Fig. 1e, a second concentric hole forms on the second layer and begins to expand in a similar fashion to the first one. Then, as displayed in Fig. 1g, a third concentric hole appears on the third exposed layer. The dark lines (edges of the hole) are typical of single atomic steps suggesting that each hole is most likely one monolayer-deep. The shape varies between holes (see Figs. S1-S3 in Supplementary Information (SI)), but all of them are elongated with long axis to short axis ratios in the 1.5-2.3 range. Using the low energy election diffraction (LEED) (Fig. 1i), we determine that the long axis is parallel to the [100] crystallographic direction for all observed holes. This contrasts with the [001] direction reported earlier18. The origin of this difference is unclear. However, we carried out a set of calibration LEEM/LEED tests on reference structures indicating that diffraction and imaging modes are correctly aligned (SI) and therefore confirming that the long axis of the expanding holes is aligned along the [100] direction. Some of the observed holes display clear facets (see Fig. 1h for example) with long straight edges parallel to the [100] direction while others are rounded. All holes display sharp tips pointing to ±[100] direction with edges composed of alternating (1,1) and (1,0) steps as confirmed in the simulations discussed below. The dynamics of the hole expansion is exhibited in Figs. 2a and 2b. The figure demonstrates that the displacements along both short and long axes vary linearly with time. The slopes correspond to sublimation velocities at 412 ± 20 °C and indicate velocities of = 1.16 ± 0.10 / along the long axis and = 0.74 ± 0.10 / in the short one. Fig. 2b shows an Arrhenius plot of the sublimation velocity recorded at different temperatures indicating an activation energy of = 1.64 ± 0.10. In the following, we discuss the basic mechanisms of P detachment throughout the sublimation process by evoking the different plausible pathways displayed in Fig. 3. The first model (M1) (Fig. 3a) assumes that the 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sublimation occurs through the detachment of single P atoms18, whereas the second model (M2) involves the detachment of pairs of P atoms (P2) (Fig. 3b). In both models, only edge atoms can sublimate and the out-of-plane interactions between different layers are neglected. We will then discuss the relevance of these models based on LEEM, DFT and KMC results. The main difference between the two models is that M1 considers the sublimation of individual P atoms, whereas M2 involves the sublimation of P2 molecules. As speculated by Liu et al.18, M1 can explain to a certain extent the anisotropy of the holes. However, as demonstrated below, our results indicate that sublimation through single atom detachment is unlikely to occur. According to M1, sublimation is easier if an atom on the edge has only one bond (process #1 in Fig. 3a)) and becomes harder if the atom has two bonds. In addition, in that work18, the authors define two possible configurations for atoms with two bonds: i) Atom with two out-of-plane bonds, labeled as process #2 in Fig. 3a, where a P atom either from the top (in gray) or from the bottom (in black) is bonded to two dissimilar P atoms; ii) Atom with two in-plane bonds, process #3 in Fig. 3a, where the two interatomic bonds involve identical P atoms either from top or bottom part of the layer. In contrast, there are five different configurations for P2 detachment considered in M2. Processes #4 and #5 (Fig. 3b) involve the detachment of P2 molecules from a zigzag edge. Processes #6 and #7 correspond to the sublimation of P2 on an armchair edge and process #8 to the sublimation of a P2 molecule on a (1,1) edge. The distinction between processes #4 and #5 and between #6 and #7 comes from the presence of nearest P2 neighbor on the zigzag and armchair edges, respectively. This will be discussed in more details later on along with the DFT results. Herein, the anisotropy of the sublimation should come from different activation energies of those processes and not directly from the anisotropy of the structure, as the sharp tip and large aspect ratio cannot be simply attributed to the anisotropic atomic density. 6 ACS Paragon Plus Environment

Page 6 of 24

Page 7 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To determine the sublimation activation energies, we estimated the energetics of the different processes described above using DFT. To model the edges, we introduce bP nanoribbon (PNR) structures for which the edge properties and stability have been extensively studied19-26. Here, we computed the energy differences, ∆, between the perfect continuous edge of a model PNR and that of the same edge after one sublimation event, which is the lower limit for the actual sublimation activation energy :

=

!"

≤

(1)

where $ is the total energy of the initial PNR with % atoms, $"& is the total energy of the PNR from which we removed two (one for each edge) '& (( = 1 )* 2), + is the total energy of an isolated ' or ' in vacuum and the 1/2 factor is to count for the fact that we remove '& on both side of the PNRs. As shown in Fig. 4, PNRs with zigzag edges (ZZ-PNRs), armchair edges (AC-PNRs), and edges along the [110] direction (11-PNR) are used to compute sublimation energy at the edges. Furthermore, a (4×4) monolayer bP supercell is used to calculate the sublimation energies at the surface. Energies of processes #4a and #6a were obtained based on a simple energy conservation argument and can be calculated as: Δ = 2 -+ − + , where -+ is the total energy per atom for monolayer bP and + the total energy of a ' molecule in vacuum. This can be understood by considering an infinite asymmetric ribbon having a width of N ' rows from position ] − ∞, −1] (positions along the ribbon) and a width of N+1 ' rows from position [0, ∞[. When we remove a ' at position 0, we obtain the exact same configuration due to infinite boundaries. Therefore, the energy variation must be the cohesive energy difference between bP and ' . We then examine the energy difference ∆ as well as the structural



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 24

relaxation to determine which atomic processes are most likely to occur during the sublimation. The results are shown in Table 1 and will be discussed below in the light of KMC simulations. Table 1 Calculated DFT and fitted KMC energies for different sublimation processes. For each model, the left column represents the energy calculated using equation (1) or (2) on the modeled PNRs. The right column represents the fitted KMC activation energies. For M1, an ellipse is fitted on the hole every n2 iteration to track the hole expansion dynamics. For M2, only the extremity of the holes are followed. The fitted KMC values allow to reproduce the 1.8 long to short axis ratio and the 21.5nm/s velocity measured at 495°C within 1% error. For shape optimization, 245 000 atoms are sublimated and 10 simulations are averaged. For velocity calculations, 40 000 atoms are sublimated and 1000 simulations are averaged Model #1

Model #2

23 (DFT)

34 (KMC fit)

23 (DFT)

34 (KMC fit)

Process #1

2.32 eV

≪ 1.62 eV

---------------

-----------------

Process #2

3.43 eV

1.57 eV

---------------

-----------------

Process #3

4.02 eV

1.62 eV

---------------

-----------------

Process #4a

----------------- -----------------

1.69 eV

1.61 eV

Process #5

----------------- -----------------

1.04 eV

1.04 eV

Process #6a

----------------- -----------------

1.69 eV

1.41 eV

Process #7

----------------- -----------------

2.30 eV

1.59 eV

Process #8

----------------- -----------------

N/A

1.66 eV

The distinction between processes #4 (a and b) and #5 can be justified as follows: the relaxed structure in Fig. 4c for the 11-PNR shows that the P-P pair in the #4a configuration (green dashed circle) relaxes and forms another bond with the neighboring ' atom, which doesn’t take place for configuration #5 (red circle in Fig. 4a). We suppose that the same thing can happen with either the left or the right neighbor in configuration #4b, which would lead to similar activation energy for both processes. Possible edge reconstructions could increase this activation energy, as shown in Ref. 20, but most likely happen on continuous [100] edges, which could explain to a certain extent why some holes have long straight [100] edges (see Fig.1h). 8 ACS Paragon Plus Environment

Page 9 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The distinction between processes #6 and #7 comes from the large difference between the stability of the two edges in Fig. 4b, which leads to ∆=1.10 eV for process #6b, which is lower than ∆=2.30 eV for process #7. This difference might be overestimated considering the short periodicity of defects (length of the supercell) along the length of PNRs. Finally, based on the obtained activation energies with KMC (explained below), no clear distinction is found between processes #6a and #6b and hence for these processes, the experimental results are simulated by KMC using the same activation energy. We note in Table 1 that energies for sublimation of single P atoms (M1) with two bonds are close to 4 eV, which is much higher than the 1.64 eV activation energy extracted from experiments. However, the sublimation energies obtained for individual processes considered in M2 are in the 1-2.3 eV range, which is in reasonable agreement with experiments. Note that the energy criterion given by Eq. (1) is not respected as our calculated energy values are slightly higher than the KMC fitted values. This difference partly comes from the PBE functional. In fact, the calculated -3.49 eV/atom cohesive energy of single layer is larger than the experimental value of -3.43eV/atom27 for bulk, which includes interlayer interactions. Moreover, some other processes, not included in our model, such as P4 detachment discussed below can modify the sublimation energies. Nevertheless, results obtained from the structural relaxation also point towards P2 sublimation, especially for the ZZ-PNR and AC-PNR. In fact, there is formation of PP dimers on the edges of the ribbons. The P-P bond lengths are 1.96 Å and 2.06 Å for ZZ-PNR and AC-PNR, respectively. This huge P-P bond contraction (9-13%) indicates that ' sublimation is likely the main sublimation pathway because ' bond lengths are 1.90 Å in gas phase and 2.26 Å in bP monolayer.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 24

The results obtained for the 11-PNR are less clear, primarily because the ribbon relaxes by more than 4 eV after P2 removal (compared to ~0.3eV for AC-NR and ZZ-NR), meaning that the ∆ calculation poorly approximates the activation energy. In addition, this ribbon orientation is not mechanically stable as it has “multiple large imaginary [vibrational] modes”19. Finally, the modelled ribbon presents an artificially short periodicity (8 atoms) of defects along the edges. Though not considered in this work, the structural relaxation in the 11-NR suggests that '6 sublimation can also take place (see blue dashed circle in Fig. 4c) on the [101] edge. '6 removal under this configuration only requires 1.16 eV, but as for ' , the PNR relaxes by >1.5eV after '6 removal and the reconstruction is complex and hard to predict for an arbitrary bP edge. Note that there is hardly an experimental method that would give direct evidence of the proposed sublimation mechanism. One may think of thermal desorption spectrometry (TDS) as a technique to elucidate the nature of species emitted throughout the sublimation of bP. However, due to their strong reactivity, P atoms should form molecules (mainly P2 and P4 to a smaller extend) immediately after their detachment from a lattice. Thus, TDS will be inconclusive about the basic mechanisms of bP sublimation. Herein, in order to validate the models discussed above, we propose to study the dynamics of sublimation for each one of the models and investigate their morphological and structural implications using KMC simulations. The results of these simulations are compared to the experimental observations and DFT calculations. For M1, we first set 7 =1eV (activation energy of process #1), as it needs to be much smaller than the other energies. We then vary and 8 to fit the experimental velocities. For M2, the activation energies were optimized to fit the velocities and to reproduce the observed morphology of the hole (sharp tip with long to short axis ratio of 1.8). Several values for activation energies reproduce the experimental results, but we restricted our search to values as 10 ACS Paragon Plus Environment

Page 11 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


close as possible to the DFT calculations. We also considered that the P-P pairs with only one bond to the bP layer can easily sublimate so they were attributed a small fixed activation energy of =1 eV. As seen in Fig. 5, KMC data demonstrate, in agreement with DFT, that M2 accurately reproduces the observed shape of the expanding hole during bP sublimation, including the characteristic sharp tip (Fig. 5b). However, the three configurations considered in M1 for single atom detachments yield systematically an ellipse-shaped hole (Fig. 5a), which is clearly not consistent with experimental observations (Fig. 5b). In addition, the activation energies allowing to reproduce the shape and velocities of the hole from LEEM data (see Table 1) are also in good agreement with DFT calculations. Note that even though 9 was not calculated in DFT work, it is reasonable to say that it has the highest activation energy because sublimation of process #8 requires in total the breaking of three P-P bonds. The good agreement between KMC simulations, DFT calculations and experimental results validates to a certain extent the sublimation pathway proposed in M2. In the following, we discuss the general trend in the thermal behavior of 2D-bP and its sensitivity to key parameters, such as thickness, size, orientation and defect density. In contrast to covalent and metallic 3D materials, one would expect that the thickness should have a limited effect on the thermal stability of a vdW material. These lamellar structures are in effect selfpassivated and do not produce dangling bonds after cleavage. Nevertheless, it has been suggested that sublimation takes place at T=550°C for bulk bP and at T=400°C for exfoliated bP, presumably due to the high surface to volume ratio18. To our knowledge, no precise studies were conducted on the sublimation of bulk bP, albeit, as mentioned earlier, a significant vapor



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 24

pressure was measured by both Bridgman16 and Smits17 starting at relatively lower temperatures of 357 °C and 445°C, respectively.

This suggests that bulk bP could well sublimate at

temperatures similar to those observed for 2D-bP. From a thermodynamic point of view, the stability of a vdW solid should not be significantly affected by its thickness due to the inherently weak out-of-plane interaction between layers. This can be understood in terms of the Gibbs or Kelvin equations28. The Kelvin equation basically states that the vapor pressure increases with the inverse of the radius of curvature, as observed in many nanoscale crystals29. In the case of an infinite flat surface like 2D-bP, this contribution is nil. In terms of the Gibbs equation, there are two contributions to the Gibbs energy, given by :; = :-?=@ + :?B;CD

28

. For films in

vacuum at 0 K, the only non-zero contribution to the Gibbs energy is the cohesive energy. Considering a cohesive energy of -3.57 eV/atom for bulk and of 11.2 meV/Å2 for the surface30, the relative contribution of :?B;CD to :; is only -0.11% for 10 layers and -1.1% for one monolayer of 2D-bP. In general, the sublimation dynamics depends on several factors such as the type of defects and the crystal size. We think that hole nucleation occurs mostly on existing defects, such as vacancy-like defects. In fact, Fig. S1 shows that all holes on the first layer have approximately the same size, meaning that all hole nuclei were already present when sublimation started. This agrees with our sublimation energy calculations (∆E=5.50 eV for ' sublimation and 5.00 eV for ' sublimation), excluding the possibility of bulk sublimation. Furthermore, we see that holes usually nucleate immediately or shortly after their nuclei are exposed to vacuum thus supporting the hypothesis of heterogeneous hole nucleation. There are plausibly different types of defects that could extend through multiple layers and serve as hole nucleation sites. Fig. S2 shows two types of multilayer defects. The first one is a defect extending across the thickness of a few 12 ACS Paragon Plus Environment

Page 13 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


layers, where multiple hole nucleate at the same point for each layer. The second one is a linear defect along the surface extending through the thickness of a few layers, where many holes nucleate on the same line at each layer. To describe the general behavior during the sublimation, a simplified model was developed consisting of a continuous film of multilayer bP exposed to the vacuum with holes nucleating on the top layer at randomly distributed point defects immediately after an abrupt increase of temperature at E = 0 to trigger the sublimation. On the layers below, holes nucleate at randomly distributed point defects immediately after they are exposed to vacuum after the sublimation of top layers. For simplicity the holes considered here are of circular shape, without affecting the general sublimation behavior. The flux of P2 molecules from the surface to vacuum is the sum of the fluxes from all layers. The flux from the first layer is given by (see SI for details): F7 = G ∗ ∗ IJ = IJ (2L EIM ) "OPQ R

S

(2)

where is the hole expansion velocity, G is the total length of the hole edges on the surface, IJ is the surface atomic density, IM is the surface defect density and E is the time. For multilayer sublimation (details of the treatment are presented in SI), we find that the flux initially depends on IM E, thus it is a much more temperature sensitive process than sublimation of a stepped surface for example (due to the exp (−2 /WX) dependence). After this initial stage, which lasts for roughly 0.5 monolayer, the flux still increases asymptotically and reaches 95% of its maximum value at about 10 monolayers. When steady state is reached, the flux depends on IMY.Z . From reflection high-energy electron diffraction (RHEED) intensity simulations (Fig. S10), we can conclude that the sublimation has an initial weak layer-by-layer mode and tends 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 24

towards a multilayer sublimation mode (SI). We can estimate that defect-assisted sublimation is the dominating mechanism for large (monolayer) films down to a lateral size of *C ≈ 1.60 IM"Y.Z ≈ 358 , with *C is the radius of a 2D crystal for which the half-life is the same whether we consider sublimation at the edge or at surface defects. We used IM = 20 _ " extracted from Fig. S1b. This *C is very large, meaning that the reduced lateral size has a much bigger impact on thermal stability than the film thickness. For smaller crystals, the stability will be mostly determined by the presence of different edges, as well as different edges reconstructions, which can greatly increase the stability as suggested recently20. In summary, in situ LEEM observations in combination with KMC and DFT calculations were carried out to elucidate the sublimation process of exfoliated bP and its dynamics as a function of temperature. We found that the exfoliated bP remains stable below ~375°C. Above this temperature, a thermal decomposition takes place through the successive sublimation of individual layers. This phenomenon involves the formation and expansion of asymmetric holes having a long axis aligned along the [100] direction and sharp tips with edges composed of alternating (10) and (11) steps. The expansion of these holes is a thermally activated process with an activation energy of ~1.6 eV. Experiments and calculations demonstrate that the sublimation implies detachments of pairs of P atoms rather than the previously hypothesized single atoms. The role of defects on hole nucleation was also discussed based on an analytical model describing the general behavior of the observed sublimation. In addition to expanding our fundamental understanding of the basic structural and morphological stability of the emerging 2D-bP, this study also provides an important piece of quantitative, detailed atomistic-level information which will serve as a strong base to understand and design


Page 15 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


experiments involving high temperatures, such as epitaxial growth, which is the key step to power up 2D-bP research and make applications possible.

Methods Layers of 2D-bP were exfoliated from bulk substrates and then transferred onto Si(100) with the scotch tape method in a glovebox using polydimethylsiloxane (PDMS). Si host substrates were previously cleaned for 10 min in piranha solution followed by a 1 min HF dip (1%) to remove the native oxide. To limit the exposure to ambient air and avoid oxidation, the samples were transferred within a few minutes from the glovebox to the LEEM chamber in a paraffin sealed container. The samples were then transferred into the LEEM chamber under a constant N2 flux which was pumped immediately. All subsequent steps were carried out under UHV environment (10-10-10-9 mbar). LEEM observations were performed in bright field mode at an electron energy between 30 and 35 eV. Due to the low energy of the electrons, beam damage is thought to be minimal, as compared with TEM or SEM. Local heating by the beam is also minimal due to the low current and low electron energy. The instrument optimal spatial resolution is ~5 nm. DFT calculations were carried out with Quantum Espresso31 using the projector augmented wave method32 with the PBE exchange-correlation functional33. To model bP lattice we introduced nanoribbons (PNRs) with supercell geometries of the following dimensions: (2×1×6), 48 atoms, for zigzag PNRs (ZZ-PNR), (7×1×2) for armchair PNRs (AC-PNRs) and (7×1×2) with a monoclinic supercell for the ribbons parallel to the [110] direction (11-PNR). 10 Å or more vacuum spacing was used in all non-periodic directions of the supercell and along with a 544eV plane wave energy cutoff. All structures were relaxed until all components of all forces were smaller than 0.01eV/Å. The criterion for electronic convergence was set at 1.36_ 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 24

and a 0.0136eV Gaussian smearing width of the Fermi surface was used. Monkhorst-Pack kpoint meshes for Brillouin zone sampling34 were optimized separately for each structure with bP monolayer equivalent densities of: (32×1×9) for AC-PNRs, (9×2×32) for ZZ-PNRs, and (12×1×9) for monolayer bP. For the 11-PNR, we used a (2×1×16) k-point mesh (actual density for the supercell). For KMC simulations, attempt frequencies were chosen to be 2.44 × 1078 Hz corresponding to the average Debye frequency as extracted from reference35. The initial structures consist of a bP monolayer with either an atomic or diatomic vacancy defect to initiate sublimation. To simulate the dynamics of sublimation, 40,000 atoms are sublimated and the resulting morphology is fitted after every n2 events to track the evolution in time. 1,000 simulations are averaged and the activation energies are used as fitting parameters to reproduce the experimental observations.

Acknowledgment. The authors thank A. Favron, É. Gaufres, P. Desjardins, S. Francoeur, A. Rochefort, and M. Hersam for fruitful discussions. OM acknowledges support from NSERC-Canada (Discovery Grants), Canada Research Chair, Fondation de l’École Polytechnique de Montréal, MRIF Québec, Calcul Québec, and Computer Canada. RM acknowledges funding from NSERC, Université de Montréal, UdeM Research Chair Program. OM and RM are grateful to l’Institut de l’Énergie Trottier for supporting research on black phosphorus. Computations were made on the supercomputer Briarée from Université de Montréal, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), ministère de l’Économie, de la Science et de l’Innovation du Québec (MESI) and the Fonds de recherche du Québec - Nature et technologies (FRQ-NT). Supporting information description Supporting information includes additional LEEM images of the holes formed by bP sublimation, as well as holes nucleating on different types of defects. Two independent calibrations of the rotational angle between LEEM and LEED modes are then presented. Finally, the derivation of the analytical model describing the general sublimation behavior is presented.


Page 17 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References 1. Liu, H.; Du, Y.; Deng, Y.; Peide, D. Y. Semiconducting black phosphorus: synthesis, transport properties and electronic applications. Chem. Soc. Rev. 2015, 44, 2732-2743. 2. Liu, H.; Neal, A. T.; Zhu, Z.; Luo, Z.; Xu, X.; Tománek, D.; Ye, P. D. Phosphorene: An Unexplored 2D Semiconductor with a High Hole Mobility. ACS Nano 2014, 8, 4033-4041. 3. Qiao, J.; Kong, X.; Hu, Z.-X.; Yang, F.; Ji, W. High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus. Nat. Comm. 2014, 5, 4475. 4. Li, L.; Yu, Y.; Ye, G. J.; Ge, Q.; Ou, X.; Wu, H.; Feng, D.; Chen, X. H.; Zhang, Y. Black phosphorus field-effect transistors. Nat. Nanotechnol. 2014, 9, 372-377. 5. Qin, G.; Yan, Q.-B.; Qin, Z.; Yue, S.-Y.; Hu, M.; Su, G. Anisotropic intrinsic lattice thermal conductivity of phosphorene from first principles. Phys. Chem. Chem. Phys. 2015, 17, 4854-4858. 6. Wang, X.; Jones, A. M.; Seyler, K. L.; Tran, V.; Jia, Y.; Zhao, H.; Wang, H.; Yang, L.; Xu, X.; Xia, F. Highly anisotropic and robust excitons in monolayer black phosphorus. Nat. Nanotechnol. 2015, 10, 517-521. 7. Butler, S. Z.; Hollen, S. M.; Cao, L.; Cui, Y.; Gupta, J. A.; Gutiérrez, H. R.; Heinz, T. F.; Hong, S. S.; Huang, J.; Ismach, A. F. et al. Challenges, and Opportunities in Two-Dimensional Materials Beyond Graphene. ACS Nano 2013, 7, 2898-2926. 8. Favron, A.; Gaufrès, E.; Fossard, F.; Phaneuf-L’Heureux, A.-L.; Tang, N. Y.; Lévesque, P. L.; Loiseau, A.; Leonelli, R.; Francoeur, S.; Martel, R., Photooxidation and quantum confinement effects in exfoliated black phosphorus. Nat. Mater. 2015, 14, 826-832. 9. Kim, J.-S.; Liu, Y.; Zhu, W.; Kim, S.; Wu, D.; Tao, L.; Dodabalapur, A.; Lai, K.; Akinwande, D. Toward air-stable multilayer phosphorene thin-films and transistors. Sci. Rep. 2015, 5, 8989. 10. Wood, J. D.; Wells, S. A.; Jariwala, D.; Chen, K.-S.; Cho, E.; Sangwan, V. K.; Liu, X.; Lauhon, L. J.; Marks, T. J.; Hersam, M. C. Effective passivation of exfoliated black phosphorus transistors against ambient degradation. Nano Lett. 2014, 14, 6964-6970. 11. Barreiro, A.; Börrnert, F.; Rümmeli, M. H.; Büchner, B.; Vandersypen, L. M. K. Graphene at High Bias: Cracking, Layer by Layer Sublimation, and Fusing. Nano Lett. 2012, 12, 1873-1878. 12. Huang, J. Y.; Qi, L.; Li, J. In situ imaging of layer-by-layer sublimation of suspended graphene. Nano Res. 2010, 3, 43-50. 13. Huang, J. Y.; Ding, F.; Yakobson, B. I.; Lu, P.; Qi, L.; Li, J. In situ observation of graphene sublimation and multi-layer edge reconstructions. Proc.Nat. Acad. Sci. 2009, 106, 10103-10108. 14. Jia, X.; Hofmann, M.; Meunier, V.; Sumpter, B. G.; Campos-Delgado, J.; Romo-Herrera, J. M.; Son, H.; Hsieh, Y.-P.; Reina, A.; Kong, J. Controlled formation of sharp zigzag and armchair edges in graphitic nanoribbons. Science 2009, 323, 1701-1705. 15. Jacobs, R. B., Phosphorus at High Temperatures and Pressures. J. Chem. Phys. 1937, 5, 945-953. 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

16. Bridgman, P. W., Two New Phenomena at Very High Pressure. Phys. Rev. 1934, 45, 844-845. 17. Beck, R. P.; Meyer, G.; Smits, A. On Black Phosphorus. I. Koninklijke Nederlandse Akademie van Wetenschappen Proceedings Series B Physical Sciences 1915, 18, 992-1007. 18. Liu, X.; Wood, J. D.; Chen, K.-S.; Cho, E.; Hersam, M. C. In Situ Thermal Decomposition of Exfoliated Two-Dimensional Black Phosphorus. J. Phys. Chem. Lett. 2015, 6, 773-778. 19. Ramasubramaniam, A.; Muniz, A. R. Ab initio studies of thermodynamic and electronic properties of phosphorene nanoribbons. Phys. Rev. B 2014, 90, 085424. 20. Liang, L.; Wang, J.; Lin, W.; Sumpter, B. G.; Meunier, V.; Pan, M. Electronic Bandgap and Edge Reconstruction in Phosphorene Materials. Nano Lett. 2014, 14, 6400-6406. 21. Li, W.; Zhang, G.; Zhang, Y.-W. Electronic Properties of Edge-Hydrogenated Phosphorene Nanoribbons: A First-Principles Study. J. Phys. Chem. C 2014, 118, 22368-22372. 22. Zhu, Z.; Li, C.; Yu, W.; Chang, D.; Sun, Q.; Jia, Y. Magnetism of zigzag edge phosphorene nanoribbons. Appl. Phys. Lett. 2014, 105, 113105. 23. Carvalho, A.; Rodin, A. S.; Neto, A. H. C., Phosphorene nanoribbons. EPL 2014, 108, 47005. 24. Xu, L.-C.; Song, X.-J.; Yang, Z.; Cao, L.; Liu, R.-P.; Li, X.-Y. Phosphorene nanoribbons: Passivation effect on bandgap and effective mass. Appl. Surf. Sci. 2015, 324, 640-644. 25. Guo, H.; Lu, N.; Dai, J.; Wu, X.; Zeng, X. C. Phosphorene Nanoribbons, Phosphorus Nanotubes, and van der Waals Multilayers. J. Phys. Chem. C 2014, 118, 14051-14059. 26. Du, Y.; Liu, H.; Xu, B.; Sheng, L.; Yin, J.; Duan, C.-G.; Wan, X. Unexpected Magnetic Semiconductor Behavior in Zigzag Phosphorene Nanoribbons Driven by Half-Filled One Dimensional Band. Sci. Rep. 2015, 5, 8921. 27. Kittel, C. Introduction to solid state physics 8th edn, Ch. 3, 49–88 (Wiley, Hoboken, NJ, USA, 2004). 28. Kaptay, G. The Gibbs equation versus the Kelvin and the Gibbs-Thomson equations to describe nucleation and equilibrium of nano-materials. J. Nanosci. Nanotechnol. 2012, 12, 26252633. 29. Buffat, P.; Borel, J. P. Size effect on the melting temperature of gold particles. Phys. Rev. A 1976, 13, 2287. 30. Shulenburger, L.; Baczewski, A. D.; Zhu, Z.; Guan, J.; Tománek, D., The Nature of the Interlayer Interaction in Bulk and Few-Layer Phosphorus. Nano Lett. 2015, 15, 8170–8175. 31. Paolo, G.; Stefano, B.; Nicola, B.; Matteo, C.; Roberto, C.; Carlo, C.; Davide, C.; Guido, L. C.; Matteo, C.; Ismaila, D. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condensed Matter 2009, 21 (39), 395502. 32.

Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953-17979.


Page 18 of 24

Page 19 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


33. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865-3868. 34. Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188-5192. 35. Zhu, Z.; Tománek, D. Semiconducting Layered Blue Phosphorus: A Computational Study. Phys. Rev. Lett. 2014, 112, 176802.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. Bright-field LEEM snapshots of hole expansion during sublimation of exfoliated bP. 2 seconds between each image from (a) to (h) recorded respectively at the following temperatures: 486°C, 488°C, 490°C, 491°C, 493°C, 495°C, 497°C and 499°C. Dashed blue line highlighting the hole in (d) with black double arrow showing its short axis and red double arrow showing its long axis. Scale bar is 200 nm. (i) The corresponding LEED pattern before sublimation. The bP reciprocal lattice vectors are also indicated. (j) Illustrations of bP direct lattice unit cell vectors and of the early stage of sublimaiton-induced hole on top layer (in black) with bottom layer (in blue) still intact. 20 ACS Paragon Plus Environment

Page 20 of 24

Page 21 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. (a) Average displacements along the long and short axis of holes recorded at T=412± 20°C; (b) Arhenius plot of the short axis velocity, c



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Illustrations of the two sublimation models considered in this work. (a) Model of single P atoms sublimation. (b) Model of P2 sublimation. The different atomic processes considered for each model are color-coded and labeled from 1 to 8. The faded atoms are the ones that sublimate and the colored brighter ones are the sublimated atoms.


Page 22 of 24

Page 23 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Structures of the PNRs models used for the DFT calculations. The left side of the ribbon is the initial relaxed configuration. The right side is the relaxed ribbon after sublimation of a P2 molecule. The red circles are the P2 to be sublimated: (a) ZZ-PNR (b) AC-PNR (c) 11PNR. 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. KMC simulated holes after 106 atoms are sublimated for: (a) M1 (b) M2. (c) LEEM image (460 × 720 nm2, see Fig.1d for a complete time evolution). Note that the simulated aspect ratios are not exactly the same as in (c), but they are within the range of experimentally recorded values.


Page 24 of 24

Dynamics and Mechanisms of Exfoliated Black Phosphorus

Recommend Documents