Phonon Anharmonicity in Few-Layer Black Phosphorus

Subscriber access provided by Duke University Libraries

Article

Phonon Anharmonicity in Few-Layer Black Phosphorus Damien Tristant, Andrew Cupo, Xi Ling, and Vincent Meunier ACS Nano, Just Accepted Manuscript • DOI: 10.1021/acsnano.9b04257 • Publication Date (Web): 22 Aug 2019 Downloaded from pubs.acs.org on August 23, 2019

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 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 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.

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 40 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

ACS Nano

Phonon Anharmonicity in Few-Layer Black Phosphorus Damien Tristant,† Andrew Cupo,† Xi Ling,‡ and Vincent Meunier∗,† †Department of Physics, Applied Physics, and Astronomy, Rensselaer Polytechnic Institute, Troy, New York 12180, United States of America ‡Department of Chemistry, Division of Materials Science and Engineering, and The Photonics Center, Boston University, Boston, Massachusetts 02215, United States of America E-mail: [email protected] Abstract We report a temperature-dependent Raman spectroscopy study of few-layer black phosphorus (BP) with varied incident polarization and sample thickness. The Ramanactive modes A1g , B2g , and A2g exhibit a frequency downshift, while their linewidth tends to increase with increasing temperature. To understand the details of these phenomena, we perform first-principles density functional theory calculations on freestanding monolayer BP. The effect of thermal expansion is included by constraining the temperature dependent lattice constant. The study of the temperature-induced shift of the phonon frequencies is carried out using ab initio molecular dynamics simulations. The normal mode frequencies are calculated by identifying the peak positions from the magnitude of the Fourier transform of the total velocity autocorrelation. Anharmonicity induces a frequency shift for each individual mode and the three- and four-phonon process coefficients are extracted. These results are compared with those obtained from many-body perturbation theory, giving access to phonon lifetimes and

1

ACS Paragon Plus Environment

ACS Nano 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

lattice thermal conductivity. We establish that the frequency downshift is primarily due to phonon-phonon scattering while thermal expansion only contributes indirectly by renormalizing the phonon-phonon scattering. Overall, the theoretical results are in excellent agreement with experiment, thus showing that controlling phonon scattering in BP could result in better thermoelectric devices or transistors that dissipate heat more effectively when confined to the nanoscale.

Keywords: Raman spectroscopy, first-principles, phonon-phonon coupling, anharmonicity, frequency shift, phonon lifetime, black phosphorus Heat transfer in non-metals is dominated by phonons, the quanta of atomic vibrations. Phonons, also called normal modes of vibration, are characterized by their frequency ω(q) where q is the quasi-momentum in the material’s first Brillouin zone (BZ). At zero temperature and neglecting the quantum nature of the ion cores, the vibrational properties of nanostructures are usually studied within the harmonic or quasi-harmonic approximations. In this case, harmonic phonons are non-interacting, they have an infinite lifetime, and, as a result, the corresponding thermal conductivity is infinite. At finite temperature, this representation is no longer an accurate description of the physical properties of a crystal lattice as an infinite thermal conductivity is inconceivable. Therefore, a more realistic approach calls for a framework where phonon anharmonicity is considered by including phonon-phonon scattering. Experimentally, Raman spectroscopy gives access to the temperature dependence of the phonon frequency shifts and linewidths. Raman is a non-destructive and powerful method to investigate and understand anharmonic contributions and it has been successfully applied to many classes of materials, notably two-dimensional (2D) van der Waals (vdW) layered and three-dimensional (3D) materials. 1–8 Among those, of recent interest is few-layer black phosphorus (BP), which is an AB-stacked layered material with in-plane anisotropy (see Fig. 1(a)). BP has generated a growing interest for applications in nanoelectronic devices due in part 2


Page 2 of 40

Page 3 of 40 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

ACS Nano

to its tunable direct electronic band gap, carrier mobility up to 55,000 cm2 /(Vs) at cryogenic temperatures, and anisotropic thermal conductivity. 9–11 One major challenge faced by the semiconductor device industry is the fabrication of nanoscale semiconductor components with high thermal conductivities that can dissipate heat effectively. Recently, F. Nasri et al., compared the thermal properties of a conventional silicon MOSFET transistor with that of another one made from BP. 12 It was demonstrated that the BP-based device featured a net temperature reduction of ∼5% compared to a conventional one. This result is due in large part to the improvement of the thermal conductivity. For example, at room temperature a 300 nm-thick BP sample was found to exhibit a thermal conductivity up to 86 W/(mK) along the zigzag (ZZ) direction. 13 Unfortunately, this value decreases with increasing temperature. 14 To understand this phenomenon from a fundamental standpoint requires a study dedicated to the behavior of phonons as a function of temperature. Like many other 2D materials, the physical properties of BP, such as current modulation, electronic mobility, photoresponsivity, and response time, depend on the number of layers. 15 For this reason it is important to accurately characterize the thickness of BP. One possibility is to examine the ultra-low frequency, Raman active, interlayer modes as it was shown that the breathing mode frequencies of BP provides an accurate measure of the number of layers. 16,17 Theoretically, X. Ling et al., 18 quantified this variation by showing that increasing the number of layers from 2 to 8 decreases the frequency associated with the main breathing mode by ∼38 cm−1 . Besides the behavior of these low-wavenumber modes, Davydov-splitting of the A2g mode around 470 cm−1 was reported in the spectra of few-layer BP (2-4 layers). 19,20 Conversely, no substantial frequency shift is observed at high frequency when the thickness of BP increases. Those high-frequency Raman-active intralayer modes are represented in Fig. 1(b). Note that an important characteristics of Raman spectroscopy is the choice of polarization for the laser beam. For symmetry reasons, some modes, such as B2g are not visible for a laser polarized in the direction parallel or perpendicular to the atomic displacements. 21–25

3


ACS Nano 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 first experimental research conducted on the temperature dependence of Raman spectra for few-layer BP was limited to the linear regime, with an extracted first-order temperature coefficient down to -0.028 cm−1 /K. 26 The temperature range was subsequently extended to where the frequency shifts could be fit to a model based on the optical phonon decays. 8,22,27 However, a lack of theoretical knowledge limits the understanding of the vibrational properties of this 2D material. This paper reports a Raman spectroscopy study of few-layer BP with varied incident polarization for two different thicknesses, 5 and 20 nm. We show that the Raman active modes A1g , B2g , and A2g exhibit a frequency downshift, while the linewidth tends to increase with increasing temperature. These experimental findings are rationalized using ab initio molecular dynamics (AIMD) simulations based on density functional theory (DFT) calculations on freestanding monolayer BP. The effect of thermal expansion is included by constraining the temperature dependent lattice constant. By fitting to an established model for the optical phonon decay, we determine the cubic (three-phonon processes) and quartic (four-phonon processes) anharmonic coefficients, which both influence the frequency shift and linewidth of the high-frequency intralayer modes. We find that the thermal expansion is not the primary reason for the observed redshift. These results are compared to data computed using many-body perturbation theory, giving access to phonon lifetimes and lattice thermal conductivity.

Harmonic, Anharmonic, and Quasi-Harmonic Effects First-principles simulations can be used to study the atomic vibrational properties numerically, in a crystal at its equilibrium state or subject to different perturbations from the environment (e.g., variations of temperature, T , or pressure, P ). Based on the BornOppenheimer approximation, the equations which determine the electronic state treat the ionic positions as parameters and are therefore decoupled from those related to the ionic

4


Page 4 of 40

Page 5 of 40

(a)

(b)

A1g

B2g

A 2g

Zigzag

(c)

! P45°-AC/ZZ

! PZZ

! PAC Armchair

(d)

-15

-10

-5

Y (um)

Y (µm)

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

ACS Nano

0

A

5

10

2 µm

B

15 -10

0 X (um) (µm)

10

Figure 1: (a) Perspective view of AB-stacked bulk black phosphorus (BP). Phosphorus atoms in different planes are represented by different colors. (b) Schematic representations of the high-frequency Raman-active intralayer modes, A1g (out-of-plane), B2g (along zigzag direction), and A2g (along armchair direction). The cross and circle indicate vibrations going into the plane of the page and coming out of it, respectively. (c) Top view of single-layer BP. Polarization of the laser beam along: zigzag (ZZ) direction, armchair (AC) direction, and 45◦ between AC and ZZ directions. (d) Typical optical image of few-layer BP. The thickness of the areas A and B are about 5 and 20 nm, respectively.

5


ACS Nano 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 6 of 40

motion. It follows that the total potential energy, U , only results from the electronic ground state for a given geometrical configuration and can be expanded in a Taylor series around its minimum value as,

U = Ueq +

∞ X X (n) 1 X ... Φα1 l1 κ1 ...αn ln κn uα1 l1 κ1 ...uαn ln κn . n! n=2 α l κ α l κ

(1)

n n n

1 1 1

Each atom (labeled by the letter κ) in the unit cell l is allowed to vibrate with a displacement uαlκ along the Cartesian directions α = {x, y, z}. The first-order derivative is zero since a stable ionic configuration corresponds to a relative minimum of the potential. The higher-order derivatives are included in the last sum on the right-hand side and contains the harmonic (n = 2) and anharmonic (n > 2) contributions. The elements of the interatomic force constant matrix, evaluated in the equilibrium configuration, are given by: (n) Φα1 l1 κ1 ...αn ln κn

∂n U . = ∂uα1 l1 κ1 ...∂uαn ln κn 0

(2)

Different approximations are obtained depending on where this sum is truncated.

Harmonicity Historically, the term harmonicity appeared in the acoustic domain and was later extended to different fields of physics. In this case, the motion of a vibrating string is described in terms of sine and cosine functions, which are referred to as harmonics. 28 This regime corresponds to the situation where the largest n value in Eq. 1 is equal to 2. This is equivalent to a situation where the atomic displacements are linearly proportional to the force. The corresponding quadratic form of the potential energy function is thus:

Uhar =

1 X X (2) Φ uα l κ uα l κ . 2 α l κ α l κ α1 l1 κ1 α2 l2 κ2 1 1 1 2 2 2 1 1 1

2 2 2

6


(3)

Page 7 of 40 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

ACS Nano

It is straightforward to convince oneself of the general validity of this expression for small displacements around the bottom of the potential, as shown in Fig. S1 in the ESI. In a crystal, even if the total energy approximates the interaction between two atoms as a collection of spring-like connections, the harmonic model provides the main part of the vibrational entropy at low temperature. Based on the finite-displacement and supercell methods, 29,30 the harmonic frequency, denoted ω0,j (q), and the j normal modes of any phonons at point q in the first BZ can be obtained and represented in a dispersion spectrum. However, this approximation has a number of limitations. For example, the harmonic approximation does not allow for thermal expansion as the expectation value of the lattice constant remains unchanged by symmetry. Further, it does not describe the vibrational properties of a crystal subjected to external changes and the thermal conductivity is found to be infinite. The main reason for this is that phonons are non-interacting and have an infinite lifetime. Therefore, this theoretical model must be extended to take into account phonon-phonon scattering.

Anharmonicity One could think that the cubic term (n = 3) is the most important contribution to anharmonicity, where a phonon decays to form two other phonons or where two phonons combine into a single phonon after a finite time. The anharmonic potential energy for three-phonon processes is expressed as follows: (3)

Uanh =

1 X X X (3) Φ uα l κ uα l κ uα l κ . 6 α l κ α l κ α l κ α1 l1 κ1 α2 l2 κ2 α3 l3 κ3 1 1 1 2 2 2 3 3 3 1 1 1

2 2 2

(4)

3 3 3

At this stage, we only expanded Eq. 1 to third-order, but no theoretical guiding principle indicates the number of terms that ought to be included from the infinite sum. In fact, recent studies only including three-phonon coupling have shown a significant overestimation of thermal conductivity as compared to experimental data. Conversely, including the quartic

7


ACS Nano 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 40

term (n = 4) was found to reduce this difference. 31 Thus, four-phonon processes are likely to play a role as significant as the three-phonon processes and must be added to the Eq. 1 as given by the following expression: (4)

Uanh =

1 X X X X (4) Φ uα l κ uα l κ uα l κ uα l κ . 24 α l κ α l κ α l κ α l κ α1 l1 κ1 α2 l2 κ2 α3 l3 κ3 α4 l4 κ4 1 1 1 2 2 2 3 3 3 4 4 4 1 1 1

2 2 2

3 3 3

(5)

4 4 4

From perturbation theory, a fourth order process includes processes such as: a phonon decay into three other phonons (or vice versa), two phonons decay into two new phonons, and even a process involving two three-phonon decays through an intermediate state. We now introduce the physical quantity that results from these anharmonic contributions: the phonon self-energy (SE), ∆(ωj (q), T ) + iΓ(ωj (q), T ), which accounts for the renormalization of the harmonic non-interacting phonon energies due to their anharmonicity-induced mutual interaction. In their seminal papers, Maradudin et al. 32 and later Wallis et al. 33 introduced the phonon SE considering the third- and the fourth-order anharmonic terms. The real part of SE, ∆(ωj (q), T ), is associated with the frequency shift and is responsible for the temperature dependence of the phonon frequencies. The damping is described by the imaginary part Γ(ωj (q), T ) and is partially responsible for the broadening of Raman spectra. The temperature dependence of the frequency shift and the linewidth can be calculated by applying many-body perturbation theory to the exactly solvable harmonic system, with the third- and fourth-order terms from Eq. 1 treated as perturbations. Here, Γ(ωj (q), T ) is computed up to third-order. When the Raman line shape is fitted to Lorentzian functions, 2Γ(ωj (0), T ) corresponds to the full-width-at-half-maximum (FWHM) and its reciprocal is the total lifetime, denoted τj (T ). 32 Besides the phonon-phonon scattering, the electron-phonon coupling is also one of the fundamental interactions between elementary excitations in solids. In particular in metals, low-energy electronic excitations are strongly affected by the coupling to lattice vibrations. In the case of semiconductors, the electron-phonon coupling relies on the deformation po-

8


Page 9 of 40 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

ACS Nano

tential form in the rigid ion model, which can predict the temperature dependence of carrier mobility. 34 In this work we neglect the electron-phonon scattering, assuming that the most important contribution to anharmonicity originates from the phonon-phonon coupling.

Quasiharmonicity So far, we have focused on the changes in vibrational properties generated only by the interactions between phonons while neglecting variations in temperature and pressure. This contribution – thermal expansion – is not directly present in the potential energy (Eq. 1) but can be included via the calculation of the volume V (T, P ), within the quasiharmonic approximation (QHA). 35,36 This model is often misrepresented as being part of the anharmonic effects: this is not strictly correct as the quasi-harmonic approximation only includes two-phonon processes. However several investigations show that it can be used as a first approximation to study the behavior of vibrational modes at low temperature. 36 The temperature and pressure dependent equilibrium volume, V0 (T, P ), is determined by minimizing the Gibbs free energy from the following equation:

G(T, P ) = min{F (T, V ) + P V }, V

(6)

with respect to volume. The Helmholtz free energy, obtained at constant pressure, is given by: 35,36 X X ~ωj (q, V ) ~ωj (q, V ) + kB T ln 1 − exp − . F (T, V ) = E0 (V ) + 2 k BT q∈BZ,j q∈BZ,j

(7)

The first term is the ground state energy (enthalpy) of the crystal. The second and last terms are the contributions from the zero point and vibrational (entropy) energies, respectively. The sums run over all phonon wave vectors in the first BZ and all modes. ~ and kB are the reduced Planck and the Boltzmann constants, respectively. Note that in this work we did

9


ACS Nano 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

not consider the effect of pressure. Based on the QHA, we have previously determined the temperature dependent lattice parameters (a(T ) and b(T )) of freestanding single-layer BP, including vdW interactions. 37 From this knowledge, we calculated the phonon dispersion of this structure at 0 K (including the zero point energy, a(0) = 4.508 ˚ A and b(0) = 3.312 ˚ A) and at 300 K (a(300) = 4.512 ˚ A and b(300) = 3.315 ˚ A), see Fig. 2. At room temperature, taking into account the effect of thermal expansion does not significantly affect the vibrational properties of freestanding monolayer BP. Only the high-frequency modes are slightly downshifted. In fact, Maradudin et al., 32 claim that it is not necessary to introduce the effect of thermal expansion into the phonon SE since it already accounts for the higher order anharmonic terms. However, we will see later that the thermal expansion indirectly modifies the anharmonic contributions to the temperature dependent frequency shift. Here we extended the definition of the phonon SE by including the effects of thermal expansion.

Results and Discussion In this paper, we focus on the specific case of BP and illustrate the influence of temperature on the three effects outlined in the previous section.

Experiment: Anharmonicity via Raman Spectroscopy Raman spectroscopy is a powerful technique to measure ion dynamics, giving access to information related to the structural and electronic properties of materials. 39,40 Combined with a cryostat, it can be used to characterize the thermal properties of a material, such as thermal expansion and thermal conductivity, induced by a change in temperature. Here, we report Raman line positions as well as FWHM of few-layer BP with varied incident polarization along three different directions: zigzag (ZZ), armchair (AC), and at a 45◦ angle between the AC and ZZ directions (Fig. 1(c)). Furthermore we consider two different

10


Page 10 of 40

Page 11 of 40

500 500

A22g A g

0K 0K 300 K 300 K

400 400 -1-1

Phonon PhononFrequency Frequency(cm (cm ))

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

ACS Nano

B2g B 2g 300 300

3-phonon decay

A11g A g B13g

200 200

100 100

LA

00Γ

Γ

X X

SS

Y Y

Γ 0.0 0.0 Γ

DOS (a.u.) (a.u.) DOS

0.3 0.3

Figure 2: Calculated phonon dispersion (left panel) and total density of states (right panel) of single-layer black phosphorus (BP). The solid and dashed lines correspond to 0 K (including zero point energy) and 300 K (including zero point and vibrational energies), respectively. Based on the quasiharmonic approximation method, the distances between high-symmetry points (ΓX, XS, SY, and YΓ) have been renormalized. The green arrows indicate the asymmetric decay channels for the zone-center optical B2g into a longitudinal acoustic (LA) mode and a optical B13g mode. The representation of the LA and B13g modes can be found in Hackney et al. 38 thicknesses: 5 and 20 nm. The positions of the laser spot are shown in Fig. 1(d). The experimental data are directly extracted from Raman spectra and all parameters are obtained by fitting the peaks using a Lorentzian/Gaussian mixed function (see Figs. S2 and S3 in the ESI for more details). In agreement with literature results, 8,22,26,27 Fig. 3 confirms that an increasing temperature leads to a nonlinear redshift of the Raman active A1g , B2g and A2g phonon modes. Moreover, we note that this behavior does not depend on laser polarization. The measured temperature dependence of each phonon mode frequency at constant pres-

11


ACS Nano

(b)

-1

Raman Shift (cm )

-1

469 468 467 466 465 464 463 462 461 460 459 458 440

439

438

438

-1

Raman Shift (cm )

439

437 436 435 434

B2g

435 434 433 432 363

362 362

362 -1 -1

360 360 359 359 358 358

A

1 g

B2g

361 360 359 358

A1g

357

357 357 356 356

A 2g

436

433

361 361

Thickness: 20 nm

437

432 363 363

Raman Shift (cm )

-1

Thickness: 5 nm 469 468 467 466 465 464 463 60° ZZ 60° 105° 2 462 45° - AC/ZZ 150° 105° Ag AC 150° 461 460 459 458 Raman Shift (cm )

469 468 467 466 465 464 463 462 461 460 459 458 440

-1

Raman Shift (cm )

-1

Raman Shift (cm )

(a)

-1 Raman Shift (cm )

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 40

100

200 300 Temperature (K)

400

356

500

100


400

500

Figure 3: Temperature/polarization dependence of A1g , B2g and A2g Raman mode positions of black phosphorus for two thicknesses: (a) 5 nm and (b) 20 nm. Different incident polarization directions are used: zigzag (ZZ) direction, armchair (AC) direction, and 45◦ between AC and ZZ directions. Experimental values are fitted with Eq. 8 taking into account anharmonic phonon-phonon coupling and thermal expansion. sure can be decomposed as

ωj (q = 0, T ) = ω0,j (0) + ∆(ωj (0), T )anh + ∆(ωj (0), T )V ,

(8)

where the last two terms on the right-hand side describe the frequency shift caused by the anharmonic phonon-phonon coupling and the volumetric thermal expansion in the case of 3D materials, respectively. The Balkanski model can be used to fit the anharmonic part of

12


Page 13 of 40

(a)

Thickness: 5 nm 469 468 A 2g 467 466 465 464 463 462 461 460 459 458

8

5 4 3 2 7

7 -1

ZZ 60° 45° - AC/ZZ 105° AC 150°

6

A 2g

5 4 3 2 7

B2g

B2g

-1

FWHM (cm )

6

-1

FWHM (cm )

5

5

4

4

3 8 363

8 363

3627

362 7

3

-1

361 6 360 5 359 4 358

A1g

3573

2 356

100

A1g

Raman Shift (cm -1 ) FWHM (cm )

-1

Thickness: 20 nm

8

FWHM (cm )

6

6

(b)

-1

Raman Shift (cm )

-1

FWHM (cm )

7

Raman Shift (cm ) -1 FWHM (cm )

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

ACS Nano

361 6 360 5 359 4 358

3 357 200 300 Temperature (K)

400

2 356

500

100


400

500

Figure 4: Temperature/polarization dependence of A1g , B2g and A2g Raman mode full-widthat-half-maxima (FWHM) of few-layer black phosphorus for two thicknesses: (a) 5 nm and (b) 20 nm. Varied incident polarizations are used: zigzag (ZZ) direction, armchair (AC) direction, and 45◦ between AC and ZZ directions. Experimental values are fitted with Eq. 13 from the Balkanski model. 41 the frequency shift. It includes three- and four-phonon processes: 41 ∆(ωj (0), T )anh = Aj

2 1+ x e −1

+ Bj 1 +

3 3 + ey − 1 (ey − 1)2

.

(9)

The two coefficients, Aj and Bj , represent the cubic and quartic anharmonic constants. This expression depends on the Bose-Einstein phonon population where x = ~ω0,j (0)/2kB T and y = ~ω0,j (0)/3kB T .

We note that energy and momentum are both conserved in

13


ACS Nano 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 40

these processes. At temperatures approaching 0 K, we define the anharmonic frequency 0 (0) = ω0,j (0) + Aj + Bj . The volumetric thermal expansion contribution to the as ω0,j

frequency shift is expressed as follows: 42

Z

∆(ωj (0), T )V = ω0,j (0) exp −γj (0)

T

dT αV (T ) − 1 , 0

0

(10)

0

where the volumetric thermal expansion and microscopic Gr¨ uneisen parameters are calculated using the equations: 1 dV (T ) V dT

(11)

V0 dωj (q) , ω0,j (q) dV

(12)

αV (T ) = and γj (q) = −

respectively. V0 is the volume of the bulk material at zero temperature. To determine αV (T ), we employed the lattice parameters reported in the literature for BP, as measured by X-ray diffraction in the 85–320 K temperature range. 43 To obtain the lattice constants, as well as the volume as a function of temperature, the extracted experimental data are fitted using a third degree polynomial function. Then αV (T ) is calculated using Eq. 11. We estimated the microscopic Gr¨ uneisen parameters based on the experimental work of G. Xiao. et al., 44 who found γj (0) coefficients equal to 0.51, 0.20, and 0.13 at the Γ-point for the A1g , B2g and A2g modes of few-layer BP at high pressure, respectively. Gr¨ uneisen parameters depend weakly on temperature and can be considered constant when analyzing quasiharmonic effects. We use Eq. 8 (which includes both phonon-phonon scattering and thermal expansion effects) to fit each experimental data set presented in Fig. 3(a-b). The calculated values of the anharmonic constants (Aj and Bj ) and the harmonic frequency (ω0,j (0) at zero temperature) are directly extracted from the fit (R2 → 1). These results are summarized in Table 1 where we also provide a comparison with published values. Most of the time, the ratio between the two anharmonic constants, Aj /Bj , is large. This indicates a higher probability for a BP phonon to decay into two phonons instead of three. The most probable decay channels are 14


Page 15 of 40 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

ACS Nano

those for which the created vibrations have a large density of states (DOS). One of the threephonon process mechanisms is known as the Ridley channel. 45 The asymmetric decay of an optical phonon occurs into a transverse optical (TO) phonon and a longitudinal acoustic (LA) phonon subject to energy (ω) conservation. By inspection of the phonon frequencies at highsymmetry points and the DOS in Fig. 2, we find that one of the three-phonon processes that obeys energy conservation is B2g (Γ, 417.6 cm−1 ) → B13g (Γ, 220.2 cm−1 ) + LA(Y, 194.5 cm−1 ). This decay channel is represented on the vibrational band structure (see Fig. 2). From D. J. Late’s raw data, 26 we extracted these anharmonic parameters for a 4 nm-thick BP sample. We find that there is a higher probability of having four-phonon processes, especially for the A1g and A2g modes. Our work confirms this result for B2g (AC-polarization, 20 nm) and A2g (ZZ-polarization, 5 nm). Thus, it seems that there is no universal rule for predicting the strength of the Aj /Bj ratio for different choices of thickness, polarization, or normal mode. In addition, BP’s measured properties are known to be extremely sensitive to the surrounding environmen as well as oxidation. 46,47 For the sake of comparison, we also analyzed the data from Fig. 3(a-b), assuming that the frequency shift is purely due to anharmonic effects, i.e., using Eq. 9. The results are summarized in Table S1 in the ESI. We see no large change, which implies that the frequency shift is mainly due to phonon-phonon coupling. Turning to the effect of thickness on frequency shifts, for the AC polarization at the temperature of liquid nitrogen (∼78 K) the frequency of the A1g mode is downshifted by ∼3 cm−1 when the thickness is increased from 5 to 20 nm, while a difference of ∼4 cm−1 is obtained at high temperature (500 K). This implies that a change in the number of layers significantly influences the vibrational frequency of the modes, especially at high temperature. This feature has already been observed for a low-frequency mode 16 as the interlayer breathing mode A3g shows a large redshift with increasing thickness. Thus our results suggest that the A1g mode could also be used to identify the number of layers in BP at high temperature.

15


ACS Nano 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 40

We also investigated the polarization dependence on the Raman spectra. By rotating the sample, we effectively change the angle between the laser polarization and the crystal orientation of BP. Raman tensor theory makes it possible to predict the Raman cross section, k

S, which gives the scattering intensity. The one associated with the B2g mode, i.e., SB2g = (2|f | sin θ cos θ)2 , establishes that no intensity is expected when the laser polarization is along the AC (θ = 0◦ ) or ZZ (θ = 90◦ ) directions. 25 For this reason, in our case, the B2g mode is only visible when the polarization is along the direction aligned at a 45◦ angle between AC and ZZ, regardless of sample thickness, as shown in Figs. 3 and 4. We will now analyze the effect of temperature on the FWHM of each mode. The Raman broadening is related to the decay of optical phonons and can be described using the Balkanski model, expanded up to four-phonon processes: 41

Γ(ωj (0), T ) = Γ0 (ω0,j (0)) + Cj

2 1+ x e −1

+ Dj 1 +

3 3 + y y e − 1 (e − 1)2

,

(13)

where Γ0 (ω0,j (0)) is the temperature-independent harmonic linewidth and Cj and Dj represent the anharmonic parameters for line broadening from the three- and four-phonon processes, respectively. We also define the anharmonic linewidth at zero temperature, as Γ0 (ωj (0), T ) = Γ0 (ω0,j (0)) + Cj + Dj . These values are summarized in Table 1. We observe that the FWHM has a strong polarization dependence. The table shows that all the active mode linewidths exhibit a nonlinear temperature dependence, and remain almost constant at low temperature. In some cases the FWHM value increases gradually as the temperature increases, which is due to the fact that Γ(ωj (0), T ) is inversely proportional to the phonon lifetime. More surprising, the FWHM becomes non-monotonic regardless of the polarization, sample thickness, and nature of the mode. One possible explanation for this effect is the role played by other scattering mechanisms taking place at different high-symmetry points that dominate the phonon-phonon coupling effect. As discussed previously for the

16


Page 17 of 40 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

ACS Nano

frequency shift, the ratio Cj /Dj is also important. In some cases, such as for the A2g mode (ZZ-polarization, 5 nm), phonon-phonon scattering is dominated by four-phonon processes, which is in agreement with our previous observations on frequency shift. Table 1: Harmonic frequency, ω0,j (0) (cm−1 ) and linewidth, Γ0 (cm−1 ), as well as anharmonic 0 (0), Cj , Dj and Γ00 (cm−1 ) of few-layer black phosphorus for different constants, Aj , Bj , ω0,j polarization angles. These values are obtained from a temperature-dependent analysis of Raman mode positions (fitted with Eq. 8) and full-width-at-half-maxima (fitted with Eq. 13). Only the Raman active modes at the Γ-point, A1g , B2g and A2g are taken into consideration for two thicknesses 5 nm and 20 nm (in parentheses). Varied incident polarization directions are used: zigzag (ZZ) direction, armchair (AC) direction and 45◦ between AC and ZZ directions. Available values from the literature are reported. Mode Polarization ω0,j (0) ZZ 361.6 A1g (360.4) ◦ 45 -AC/ZZ 362.4 (362.0) AC 361.6 (361.3) Ref. 8 366.5 365.8 26 Ref. 365.5 Ref. 22 366.1 ◦ B2g 45 -AC/ZZ 440.4 (440.5) Ref. 8 448.9 447.0 Ref. 26 443.7 22 Ref. 442.2 2 ZZ 468.5 Ag (465.2) ◦ 45 -AC/ZZ 469.7 (468.2) AC 469.0 (467.4) Ref. 8 476.4 475.9 26 Ref. 471.6 Ref. 22 474.0

0 (0) ω0,j 362.0 (360.5) 362.3 (361.7) 361.8 (361.6) 365.5 365.6 365.3 363.8 439.7 (438.9) 444.6 444.4 443.5 440.9 468.3 (464.6) 468.4 (466.3) 467.8 (466.5) 472.9 472.3 471.2 470.5

Aj 0.433 (0.201) -0.080 (-0.271) 0.233 (0.359) -0.81 -0.28 -0.090 -2.088 -0.569 (-1.487) -4.21 -2.29 0.045 -0.719 0.041 (-0.378) -1.221 (-1.784) -1.148 (-0.662) -3.26 -3.37 -0.078 -3.051

17

Bj -0.081 (-0.100) -0.012 (-0.052) -0.068 (-0.100) -0.188 0.035 -0.103 -0.250 -0.122 (-0.097) -0.06 -0.36 -0.538 -0.268 (-0.247) -0.061 (-0.114) -0.094 (-0.227) -0.28 -0.26 -0.325 -0.400


Γ0 3.723 (5.047) 6.287 (2.323) 6.764 (2.744)

Γ00 3.915 (4.380) 6.376 (2.677) 7.166 (2.654)

4.526 (3.118)

4.312 -0.316 (3.701) (0.530)

4.637 4.743 (11.570) (6.645) 5.093 5.076 (2.642) (3.197) 1.632 4.491 (3.490) (3.140) 0.51

Cj 0.174 (-0.800) 0.115 (0.351) 0.467 (-0.152)

0.033 (-5.741) -0.077 (0.513) 3.330 (-0.504) 4.39

Dj 0.018 (0.132) -0.026 (0.003) -0.064 (0.062)

0.102 (0.054)

0.073 (0.816) 0.060 (0.041) -0.471 (0.154)

ACS Nano 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

Calculating anharmonicity from ab initio molecular dynamics simulations So far, we have used a decay model to interpret Raman spectroscopy data and to account for the temperature dependence of both frequency shift and line broadening. Theoretically, to analyze in detail the anharmonic contributions, we performed AIMD calculations on freestanding monolayer BP for temperature values ranging from 4 to 300 K. The normal mode frequencies were calculated by identifying peak positions from the magnitude of the Fourier transform of the total velocity autocorrelation (see the computational details section). 48 An example of a full power spectrum calculated at room temperature is given in Fig. S4 in the ESI. To take into account the effect of thermal expansion of the 2D material, we first define both the surfacic thermal expansion, αS (T ), and the microscopic Gr¨ uneisen parameters, γj,S (0), based on the variables previously introduced in Eqs. 11 and 12, respectively. Simply, to complete this step, we no longer consider the volume of the material but its in-plane area. αS (T ) was calculated in our previous work. 37 We found γj,S (0) coefficients equal to 0.510, 0.514, and 0.457 at the Γ-point for the A1g , B2g , and A2g phonons of the freestanding monolayer BP, respectively. These values are close to those obtained for bulk BP. 49 The power spectra, P (ω), calculated for freestanding monolayer BP (including thermal expansion) are shown in Fig. 5. Only the three Raman active modes, A1g , B2g , and A2g are represented. The peaks were fitted by Lorentzian functions and the peak intensities were normalized relative to the A2g mode. In these calculations the intensity does not correspond to a measurable physical property. However, P (ω) is an excellent tool for predicting the frequency position of each mode, including anharmonic effects. The harmonic frequencies, calculated at zero temperature based on the finite-displacement and supercell methods, 29,30 are marked by the dashed black line. For comparison, we also performed all these calculations (static and dynamic) without including thermal expansion (a(T ) = 4.509 ˚ A and b(T ) = 3.305 ˚ A). In both cases, with and without thermal expansion effects, we note that temperature induces a redshift in frequency for each mode, which is in agreement with the 18


Page 18 of 40

Page 19 of 40

1.0

A1g

B2g

A 2g

0.8 Intensity (a.u.)

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

ACS Nano

300 K 250 K

0.6

200 K

0.4

150 K 100 K

0.2 0.0 300

50 K 4K

350

400 -1 Frequency (cm )

450

500

Figure 5: Power spectra of freestanding monolayer black phosphorus calculated in the temperature range of 4-300 K, including thermal expansion. Dashed lines correspond to the calculated harmonic frequencies. experimental data where both effects are included. To characterize in detail the harmonic, quasiharmonic, and anharmonic contributions that influence the frequency values, the temperature dependence of the extracted frequency shifts from the power spectra are presented in Figs. 6(a-b). In the first case (Fig. 6(a)), only the anharmonic effects are considered. To fit the extracted data (black dots), we use Eq. 8 taking into account only the first two terms, i.e., the harmonic frequency and the Balkanski expression (Eq. 9). The calculated ω0,j (0) at the Γ-point and fitting parameters Aj and Bj are listed in Table 2. From the fitting results, regardless of the temperature we observe for the A1g and B2g modes that the cubic contribution to anharmonicity dominates, while for the A2g mode the quartic term is the most important. In the case of the A2g mode, we note that the cubic contribution remains constant. When thermal expansion is included, the frequency shift is fitted by adding contributions from Eq. 9 and Eq. 10 (see Fig. 6(b)). For all modes, we find that the quasiharmonic 19


ACS Nano

(a)

w/o Thermal Expansion 445

-1

Frequency (cm )

2 443 A g

441 419 417 B 2g 415 350 348

A1g

346 0

(b)

-1

Frequency (cm )

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 20 of 40

50

100 200 150 Temperature (K)

250

300

250

300

Thermal Expansion 446 444 442 A 2 g 440 418 416 414 B 2g 412 350 348 346 0

A1g 50

100 200 150 Temperature (K)

Figure 6: Temperature dependence of frequency shift for freestanding monolayer black phosphorus, (a) without and (b) including the thermal expansion. The black dots and the red squares represent the calculated data extracted from power spectra and the calculated harmonic frequencies, respectively. The solid black line is a theoretical fit including the three-phonon, four-phonon and thermal expansion contributions, shown as dashed blue, dash/dotted green and dotted red lines, respectively. contribution is negligible, but it indirectly modifies the contribution of the anharmonic effects, in agreement with the computational study presented in Fig. 2. For the A1g mode, the temperature dependence is found to be mainly due to three-phonon processes. Interestingly, the four-phonon processes and the effect of thermal expansion vary linearly with

20


Page 21 of 40 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

ACS Nano

Table 2: Calculated harmonic frequencies, ω0,j (0) (cm−1 ), at the Brillouin zone center, and 0 (0) (cm−1 ), of single-layer black phosphorus. anharmonic frequencies at zero temperature, ω0,j −1 Fitted anharmonic parameters Aj and Bj (cm ) for the temperature dependent frequency shift. Values in parentheses include effect of thermal expansion. A1g B2g A2g

ω0,j (0) 349.2 (349.5) 419.7 (417.8) 445.4 (445.2)

0 (0) ω0,j 348.3 (348.2) 418.0 (416.2) 444.7 (443.8)

Aj -0.882 (-1.267) -1.533 (-1.241) -0.142 (-0.923)

Bj -0.057 (-0.007) -0.165 (-0.392) -0.592 (-0.478)

temperature and they are found to cancel each other. In the case of the B2g and A2g modes, a crossing of the curves around 272 and 166 K, respectively, is observed between the threeand four-phonon processes. At high temperature, the quartic phonon decay process provides the largest contribution. At zero temperature the harmonic effects dominate, however a discrepancy of about 1.5 cm−1 appears between the harmonic frequency (red squares) and the fit including all anharmonic contributions. This gap would be reduced if we included the anharmonic effects near 0 K. As already mentioned in the literature, 37 we notice that all calculated harmonic frequencies, ω0,j (0), for BP are underestimated compared to experimental bulk values. However, the anharmonic coefficients, Aj and Bj , are within the range of values obtained from Raman spectroscopy. To verify that the frequency shift results mainly from phonon-phonon coupling, we also calculated the power spectra of the isolated modes (by properly initializing the velocity of each atom) at room temperature (no thermal expansion). By comparing these frequency values with those obtained at zero temperature (purely harmonic), we find that the frequency remains unchanged regardless of the nature of the mode. These calculations are in good agreement with the experimental results and confirm that the theoretical method employed here makes it possible to predict the behavior of the modes of BP at different temperatures.

21


ACS Nano 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

Anharmonicity from many-body perturbation theory The AIMD simulation technique discussed so far is well suited for determining the normal mode frequencies, although at a very high computational cost. To address this issue and enable a more routine use of the proposed approach, we only evaluate the cubic force constants using the finite-displacement method for freestanding monolayer BP and then carry out many-body perturbation theory to get access to the phonon lifetime, τj (q), and thermal conductivity, κ(T ), (and indirectly to the frequency shift). This process is computationally demanding but it is also particularly well-parallelized, thus considerably reducing the computation time compare to AIMD calculations. The thermal conductivity of bulk and few-layer BP has already been investigated theoretically using this method. 50–60 However the reported results differ widely from each other, e.g., κ(300 K) along the ZZ direction varies from 15.33 to 152.7 W/(mK), while along the AC direction ranges from 4.59 to 63.9 W/(mK). We found that to reduce this variability and get precise results comparable to experiments, it is important to perform several preliminary tests to ensure numerical convergence. Based on the single mode relaxation time approximation (RTA), we tested the q-point sampling mesh and the cutoff pair distance (see Fig. S5 in the ESI). The thermal conductivity of monolayer BP at room temperature along the ZZ (κZZ ) and AC (κAC ) directions converge to 71.75 and 31.11 W/(mK), respectively. A ratio κZZ (300 K)/κAC (300 K) ∼ 2.4 is obtained, due to the anisotropy of the material, indicating an enhanced thermal conductivity along the ZZ direction, in agreement with experiment. 11 To gain more information about the contribution of different phonon modes to κ(T ), we extract the phonon lifetime distributions of the optical A1g , B2g , and A2g modes at 300 K, as represented in Fig. 7 (here including thermal expansion). Four high-density regions of phonon modes centered at around 349, 403, 418, and 423 cm−1 with corresponding phonon lifetimes of about 3.3, 2.7, 1.7, and 2.8 ps can be observed. From the DOS shown in Fig. 2, we identify the peaks with the largest amplitudes at high frequency, which are located at frequencies of 347, 400, 416, and 423 cm−1 . These peaks coincide perfectly with the centroid regions in the phonon lifetime distributions. This 22


Page 22 of 40

Page 23 of 40 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

ACS Nano

means that most of the phonon-phonon interactions occur in the dense areas shown in red in Fig. 7. Since they have short lifetimes (∼ 1.7 − 3.3 ps), these phonons play a key role in the degradation of the lattice thermal conductivity. The temperature dependence of the calculated average phonon lifetimes for freestanding monolayer BP including the thermal expansion is depicted in Fig. 8(a). These values are compared with the experimental ones extracted from Fig. 4. As noticed previously, a rise in temperature implies an increase in the scattering processes which results in a decrease in phonon lifetime. We also note that the phonon lifetime varies with the number of layers. We found that at 50 K, the average phonon lifetime value for the A1g mode does not exceed 16.5 ps. This indicates that our AIMD simulations (60 ps) were sufficiently long. Next, we calculated the real part of the phonon SE, ∆(ωj (0), T ), from its imaginary part, Γ(ωj (0), T ), using the Kramers-Kronig relations at the Γ-point for several temperatures (see the computational details section). These values are compared with those obtained exper0 imentally in Fig. 8(b) where ∆ωj (0, T ) = ωj (0, T ) − ω0,j (0) (including the anharmonicity

coefficients at 0 K) is plotted. These results are also compared with those obtained via AIMD, taking into account the three- and four-phonon processes. Compared to our previous AIMD simulations, the many-body perturbation theory method, which only includes three-phonon coupling, gives access to a similar result where the frequencies of the three Raman active modes are found to be redshifted with increasing temperature. Both theoretical methods tend to agree with experimental data obtained by Raman spectroscopy. However, it is important to mention that taking into account 4-phonon processes (see AIMD results) slightly modifies the frequency shift of a freestanding monolayer BP. A1g and A2g modes exhibit a larger frequency downshift while the frequency of the B2g mode is reduced, which directly affects the lattice thermal conductivity. 31

23


ACS Nano

Phonon Lifetime (ps)

6

5

1 g

6

5

5

B2g mode

A mode 56

4 3 2 1 330



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

6


6

5

4

4

33

Page 24 of 40

4

4

3

3

2

2

1

1 390 395350

4

4

3

3

2

2

A 2g mode

2

21

370 330 390 370

330 350 350

-1

1 330

1 400 350 425 390 -1

Phonon Frequency (cm ) 350

1 370 415330

405

Phonon Frequency (cm ) 370

420

370 440

460 390 -1

Phonon Frequency (cm )

390

2 1 Figure 7: Calculated phonon lifetime Phonon Frequency (cm ) distribution of Ag , B2g , and Ag modes at 300 K with thermal expansion. Low (high) density of phonon modes are represented by blue (red) color. -1

Conclusion Using Raman scattering spectroscopy, the vibrational properties of few-layer BP were obtained beyond the harmonic regime, with varied incident polarization and sample thickness. The frequencies of the modes A1g , B2g , and A2g are redshifted, while the linewidths tend to increase with increasing temperature. Applying a decay model including thermal expansion effects, we found that these phenomena are mainly due to phonon-phonon coupling while thermal expansion has only a modest effect on the shift. There is a thickness-dependent high frequency intralayer mode, specifically the out-of-plane A1g mode, which could be used to determine the number of layers at finite temperature. More surprisingly, the polarization of the laser beam strongly influences the FWHM, which sometimes evolves non-monotonically during a temperature increase. We also performed AIMD simulations based on DFT calculations on freestanding single-layer BP. The mode frequencies were calculated by identifying peak positions from the magnitude of the Fourier transform of the total velocity autocorrelation. The extracted cubic and quartic anharmonic terms are influenced by thermal expansion effects, where the four-phonon processes can be dominant at high temperature. These results are also obtained using many-body perturbation theory, giving access to phonon lifetimes but also, based on the phonon self-energy, to the frequency shift. A rise in temperature implies an increase in the scattering processes which result in a decrease in phonon lifetime. 24


Page 25 of 40

Average Phonon Lifetime (ps)

(a)

(b)

-1

∆ω (cm (cm-1)) ∆ω

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

ACS Nano

15 10 5 0 10

A 2g

B2g

5 0 15 10 5 0

00 -2 -2 -4 -4 -6 -6 00 -2 -2 -4 -4 -6 -6 00 -2 -2 -4 -4

A1g

100

200


5 5nm nm- -60° ZZ o 5 5nm nm- -105° 45 - AC/ZZ 5 5nm 150° nm - AC

A 2g

20 nm - 60° 2020nm nm- -105° ZZ o 2020nm nm- -150° 45 - AC/ZZ

B2g

500

20 nm - AC

MD MDcalculations calculations Perturbation Perturbationtheory theory

100 100

A1g

200 300 400 200 300 400 Temperature Temperature(K) (K)

500 500

Figure 8: Temperature dependence of (a) average phonon lifetime and (b) frequency shift, 0 ∆ωj (0, T ) = ωj (0, T ) − ω0,j (0). These calculations include the contribution of the volumedependent quasiharmonicity. In the case of BP, an improvement in thermal conductivity can be obtained by decreasing the anharmonic effects. One practical solution would be to control the number of layers which composes the material to improve the thermal properties of thermoelectric devices or transistors. 25


ACS Nano 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

Methods Sample Preparation Mechanical exfoliation method was used to prepare BP flakes with different thicknesses from bulk BP crystal (purchased from Smart Elements). The BP flakes were further transferred onto a SiO2 /Si substrate from the scotch-tape. Optical microscope was used to determine the location of the interested BP flakes. The thickness of the BP flakes were identified using atomic force microscope (AFM).

Raman Spectroscopy Raman measurements were performed on a Jobin-Yvon Horiba HR800 spectrometer equipped with a He-Ne 633 nm laser line (red wavelength). For the temperature dependence study, we placed the BP sample in a cryostat (THMS600, Linkam Scientific), where the temperature was properly controlled from 78 to 483 K. For the polarization dependence measurement, the sample was placed on a rotation stage and was rotated during the measurement. The laser power was kept at below 1 mW on the sample, in order to avoid thermal heating caused by the laser beam. A 100× magnification objective was used to focus the laser beam into about 1 micrometer. All the spectra parameters were obtained by fitting the peaks using a Lorentzian/Gaussian mixed function.

Computational Details We performed first-principles calculations based on DFT. The structural and vibrational properties, and the AIMD simulations were carried out using the Vienna Ab Initio Simulation Package (VASP). 61–64 Ion cores were modeled with projector augmented wave (PAW) pseudopotentials. 65 The valence 2s and 3p states of phosphorus were treated explicitly. A plane-wave basis energy cutoff of 500 eV and a Gaussian smearing of 0.005 eV were found to yield converged total energy and forces. To include vdW corrections in our calculations, we 26


Page 26 of 40

Page 27 of 40 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

ACS Nano

used the optB86b-vdW scheme. 66,67 This particular choice of exchange-correlation functional is based on a detail comparative study performed in a previous work. 37 After several tests, we added a vacuum region of 12 ˚ A to minimize interactions between periodic images in the z direction. The k-point sampling was based on a Γ-centered grid for all calculations. For the geometric optimization, the primitive cell and all atoms were relaxed to a force cutoff ˚−1 and a (10 × 15 × 1) grid was used. Concerning the AIMD simulations, a of 10−4 eV A (3 × 3 × 1) supercell of freestanding BP and a (3 × 4 × 1) grid were used. A time step of ∆t = 0.5 fs was chosen and the temperature was set from 4 to 300 K using a Nosé-Hoover chain thermostat. 68–70 Starting from the optimized geometry, 65 ps AIMD runs were performed. To ensure proper thermalization, we only used the last 60 ps of the dynamics to perform the analysis. The velocities were evaluated numerically from the atomic positions by approximating the time derivative with a centered-difference stencil method with 7 points. By partitioning the supercell into primitive cells (each indexed by ρ) with corresponding atoms between primitive cells indexed by N , the velocities can be written as vρN (t). To project onto the Γ-point we first carried out the sum:

vN (t) =

X

vρN (t).

(14)

ρ

For each atom in the effective primitive cell, the autocorrelation of the velocity is expressed as: hvN (t) · vN (t + τ )i .

(15)

One way to calculate a vibrational spectrum from the AIMD trajectory is to use the power spectrum: 71 Z P (ω) = m

+∞

−∞

X

hvN (t) · vN (t + τ )i e

N

−iωτ

dτ ,

(16)

where m is the reduced mass of the system. The harmonic and quasiharmonic (thermal expansion contribution) phonon calculations were performed using the Phonopy code. 72 The

27


ACS Nano 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 28 of 40

details of these calculations are described in a previous paper. 37 We calculate the DOS at room temperature with a (100 × 100 × 120) q-point mesh. By considering only three-phonon interactions, under the conservation of momentum and energy, the lattice thermal conductivity tensor, κ(T ), can be calculated by solving the linearized phonon Boltzmann equation within the single mode relaxation time approximation (RTA), as implemented in the Phono3py code: 73 It can be expressed as: 74

κ(T ) =

1 X τλ (T )vλ ⊗ vλ Cλ (T ), Nq V λ

(17)

where Nq is the total number of q-points sampling the BZ, V denotes the unit cell volume, vλ = ∇q ωλ is the group velocity, Cλ (T ) represents the specific heat capacity and τλ (T ) is the phonon lifetime. We note that κ(T ) depends on the phonon mean free path, lλ (T ) = τλ (T )vλ . For clarity, we introduce the letter λ = {j, q}. This equation is only valid for 3D materials and in the case of a 2D material the vacuum must be removed. For this reason, it is necessary to normalize κ(T ) by multiplying by Lz /d, in which Lz and d are the length of the unit cell along the z direction and the thickness of 2D material, respectively. In this work, we set d = 5.3 ˚ A, as being equal to half of Lz for the AB-stacked bulk BP. 37 The phonon lifetime is obtained as the inverse of the phonon linewidth, such as:

τλ (T ) =

1 . 2Γλ (T )

(18)

The phonon self-energy takes a form analogous to Fermi’s golden rule: 73 Γλ (ω) =

18π X 2 0 λ00 | |φ (nλ0 + nλ00 + 1)δ(ω − ωλ0 − ωλ00 )+ −λλ ~2 λ0 λ00

(nλ0 − nλ00 ) δ(ω + ωλ0 − ωλ00 ) − δ(ω − ωλ0 + ωλ00 ) ,

(19)

where nλ represents the phonon equilibrium occupancy and φ−λλ0 λ00 is the strength of interaction between the three phonons λ, λ0 and λ00 involved in the scattering, obtained from the 28


Page 29 of 40 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

ACS Nano

third-order interactomic force constant (IFC). Note that the energy conservation is described by δ(ω ± ωλ0 ± ωλ00 ), where the signs depend on the scattering process. The IFCs were calculated using the finite-displacement and supercell methods. For purpose of accuracy, we used two different supercells to calculate the second-order harmonic and the third-order anharmonic IFCs using (8 × 8 × 1) and (5 × 5 × 1) supercells, respectively. The displacement amplitude was set to 0.03 ˚ A. At this stage, VASP was used to compute the forces with an electronic energy convergence threshold of 10−8 eV and a k-point grid equal to (2 × 2 × 1) and (4 × 4 × 1) for the largest and the smallest supercells, respectively. For the lattice thermal conductivity calculations, a tetrahedron method was used for BZ integration with σ = 0.1 THz. The convergence of κ(300 K) was studied with respect to the q-point sampling mesh as well as the cutoff interaction distance between atoms, as shown in Fig. S5 in the ESI. Based on the results, κ(300 K) values, along the ZZ and AC directions, start to converge from a (30 × 30 × 1) q-point mesh and a cutoff distance of 7.28 ˚ A, which corresponds to the 19th nearest neighbor. After convergence, the average of the ratio κZZ (300 K)/κAC (300 K) is approximately equal to 2.4, in agreement with the experimental data. 11 For the sake of accurancy, we subsequently used a (50 × 50 × 1) grid mesh and set a cutoff pair distance of 7.28 ˚ A. Then the phonon lifetime is determined as stated in Eq. 18. To obtain the real part of phonon self-energy, ∆λ (ω), we used the Kramers-Kronig relation: 2 ∆λ (ω) = P π

∞

Z 0

ω 0 Γλ (ω 0 ) 0 dω , ω 02 − ω 2

(20)

where P denotes the Cauchy principal value. Numerically this equation can be approximated as: 2 ∆λ (ω) = π

Z 0

∞

ω 0 Γλ (ω 0 ) dω 0 , ω 02 − ω 2 + iη

(21)

with η → 0. A Savitzky-Golay filter 75 was applied to the Γλ (ω) distribution in the frequency domain to make the function integrable.

29


ACS Nano 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

Supporting Information The Electronic Supporting Information (ESI) is available free of charge on the ACS Publications website at DOI: . The ESI contains: Harmonic vs anharmonic oscillations, Raman spectroscopy analyses, calculated power spectrum, and lattice thermal conductivity (PDF).

Acknowledgements Part of this work was performed using supercomputing resources provided by the Center for Computational Innovations (CCI) at Rensselaer Polytechnic Institute. DT acknowledges support from the Office of Naval Research. AC was supported by NSF Grant EFRI 2-DARE (EFRI-1542707) and by a DOE SCGSR award. XL thanks the support of Boston University (BU) and the Photonics Center at BU. “This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education (ORISE) for the DOE. ORISE is managed by ORAU under contract number DE-SC0014664. All opinions expressed in this paper are the author’s and do not necessarily reflect the policies and views of DOE, ORAU, or ORISE.”

30


Page 30 of 40

Page 31 of 40 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

ACS Nano

References 1. Tang, H.; Herman, I. P. Raman Microprobe Scattering of Solid Silicon and Germanium at the Melting Temperature. Phys. Rev. B 1991, 43, 2299. 2. Bertheville, B.; Bill, H.; Hagemann, H. Experimental Raman Scattering Investigation of Phonon Anharmonicity Effects in Li2 S. J. Phys. Condens. Matter 1998, 10, 2155. 3. Link, A.; Bitzer, K.; Limmer, W.; Sauer, R.; Kirchner, C.; Schwegler, V.; Kamp, M.; Ebling, D. G.; Benz, K. W. Temperature Dependence of the E2 and A1 (LO) Phonons in GaN and AlN. J. Appl. Phys. 1999, 86, 6256–6260. 4. Chen, Y.; Peng, B.; Wang, B. Raman Spectra and Temperature-Dependent Raman Scattering of Silicon Nanowires. J. Phys. Chem. C 2007, 111, 5855–5858. 5. Lin, J.; Guo, L.; Huang, Q.; Jia, Y.; Li, K.; Lai, X.; Chen, X. Anharmonic Phonon Effects in Raman Spectra of Unsupported Vertical Graphene Sheets. Phys. Rev. B 2011, 83, 125430. 6. Maczka, M.; Ptak, M.; Da Silva, K. P.; Freire, P. T. C.; Hanuza, J. High-Pressure Raman Scattering and an Anharmonicity Study of Multiferroic Wolframite-Type Mn0.97 Fe0.03 WO4 . J. Phys. Condens. Matter 2012, 24, 345403. 7. Zhao, Z.; Elwood, J.; Carpenter, M. A. Phonon Anharmonicity of PdO Studied by Raman Spectrometry. J. Phys. Chem. C 2015, 119, 23094–23102. 8. Lapi´ nska, A.; Taube, A.; Judek, J.; Zdrojek, M. Temperature Evolution of Phonon Properties in Few-Layer Black Phosphorus. J. Phys. Chem. C 2016, 120, 5265–5270. 9. Long, G.; Maryenko, D.; Shen, J.; Xu, S.; Hou, J.; Wu, Z.; Wong, W. K.; Han, T.; Lin, J.; Cai, Y.; Lortz, R.; Wang, N. Achieving Ultrahigh Carrier Mobility in Two-Dimensional Hole Gas of Black Phosphorus. Nano Lett. 2016, 16, 7768–7773.

31


ACS Nano 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

10. Tran, V.; Soklaski, R.; Liang, Y.; Yang, L. Layer-Controlled Band Gap and Anisotropic Excitons in Few-Layer Black Phosphorus. Phys. Rev. B 2014, 89, 235319. 11. Wang, T.; Han, M.; Wang, R.; Yuan, P.; Xu, S.; Wang, X. Characterization of Anisotropic Thermal Conductivity of Suspended nm-Thick Black Phosphorus with Frequency-Resolved Raman Spectroscopy. J. Appl. Phys. 2018, 123, 145104. 12. Nasri, F.; Aissa, M. F. B.; Belmabrouk, H. Nanoheat Conduction Performance of Black Phosphorus Field-Effect Transistor. IEEE Trans. Electron Devices 2017, 64, 2765–2769. 13. Jang, H.; Wood, J. D.; Ryder, C. R.; Hersam, M. C.; Cahill, D. G. Anisotropic Thermal Conductivity of Exfoliated Black Phosphorus. Adv. Mater. 2015, 27, 8017–8022. 14. Lee, S.; Yang, F.; Suh, J.; Yang, S.; Lee, Y.; Li, G.; Choe, H. S.; Suslu, A.; Chen, Y.; Ko, C.; Park, P.; Liu, K.; Li, J.; Hippalgaonkar, K.; Urban, J. J.; Tongay, S.; Wu, J. Anisotropic In-Plane Thermal Conductivity of Black Phosphorus Nanoribbons at Temperatures Higher than 100 K. Nat. Commun. 2015, 6, 8573. 15. Liu, S.; Huo, N.; Gan, S.; Li, Y.; Wei, Z.; Huang, B.; Liu, J.; Li, J.; Chen, H. ThicknessDependent Raman Spectra, Transport Properties and Infrared Photoresponse of FewLayer Black Phosphorus. J. Mater. Chem. C 2015, 3, 10974–10980. 16. Luo, X.; Lu, X.; Koon, G. K. W.; Castro N., A. H.; Ozyilmaz, B.; Xiong, Q.; Quek, S. Y. Large Frequency Change with Thickness in Interlayer Breathing Mode-Significant Interlayer Interactions in Few Layer Black Phosphorus. Nano Lett. 2015, 15, 3931–3938. 17. Ribeiro, H. B.; Pimenta, M. A.; De Matos, C. J. S. Raman Spectroscopy in Black Phosphorus. J. Raman Spectrosc. JRS 2018, 49, 76–90. 18. Ling, X.; Liang, L.; Huang, S.; Puretzky, A. A.; Geohegan, D. B.; Sumpter, B. G.; Kong, J.; Meunier, V.; Dresselhaus, M. S. Low-Frequency Interlayer Breathing Modes in Few-Layer Black Phosphorus. Nano Lett. 2015, 15, 4080–4088. 32


Page 32 of 40

Page 33 of 40 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

ACS Nano

19. Favron, A.; Gaufrès, E.; Fossard, F.; Phaneuf-LHeureux, A.-L.; Tang, N. Y. W.; Lévesque, P. L.; Loiseau, A.; Leonelli, R.; Francoeur, S.; Martel, R. Photooxidation and Quantum Confinement Effects in Exfoliated Black Phosphorus. Nature Mater. 2015, 14, 826. 20. Phaneuf-L'Heureux, A.-L.; Favron, A.; Germain, J.-F.; Lavoie, P.; Desjardins, P.; Leonelli, R.; Martel, R.; Francoeur, S. Polarization-Resolved Raman Study of BulkLike and Davydov-Induced Vibrational Modes of Exfoliated Black Phosphorus. Nano Lett. 2016, 16, 7761–7767. 21. Xia, F.; Wang, H.; Jia, Y. Rediscovering Black Phosphorus as an Anisotropic Layered Material for Optoelectronics and Electronics. Nat. Commun. 2014, 5, 4458. 22. Zhang, S.; Yang, J.; Xu, R.; Wang, F.; Li, W.; Ghufran, M.; Zhang, Y.-W.; Yu, Z. F.; Zhang, G.; Qin, Q.; Lu, Y. Extraordinary Photoluminescence and Strong Temperature/Angle-Dependent Raman Responses in Few-Layer Phosphorene. ACS Nano 2014, 8, 9590–9596. 23. Wu, J.; Mao, N.; Xie, L.; Xu, H.; Zhang, J. Identifying the Crystalline Orientation of Black Phosphorus Using Angle-Resolved Polarized Raman Spectroscopy. Angew. Chem. Int. Ed. 2015, 54, 2366–2369. 24. Kim, J.; Lee, J.-U.; Lee, J.; Park, H. J.; Lee, Z.; Lee, C.; Cheong, H. Anomalous Polarization Dependence of Raman Scattering and Crystallographic Orientation of Black Phosphorus. Nanoscale 2015, 7, 18708–18715. 25. Ribeiro, H. B.; Pimenta, M. A.; De Matos, C. J. S.; Moreira, R. L.; Rodin, A. S.; Zapata, J. D.; De Souza, E. A. T.; Castro N., A. H. Unusual Angular Dependence of the Raman Response in Black Phosphorus. ACS Nano 2015, 9, 4270–4276. 26. Late, D. J. Temperature Dependent Phonon Shifts in Few-Layer Black Phosphorus. ACS Appl. Mater. Interfaces 2015, 7, 5857–5862. 33


ACS Nano 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

27. Luo, Z.; Maassen, J.; Deng, Y.; Du, Y.; Garrelts, R. P.; Lundstrom, M. S.; Peide, D. Y.; Xu, X. Anisotropic In-Plane Thermal Conductivity Observed in Few-Layer Black Phosphorus. Nat. Commun. 2015, 6, 8572. 28. Axler, S.; Bourdon, P.; Wade, R. Harmonic Function Theory; Springer-Verlag New York, 2013; Vol. 137; pp 25–26. 29. Kresse, G.; Furthm¨ uller, J.; Hafner, J. Ab Initio Force Constant Approach to Phonon Dispersion Relations of Diamond and Graphite. EPL 1995, 32, 729. 30. Parlinski, K.; Li, Z. Q.; Kawazoe, Y. First-Principles Determination of the Soft Mode in Cubic ZrO2 . Phys. Rev. Lett. 1997, 78, 4063. 31. Feng, T.; Lindsay, L.; Ruan, X. Four-Phonon Scattering Significantly Reduces Intrinsic Thermal Conductivity of Solids. Phys. Rev. B 2017, 96, 161201. 32. Maradudin, A. A.; Fein, A. E. Scattering of Neutrons by an Anharmonic Crystal. Phys. Rev. 1962, 128, 2589. 33. Ipatova, I. P.; Maradudin, A. A.; Wallis, R. F. Temperature Dependence of the Width of the Fundamental Lattice-Vibration Absorption Peak in Ionic Crystals. II. Approximate Numerical Results. Phys. Rev. 1967, 155, 882. 34. Liao, B.; Zhou, J.; Qiu, B.; Dresselhaus, M. S.; Chen, G. Ab Initio Study of ElectronPhonon Interaction in Phosphorene. Phys. Rev. B 2015, 91, 235419. 35. Pavone, P.; Karch, K.; Sch¨ utt, O.; Strauch, D.; Windl, W.; Giannozzi, P.; Baroni, S. Ab Initio Lattice Dynamics of Diamond. Phys. Rev. B 1993, 48, 3156. 36. Mounet, N.; Marzari, N. First-Principles Determination of the Structural, Vibrational and Thermodynamic Properties of Diamond, Graphite, and Derivatives. Phys. Rev. B 2005, 71, 205214.

34


Page 34 of 40

Page 35 of 40 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

ACS Nano

37. Tristant, D.; Cupo, A.; Meunier, V. Finite Temperature Stability of Single-Layer Black and Blue Phosphorus Adsorbed on Au(111): a First-Principles Study. 2D Mater. 2018, 5, 035044. 38. Hackney, N. W.; Tristant, D.; Cupo, A.; Daniels, C.; Meunier, V. Shell Model Extension to the Valence Force Field: Application to Single-Layer Black Phosphorus. Phys. Chem. Chem. Phys. 2019, 21, 322–328. 39. Tristant, D.; Zubair, A.; Puech, P.; Neumayer, F.; Moyano, S.; Headrick, R. J.; Tsentalovich, D. E.; Young, C. C.; Gerber, I. C.; Pasquali, M.; Kono, J.; Leotin, J. Enlightening the Ultrahigh Electrical Conductivities of Doped Double-Wall Carbon Nanotube Fibers by Raman Spectroscopy and First-Principles Calculations. Nanoscale 2016, 8, 19668–19676. 40. Zubair, A.; Tristant, D.; Nie, C.; Tsentalovich, D. E.; Headrick, R. J.; Pasquali, M.; Kono, J.; Meunier, V.; Flahaut, E.; Monthioux, M.; Gerber, I. C.; Puech, P. Charged Iodide in Chains Behind the Highly Efficient Iodine Doping in Carbon Nanotubes. Phys. Rev. Mater. 2017, 1, 064002. 41. Balkanski, M.; Wallis, R. F.; Haro, E. Anharmonic Effects in Light Scattering due to Optical Phonons in Silicon. Phys. Rev. B 1983, 28, 1928. 42. Perakis, A.; Sarantopoulou, E.; Raptis, Y. S.; Raptis, C. Temperature Dependence of Raman Scattering and Anharmonicity Study of MgF2 . Phys. Rev. B 1999, 59, 775. 43. Joseph, B.; Demitri, N.; Lotti, P.; Lausi, A.; Dore, P. Unraveling the Peculiarities in the Temperature-Dependent Structural Evolution of Black Phosphorus. Condens. Matter 2017, 2, 11. 44. Xiao, G.; Cao, Y.; Qi, G.; Wang, L.; Zeng, Q.; Liu, C.; Ma, Z.; Wang, K.; Yang, X.; Sui, Y.; Zheng, W.; Zou, B. Compressed Few-Layer Black Phosphorus Nanosheets from

35


ACS Nano 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

Semiconducting to Metallic Transition with the Highest Symmetry. Nanoscale 2017, 9, 10741–10749. 45. Ridley, B. K. The LO Phonon Lifetime in GaN. J. Phys. Condens. Matter 1996, 8, L511. 46. Pei, J.; Gai, X.; Yang, J.; Wang, X.; Yu, Z.; Choi, D.-Y.; Luther-Davies, B.; Lu, Y. Producing Air-Stable Monolayers of Phosphorene and their Defect Engineering. Nat. Commun. 2016, 7, 10450. 47. Korolkov, V. V.; Timokhin, I. G.; Haubrichs, R.; Smith, E. F.; Yang, L.; Yang, S.; Champness, N. R.; Schröder, M.; Beton, P. H. Supramolecular Networks Stabilise and Functionalise Black Phosphorus. Nat. Commun. 2017, 8, 1385. 48. Cupo, A.; Tristant, D.; Rego, K.; Meunier, V. Theoretical Analysis of Spectral Lineshapes From Molecular Dynamics. Npj Comput. Mater. 2019, 5, 82. 49. Villegas, C. E. P.; Rocha, A. R.; Marini, A. Anomalous Temperature Dependence of the Band Gap in Black Phosphorus. Nano Lett. 2016, 16, 5095–5101. 50. Zhu, L.; Zhang, G.; Li, B. Coexistence of Size-Dependent and Size-Independent Thermal Conductivities in Phosphorene. Phys. Rev. B 2014, 90, 214302. 51. 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. PCCP 2015, 17, 4854–4858. 52. Jain, A.; M., A. J. H. Strongly Anisotropic In-Plane Thermal Transport in Single-Layer Black Phosphorene. Sci. Rep. 2015, 5, 8501. 53. Qin, G.; Zhang, X.; Yue, S.-Y.; Qin, Z.; Wang, H.; Han, Y.; Hu, M. Resonant Bonding Driven Giant Phonon Anharmonicity and Low Thermal Conductivity of Phosphorene. Phys. Rev. B 2016, 94, 165445. 36


Page 36 of 40

Page 37 of 40 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

ACS Nano

54. Smith, B.; Vermeersch, B.; Carrete, J.; Ou, E.; Kim, J.; Mingo, N.; Akinwande, D.; Shi, L. Temperature and Thickness Dependences of the Anisotropic In-Plane Thermal Conductivity of Black Phosphorus. Adv. Mater. 2017, 29, 1603756. 55. Zhang, J.; Liu, H. J.; Cheng, L.; Wei, J.; Liang, J. H.; Fan, D. D.; Jiang, P. H.; Shi, J. Thermal Conductivities of Phosphorene Allotropes from First-Principles Calculations: a Comparative Study. Sci. Rep. 2017, 7, 4623. 56. Zhu, J.; Park, H.; Chen, J.-Y.; Gu, X.; Zhang, H.; Karthikeyan, S.; Wendel, N.; Campbell, S. A.; Dawber, M.; Du, X.; Li, M.; Wang, J.; Yang, R.; Wang, X. Revealing the Origins of 3D Anisotropic Thermal Conductivities of Black Phosphorus. Adv. Electron. Mater. 2016, 2, 1600040. 57. Fei, R.; Faghaninia, A.; Soklaski, R.; Yan, J.-A.; Lo, C.; Yang, L. Enhanced Thermoelectric Efficiency via Orthogonal Electrical and Thermal Conductances in Phosphorene. Nano Lett. 2014, 14, 6393–6399. 58. Hong, Y.; Zhang, J.; Huang, X.; Zeng, X. C. Thermal Conductivity of a Two-Dimensional Phosphorene Sheet: a Comparative Study with Graphene. Nanoscale 2015, 7, 18716– 18724. 59. Xu, W.; Zhu, L.; Cai, Y.; Zhang, G.; Li, B. Direction Dependent Thermal Conductivity of Monolayer Phosphorene: Parameterization of Stillinger-Weber Potential and Molecular Dynamics Study. J. Appl. Phys. 2015, 117, 214308. 60. Zhang, Y.-Y.; Pei, Q.-X.; Jiang, J.-W.; Wei, N.; Zhang, Y.-W. Thermal Conductivities of Single-and Multi-Layer Phosphorene: a Molecular Dynamics Study. Nanoscale 2016, 8, 483–491. 61. Kresse, G.; Hafner, J. Ab-Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561.

37


ACS Nano 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

62. Kresse, G.; Hafner, J. Ab-Initio Molecular-Dynamics Simulation of the Liquid-MetalAmorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251– 14269. 63. Kresse, G.; Furthm¨ uller, J. Efficient Iterative Schemes for Ab-Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. 64. Kresse, G.; Furthm¨ uller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15 – 50. 65. Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. 66. Klimeˇs, J.; Bowler, D. R.; Michaelides, A. Chemical Accuracy for the Van der Waals Density Functional. J. Phys. Condens. Matter 2009, 22, 022201. 67. Klimeˇs, J.; Bowler, D. R.; Michaelides, A. Van der Waals Density Functionals Applied to Solids. Phys. Rev. B 2011, 83, 195131. 68. Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511–519. 69. Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. 70. Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nosé–Hoover Chains: the Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97, 2635–2643. 71. Wendler, K.; Brehm, M.; Malberg, F.; Kirchner, B.; Delle Site, L. Short Time Dynamics of Ionic Liquids in AIMD-Based Power Spectra. J. Chem. Theory Comput. 2012, 8, 1570–1579. 72. Togo, A.; Tanaka, I. First Principles Phonon Calculations in Materials Science. Scr. Mater. 2015, 108, 1–5. 38


Page 38 of 40

Page 39 of 40 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

ACS Nano

73. Togo, A.; Chaput, L.; Tanaka, I. Distributions of Phonon Lifetimes in Brillouin Zones. Phys. Rev. B 2015, 91, 094306. 74. Srivastava, G. P. The Physics of Phonons; Taylor & Francis, New York, 1990; pp 122– 136. 75. Savitzky, A.; Golay, M. J. E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Anal. Chem. 1964, 36, 1627–1639.

39


ACS Nano 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

Graphical TOC Entry 3 Phonon Processes

4 Phonon Processes

40


Page 40 of 40

Phonon Anharmonicity in Few-Layer Black Phosphorus

Recommend Documents