Mapping Phonon Modes from Reduced-Dimensional to Bulk Systems
MMapping Phonon Modes fromReduced-Dimensional to Bulk Systems
Hyun-Young Kim, Kevin D. Parrish, and Alan J. H. McGaughey ∗ Department of Mechanical Engineering,Carnegie Mellon University, Pittsburgh, PA 15213 USA (Dated: July 5, 2019)
Abstract
An algorithm for mapping the true phonon modes of a film, which are defined by a two-dimensional (2D) Brillouin zone, to the modes of the corresponding bulk material, which aredefined by a three-dimensional (3D) Brillouin zone, is proposed. The algorithm is based on normalmode decomposition and is inspired by the observation that the atomic motions generated by the2D eigenvectors lead to standing-wave-like behaviors in the cross-plane direction. It is applied tofilms between two and ten unit cells thick built from Lennard-Jones (LJ) argon, whose bulk isisotropic, and graphene, whose bulk (graphite) is anisotropic. For LJ argon, the density of statesdeviates from that of the bulk as the film gets thinner due to phonon frequencies that shift tolower values. This shift is a result of transverse branch splitting due to the film’s anisotropy andthe emergence of a quadratic acoustic branch. As such, while the mapping algorithm works wellfor the thicker LJ argon films, it does not perform as well for the thinner films as there is a weakercorrespondence between the 2D and 3D modes. For graphene, the density of states of even thethinnest films closely matches that of graphite due to the inherent anisotropy, except for a smallshift at low frequency. As a result, the mapping algorithm works well for all thicknesses of thegraphene films, indicating a strong correspondence between the 2D and 3D modes. ∗ [email protected] a r X i v : . [ c ond - m a t . m e s - h a ll ] J u l . INTRODUCTION Phonons are quanta of energy associated with the atomic vibrations in a crystalline solid.They provide a fundamental framework that is used for predicting and understanding avariety of phenomena, including thermal energy storage, thermal expansion, and thermaltransport [1, 2]. Every phonon mode is defined by a wave vector, which specifies its locationin reciprocal space [i.e., the Brillouin zone (BZ)], polarization vector, and frequency. Thephonon modes represent a Fourier transform of the atomic positions to a coordinate systemcomposed of independent harmonic oscillators [3].Advances in fabrication technologies have enabled the nanostructuring of three-dimensional(3D) bulk materials into lower-dimensional objects such as films [two-dimensional (2D)],wires [one-dimensional (1D)], and dots [zero-dimensional (0D)] [4–6]. Changes in the ma-terial dimensionality also affect the dimensionality of the BZ, thereby limiting the allowedwave vectors. This reduction in BZ dimensionality has been exploited to manipulate theproperties of electrons in applications such as field-effect transistors, solar cells, and lasers[7, 8].Our objective is to probe how the dimensionality change impacts phonons in films, whichare periodic in their two in-plane directions. For example, it is known that as a film’sthickness is decreased, phonons begin to scatter with its boundaries in addition to intrin-sic scattering mechanisms (i.e., with other phonons), which reduces their mean free paths[9–12]. Films also limit the allowed cross-plane wave vectors [10, 13] and generate surfacephonons [14]. In contrast, flexural modes in single-layer graphene have a dramatically re-duced scattering phase space compared to multi-layer graphene due to the emergence of newsymmetries that limit their interactions [15, 16].The phonon properties of a film can be predicted using the slab method, which uses a 2DBZ with a unit cell that spans the film thickness, as shown in Fig. 1 [14]. The slab methodis commonly used to model surface phonons [14] and to predict the phonon properties of2D materials such as multi-layer graphene [17]. The slab method stands in contrast to acommon practice of predicting film phonon properties by using the modes corresponding tothe 3D bulk unit cell (see Fig. 1) and assumptions related to the allowed cross-plane wavevectors and/or boundary scattering [10, 11].While the 2D slab method is rigorous, modeling phonon scattering in this framework2
D Treatment 2D Treatment z (Cross-plane) x (In-plane) FIG. 1: Model thin film. The horizontal lines denote the film boundaries. The 3Dtreatment (red) uses the bulk unit cell, which allows cross-plane wave vectors. The 2Dtreatment (blue) has a unit cell that spans the film thickness, which disallows cross-planewave vectors but creates additional polarizations.is computationally challenging. For example, if three-phonon scattering events are to beincluded, the required calculations scale with the fourth power of the number of atomsin the unit cell. Furthermore, any required ab initio calculations with density functionaltheory will scale with the third power of the number of atoms in the computational cell.This challenge leads to the question of when an approximate 3D representation is sufficientto describe the phonon modes in a film.We will investigate this question by developing and applying a mapping algorithm thatlinks the 2D and 3D modes. The majority of previous work on mapping phonon modes,which is also called ”unfolding,” has focused on the link between an alloy and its virtualcrystal representation [18–20], both of which are 3D systems. Allen et al. developed a generalunfolding algorithm that they applied to (i) phonons in a diatomic 1D chain described byan exact two-atom basis and an approximate one-atom basis (i.e., the corresponding virtualcrystal) and (ii) electrons in structures that are finite in one dimension (i.e., slabs) [21]. Theirslab calculation, however, requires the assumption of an artificial periodicity in the cross-plane direction to accommodate their unfolding formulation, which is periodic in nature. Ourproposed mapping algorithm eliminates the cross-plane periodicity requirement by searchingfor signatures of 3D phonon modes in the 2D phonon mode eigenvectors using normal modedecomposition [22, 23].The proposed algorithm is applied to films of Lennard-Jones (LJ) argon (whose bulkform is isotropic) and graphene (whose bulk form, graphite, is anisotropic). The materialdetails, underlying harmonic lattice dynamics theory, and mapping algorithm are described3n Sec. II. The 2D and 3D density of states (DOS) and heat capacities are compared inSec. III A. In Sec. III B, the 2D film dispersions are mapped onto the 3D bulk dispersionsto assess the algorithm and to identify similarities and differences between the 2D and3D modes. The differences for the LJ argon system reveal a transition from an isotropicbulk material to an anisotropic film. This transition results in the separation of pairs ofdegenerate transverse branches into a cross-plane flexural branch and an in-plane transversebranch. This separation is naturally present in graphene, which is inherently anisotropic.
II. METHODOLOGYA. Materials
We model argon using the 12-6 LJ potential φ ( r ij ) = 4 (cid:15) LJ (cid:34)(cid:18) σ LJ r ij (cid:19) − (cid:18) σ LJ r ij (cid:19) (cid:35) , (1)where φ ( r ij ) is the potential energy between atoms i and j , r ij is the distance betweenthem, and σ LJ and (cid:15) LJ are the LJ length and energy scales ( σ LJ = 3 . × − m and (cid:15) LJ =1 . × − J [24]). A cutoff of 2 . σ LJ is applied to limit the interaction range. Bulk argonis a face-centered cubic crystal with a one-atom primitive unit cell and a zero-temperature,zero-pressure lattice constant of 5 .
269 ˚A [25]. We use the four-atom conventional unit celldue to its simple-cubic structure, which has three orthogonal lattice vectors. This choiceresults in the cross-plane direction of the films aligning with one of the lattice vectors and theother two lattice vectors lying in the same plane as the in-plane directions. The alignment ofthe lattice vectors and film directions allows for straightforward visualizations and matchesthe setup of the graphene system. The in-plane directions are labeled with the Cartesiancoordinates x and y , and the cross-plane direction is labeled as z . The lattice vectors, basisvectors, and unit cell are shown in Table I. In subsequent sections, the LJ argon films areidentified by their thickness in bulk conventional unit cells. An N -unit cell film thereforecontains 4 N atoms in its unit cell.The atomic interactions within a graphene layer are modeled using the optimized Tersoffpotential developed by Lindsay and Broido [26]. The in-plane carbon atoms take on ahoneycomb structure with lattice constant of 2 .
50 ˚A [26]. The interlayer coupling in multi-4ayer graphene and graphite is modeled by a 12-6 LJ potential with σ LJ = 3 . × − mand (cid:15) LJ = 7 . × − J [17]. The layer separation is 3 .
35 ˚A [17] and a cutoff of 6 . N -unit cell graphene film thereforecontains 4 N atoms in its unit cell. This choice of unit cell allows the same labeling schemefor the in-plane and cross-plane directions as for LJ argon. For both the LJ argon andmulti-layer graphene films, a vacuum region with a size greater than the interaction cutoffis placed above and below the top and bottom surfaces.The mapping algorithm (Sec. II C) was designed with the assumption of a uniform cross-plane layer separation. Relaxing a film results in a non-uniform layer separation due tothe non-symmetrical forces near its surfaces [27]. We determined the variation in layerseparation by relaxing films using energy minimization. For the five unit cell LJ argon film,the separation of the two outermost layers is 2% larger than that at the film center. For thefive unit cell graphene film, the separation of the two outermost layers is 1% smaller thanthat at the film center. We believe that these variations are small enough that the impacton the phonon modes will not be significant. As such, the LJ argon and graphene films werenot relaxed. B. Harmonic Lattice Dynamics
The frequencies ω and polarization vectors e of the phonon modes with the wave vector κκκ can be found by solving the eigenvalue problem ω ( κκκ, ν ) e ( κκκ, ν ) = D ( κκκ ) e ( κκκ, ν ) , (2)where ν labels the polarization [2]. D is the dynamical matrix, whose entries are D j − α, j (cid:48) − β ( jj (cid:48) , κκκ ) =1( m j m j (cid:48) ) / (cid:88) k (cid:48) Φ αβ ( jk ; j (cid:48) k (cid:48) )exp( iκκκ · [ r ( j (cid:48) k (cid:48) ) − r ( jk )]) . (3)Here, m is the atomic mass and r is the equilibrium position of basis atom j in unit cell k , with k being the central unit cell. The Φ αβ terms are the second-order (i.e., harmonic)5ABLE I: Crystal structures for LJ argon (conventional unit cell) and graphite. LJ Argon Graphene and GraphiteLattice vectors v = [ a, , v = [0 , a, v = [ a, , v = [ a , a √ , v = [0 , , a ] v = [0 , , c ]Basis vectors b = [0 , , b = [0 , a , a ] b = [0 , , b = [ a , a √ , b = [ a , a , b = [ a , , a ] b = [0 , a √ , c ], b = [ a , a √ , c ]Lattice constant(s) a = 5 .
269 ˚A [28] a = 2 .
50 ˚A [26], c = 6 .
70 ˚A [17]Unit cell v v v axyz v v v acxyzκ z = 0 BZ Slice κ x κ y Γ XMΓ = [0 , , , X = [ πa , , , M = [ πa , πa , κ x κ y Γ KMΓ = [0 , , , K = [ π a , , , M = [ πa , π √ a , force constants, where α and β denote the Cartesian directions [1, 2, and 3 for x , y , and z in Eq. (3)]. They are defined as the second derivative of the system potential energy U withrespect to displacements u , or equivalently as the negative of the first derivative of the force F αjk , as Φ αβ ( jk ; j (cid:48) k (cid:48) ) = ∂U∂u αjk ∂u βj (cid:48) k (cid:48) = − ∂F αjk ∂u βj (cid:48) k (cid:48) . (4)We obtain the harmonic force constants with a four-point central difference formula on theforce. Translational invariance is enforced using a Lagrangian approach [29]. The harmonic6attice dynamics calculation can be performed either on a set of wave vectors evenly dis-tributed in the BZ to obtain the DOS or along high-symmetry directions to obtain thedispersion. C. Mapping Algorithm
Our mapping algorithm makes use of the normal mode decomposition technique, whichmaps atomic trajectories onto phonon (i.e., normal) modes [22, 23]. While the atomicmotions in a 2D film cannot form a traveling wave in the cross-plane direction due to a lackof periodicity, they may form standing waves. These standing waves can then be interpretedas a superposition of two traveling waves with opposing cross-plane wave vectors that areassociated with 3D phonons.Each 2D phonon mode has a polarization vector that describes the displacements fromequilibrium of the atoms in the unit cell. We calculate the atomic motions of all atoms inthe unit cell, u , at a temporal phase t as u ( κκκ , ν , t ) = Re { e D ( κκκ , ν )exp[ iω ( κκκ , ν ) t ] } . (5)The temporal phase is varied in the range of 0 ≤ t < πω in intervals of π ω to cover a fullperiod. Here, time is not a dynamic variable. Shifting the polarization vector by a temporalphase allows for the full range of atomic motions in that phonon mode to be explored, whichhelps to reduce numerical errors in the mapping algorithm.Normal mode decomposition is then applied to project the atomic displacements at each t onto all 3D wave vectors ( κκκ ) and their polarization vectors ( e ). As we are working inorthogonal lattices and only mapping to determine the cross-plane wave vector component κ z , the in-plane components κ x and κ y are assumed to be the same in the 2D and 3D BZs.The amplitude of the 3D normal mode coordinate q is calculated at each t from q ( κκκ , ν D , ν , t ) = (cid:88) α,j,l (cid:114) m j N u ,α,jl ([ κ ,x , κ ,y ] , ν D , t ) e † ,α,jl ( κκκ , ν )exp( − iκκκ · [ r ( jl ) − r ( jl )]) , (6)where u ,α,jl is the α -component of the displacement of basis atom j in the 2D unit cell layer l . The unit cell layer is defined as the index of the 3D unit cell stacked in the cross-plane7irection in the 2D unit cell, with l referring to the bottom of the 2D unit cell. Similarly, e ,α,jl refers to the component of the polarization vector corresponding to movement ofbasis atom j and unit cell layer l in the α -direction in the 3D description, with the † symbolindicating its complex conjugate.The sampling of κ z is limited to κ z = nπN a cross , n = ± , , , ..., N (7)to only consider modes that can be described by the discrete number of bulk unit cells inthe film thickness. The cross-plane lattice constant a cross is a for argon and c for graphene(Table I). The wavelength, 2 π/κ z , can be at most twice the film thickness (i.e., at least halfof the wavelength must fit into the film). Eq. (7) thus results in N + 1 mapping destinationsfor the N polarizations for a given κκκ . As such, there will not be a complete one-to-onecorrespondence, as will be discussed in Sec. III B.The normal mode coordinate magnitudes are then summed over the entire period foreach 2D mode for every possible 3D cross-plane wave vector and polarization, generatingquantities of the form q ( κκκ , ν , ν ) = (cid:88) t | q ( κκκ , ν , ν , t ) | . (8)For each 2D mode, the summed normal-mode coordinate is plotted versus the cross-planewave vector. The κ z with the maximum q is used as the mapping destination.This application of normal mode decomposition to the 2D atomic trajectories is similarto the Bloch-wave unfolding algorithm described by Allen et al. [21]. Their unfolding ap-proach, however, relies on comparing systems with equal dimensionality (i.e., 3D to 3D).This reliance is due to their unfolding algorithm comparing a polarization vector to trans-lated versions of itself, with periodicity necessary to wrap around the edges of a supercell.Allen et al. mitigated this limitation by introducing a midline to represent the surfaces of afilm, with the polarization vectors reflecting across it with some corrections. This strategyallowed them to unfold electrons in a silicon film to the primitive bulk unit cell. For phonons,however, normal mode decomposition enables a direct comparison between the 2D and 3Dpolarization vectors, such that the midline approximation is unnecessary.8 II. RESULTSA. Phonon Density of States
Before applying the mapping algorithm, we first compare the 2D and 3D DOS for the LJargon [Fig. 2(a)] and graphene [Fig. 2(b)] systems. The in-plane wave vector sampling wasthe same for the 2D and 3D systems, at 100 ×
100 for LJ argon and 44 ×
44 for graphene. The3D BZ was further sampled with 100 cross-plane wave vectors for LJ argon and 44 cross-plane wave vectors for graphite. The DOS were normalized to have the same integrated areato allow for a direct comparison between different thicknesses.For LJ argon, the two unit cell film DOS [Fig. 2(a)] differs significantly from the 3D DOSwith a shift to lower frequencies. The ten unit cell film DOS, in comparison, is similar tothe 3D DOS. We quantify the difference between the 2D and 3D DOS by calculating theEuclidean distance (cid:107)
DOS − (cid:107) between them as (cid:107) DOS − (cid:107) = (cid:115)(cid:88) i ( h D,i − h D,i ) , (9)where i indicates the frequency bin and h i is the DOS for that bin. A hundred uniform binswere used to divide the frequency range. This quantity was calculated for two to ten unit cellLJ argon films and the results are plotted in Fig. 2(c). As the film thickness increases, thedifference between the 2D and 3D DOS steadily decreases. This result suggests that whilethe phonon modes of the two unit cell film may differ from those of the 3D system, at largerthicknesses the 2D and 3D modes may demonstrate a close to one-to-one correspondence.In contrast, the 2D DOS of a graphene film is nearly identical to the bulk graphite DOSeven at two unit cells thick, with small differences only at low frequencies. The cross-planebonding is much weaker than the in-plane bonding. As such, the in-plane forces are notaffected by the removal of layers and the majority of the frequencies are unchanged. Theacoustic flexural phonon modes are the most affected, which is reflected in the deviationof the low frequency region of the DOS where they reside. In contrast, the LJ argon filmsintroduce anisotropy that is not present in the bulk system, which leads to differencesbetween the 2D and 3D phonon modes.The differences between the 2D and 3D DOS can be further examined by calculating the9 a) (b)(c) (d) FIG. 2: (a) DOS of the two unit cell film, ten unit cell film, and bulk LJ argon systems.As a result of anisotropy, the thinner film DOS differs significantly from the bulk DOS. (b)DOS of the two unit cell graphene film and bulk graphite systems. Even at this smallthickness, the two unit cell DOS matches the bulk DOS due to the weak cross-planeinteractions. (c) Difference between the 2D and 3D DOS as a function of film thickness forLJ argon. As the thickness increases, the 2D DOS converges to the 3D DOS. (d) Specificheat as a function of film thickness, calculated at a temperature of 10 K. The 2D and 3Dspecific heats diverge as thickness decreases due to differences in the DOS at lowfrequencies.specific heat, C ph , from [30] C ph = 1 V (cid:88) κκκ,ν k B x e x ( e x − . (10)Here, x = (cid:126) ω ( κκκ,ν ) k B T , T is the temperature, k B is the Boltzmann constant, (cid:126) is the reduced Planckconstant, and V is the crystal volume. The results for both materials at a temperature of 10K are plotted as a function of film thickness in Fig. 2(d), where they are normalized by thecorresponding bulk value. At this low temperature, the specific heat is mainly determinedby the low-frequency modes. Since the films of both materials have DOS that deviate from10ulk at low frequencies, their specific heats should also deviate. This result is seen in bothLJ argon (with a specific heat up to 1.45 times the bulk value) and graphene (with a specificup to 1.75 times the bulk value). B. Phonon Dispersion
1. LJ Argon
The mapping algorithm was applied to LJ argon films with thicknesses of two to ten unitcells. Examples of the q sweep through κ z [Eq. (8)] are plotted in Figs. 3(a)-3(d) for a tenunit cell LJ argon film at the 2D Γ -point (i.e., the center of the BZ, where κ x = κ y = 0) forfour polarizations. For simplicity, we define κ ∗ z = κ z π/a cross . The q sweeps are symmetricalabout κ ∗ z = 0 for each frequency. This result is expected, as propagating cross-plane modescannot exist in a 2D film. Any cross-plane periodicity must be captured by a standingwave built from two propagating waves that are traveling in opposite directions. There areclear maxima in q at each polarization, which are shown with dashed vertical lines and redmarkers. These peaks identify the mapping location in the 3D BZ.To visualize the mapping algorithm, the 2D dispersion for the ten unit cell LJ argon filmwith varying κ x and κ y = 0 (i.e., along the Γ − X direction) is plotted in Fig.!4. Plottedbelow the 2D dispersion are corresponding bulk 3D dispersions at κ ∗ z values of 0.15, 0.30,and 0.40. Portions of the 2D dispersion are highlighted in red, green, and blue. Applyingthe mapping algorithm matches these modes to κ ∗ z values of 0 .
15 (red), 0 .
30 (green), and0 .
40 (blue). The red and blue ΓΓΓ-point modes correspond to the data from Figs. 3(b) and3(d).If there are strong similarities between the 2D and 3D phonon modes, then the mapped2D dispersions should trace out the 3D branches. The results for the ten unit cell LJ argonfilm are shown in Figs. 5(a)-5(c) from [0 , , κ z ] to [ π/a, , κ z ] at κ ∗ z values of 0, 0.05, and0.30. With the exception of κ ∗ z = 0 .
05, the mapped 2D modes mostly trace the 3D branches.This result demonstrates that 2D modes can be mapped to 3D modes at specific cross-planewave vectors. There are some deviations in the acoustic branches for κ ∗ z = 0 and κ ∗ z = 0 . a) (b)(c) (d) FIG. 3: q sweeps for four Γ -point LJ argon modes in a ten unit cell film (i.e., κ x = κ y = 0).The mapping algorithm is run for all 3D polarizations at each κ z value and thenormalization is based on the maximum value. The peaks are marked red and are at (a) | κ ∗ z | = 0, (b) | κ ∗ z | = 0 .
15, (c) | κ ∗ z | = 0 .
50, and (d) | κ ∗ z | = 0 .
40. The peak locations are usedto map between the 2D and 3D modes.branches, such that multiple branches are required to map across the BZ. (ii) A mapped 2Dbranch has frequencies lower than any 3D branch.12IG. 4: Mapping of selected LJ argon phonon modes from the 2D dispersion (top) tomultiple 3D dispersions (bottom) in a 10 unit cell film. The modes from Figs. 3(b) and3(d) are plotted as blue and red squares. 13 a)(b)(c)
FIG. 5: LJ argon dispersions along κ x for a ten unit cell film with κ y = 0 and (a) κ ∗ z = 0,(b) κ ∗ z = 0 .
05, and (c) κ ∗ z = 0 .
30. 14IG. 6: Dispersion of the two unit cell LJ argon film compared to the 3D dispersion at κ ∗ z = 0. The 2D dispersion has a downward frequency shift for all branches and a splittingof the degenerate bulk transverse branches due to its anisotropy.The second discrepancy can be examined by considering the mapped dispersion along κ x for a two unit cell film for κ y = κ z = 0, as shown in Fig. 6. The downward frequency shift ispresent in all branches and can be attributed to two phenomena. First, there are additionalbranches (e.g., the optical branch with a ΓΓΓ-point frequency near 6 Trad/s). Such branchesare a result of the film anisotropy. As the film thickness decreases, the degenerate bulktransverse branches differentiate. The atoms in one branch move in the in-plane direction,while the atoms in the other move in the cross-plane direction. Due to the free surfaces atthe top and bottom of the film, the restoring forces in the in-plane and cross-plane directionswill be different, which leads to the branch splitting. This effect will occur for all degenerateacoustic and optical branches. Second, one of the acoustic branches is quadratic, which isexpected in a two-dimensional system. A similar behavior will be present in the ten unitcell film and leads to a reduction in frequency.The success of the mapping algorithm is also compromised by the existence of surfacephonon modes, which can be identified by their low participation ratios, P . This quantity,which has a maximum value of unity, is a measure of the fraction of atoms that participatein a given mode and is defined as [31] P ( κκκ , ν ) = (cid:88) j,k (cid:34)(cid:88) α e ,α,jk ( κκκ , ν ) e † ,α,jk ( κκκ , ν ) (cid:35) − . (11)The participation ratios for the 10 unit cell LJ argon film at the ΓΓΓ-point are plotted in Fig.15 a)(b) FIG. 7: Participation ratios for all ΓΓΓ-point modes for the 10 unit cell LJ argon film.Surface modes were identified with participation ratio lower than 0.1 after examining theirpolarization vectors. The surface mode with a frequency of 5.72 Trad/s is colored red. (b) q sweep for the surface mode with a frequency of 5.72 Trad/s. There is no dominant peak,resulting in difficulties for the mapping algorithm.7(a), where four modes with P < . q sweep plottedin Fig. 7(b), which shows no distinct peaks.The presence of the above effects at low thicknesses indicates an inherent differencebetween the 2D phonon modes from their 3D counterparts in LJ argon, going further thana simple limitation on the possible cross-plane wave vectors. In such cases, the rigorous 2D16alculation is necessary to obtain accurate phonon properties.There is also a source of numerical error that is separate from the anisotropy and the sur-face modes. For mapping destinations with the smallest cross-plane wave vector componentfor a given thickness [i.e., n = ± n . This behavior is evident in Figs. 5(b) and 5(c), which show the mapped2D dispersions for a ten unit cell film at κ ∗ z = 0 .
05 ( n = 1) and at κ ∗ z = 0 .
30 ( n = 6). The κ ∗ z = 0 .
30 dispersion has a larger number of 2D mapped modes tracing the 3D branches, butthe κ ∗ z = 0 .
05 dispersion has some 3D branches with almost no 2D mapped modes. Thesemissing modes are a result of the difficulty of fitting the largest possible wavelength intothe thickness of the 2D system. There are thus mis-mapped modes at other κ ∗ z because themodes that should have mapped to κ ∗ z = 0 .
05 must end up somewhere. This effect is presentfor all film thicknesses at the smallest allowed wave vector, but the increasing number ofmapping destinations as the film gets thicker means that the overall effect for the entire 2Ddispersion is reduced.
2. Graphene
The isotropic to anisotropic transition in the LJ argon films can be contrasted with multi-layer graphene and graphite, which are both inherently anisotropic. Mapped 2D dispersionsfor the two and five unit cell graphene films are plotted in Figs. 8(a) and 8(b). Only the lowfrequency branches are shown, as higher frequency branches do not vary greatly with respectto κ ∗ z and are difficult to differentiate. Also plotted is the bulk graphite ΓΓΓ − M dispersion at κ ∗ z = 0. The mapped 2D modes in both the two and five unit cell films closely trace the 3Dbranches. The transverse branches are split in both 2D and 3D due to the anisotropy. Asmall downward shift in frequency is only observed in the two unit cell film, consistent withthe DOS [Fig. 2(b)]. The 2D branches tend to completely map the 3D branches, rather thanin the piece-wise manner seen for LJ argon. This result indicates that the shapes of the 2Dbranches are very similar to those of the 3D branches, which can be attributed to the weakinteractions between the graphene layers. 17 a)(b) FIG. 8: Graphene dispersions along ΓΓΓ- M for (a) two unit cell and (b) five unit cell filmswith κ ∗ z = 0. Only the low-frequency modes are shown. Both the two and five unit cellmapped dispersions are in good agreement with the 3D dispersion. IV. SUMMARY
We applied lattice dynamics calculations to study the relationship between phonon modesin corresponding 2D and 3D systems. Specifically, we developed a mapping algorithm basedon normal mode decomposition to link 2D film and 3D bulk dispersions (Sec. II C). LJargon and graphene systems were analyzed due to their respective isotropy and anisotropyin bulk.As shown in Figs. 2(a) and 2(c), the DOS of LJ argon films converges to the bulk DOS asthe thickness increases. At low thicknesses, however, the 2D and 3D DOS diverge, showingthe effect of anisotropy in the films that is not present in bulk. In contrast, the grapheneDOS is very similar to the bulk DOS even at small thicknesses, as shown in Fig. 2(b),which is a result of the inherent anisotropy in the bulk material due to weak inter-layer18nteractions.Upon application of the mapping algorithm, differences between the 2D and 3D phononmodes were observed in LJ argon. At larger thicknesses, the differences are small and the 2Dmodes correspond strongly to 3D modes with specific cross-plane wave-vector components.Some differences, however, exist in the low-frequency acoustic branches [Figs. 5(a) and 5(b)].At smaller thicknesses, there are significant differences such as downward frequency shiftsand branch splitting (Fig. 6). These effects can be attributed to the emergence of anisotropyand surface modes in the LJ argon films. The graphene systems, with their weak inter-layerinteractions, do not show the same deviations. Instead, their 2D phonon modes are stronglymatched to the 3D phonon modes at specific cross-plane wave-vectors, with only smalldownward frequency shifts at low frequencies [Figs. 8(a) and 8(b)].The observed differences between mapped 2D modes and their 3D counterparts are rele-vant to predictions of film thermal conductivity. The majority of previous calculations haveassumed the existence of 3D modes whose mean free paths are modified with a boundaryscattering model [11–13, 32]. In some studies, the number of cross-plane wave vectors waslimited by the film thickness [10, 13]. While a treatment based on 3D modes may be appro-priate for thicker films or films with weak cross-plane interactions, as we have shown here,they will break down in some very thin films due to the emergence of anisotropy and surfacemodes. In such cases, an accurate thermal conductivity prediction requires use of the 2Dslab unit cell.
ACKNOWLEDGEMENTS
This work was supported by NSF Award DMR-1507325.The following article has been submitted to the Journal of Applied Physics. After it ispublished, it will be found at Link. [1] G. Chen,
Nanoscale Energy Transport and Conversion: A Parallel Treatment of Electrons,Molecules, Phonons, and Photons (OUP Oxford, 2005).[2] M. T. Dove,
Introduction to Lattice Dynamics (Cambridge University Press, 2010).
3] J. M. Ziman,
Electrons and Phonons: The Theory of Transport Phenomena in Solids (OUPOxford, 2001).[4] N. Kaiser, Applied Optics , 3053 (2002).[5] Y. Xia, P. Yang, Y. Sun, Y. Wu, B. Mayers, B. Gates, Y. Yin, F. Kim, and H. Yan, AdvancedMaterials , 353 (2003).[6] Y. Pu, F. Cai, D. Wang, J.-X. Wang, and J.-F. Chen, Industrial & Engineering ChemistryResearch , 1790 (2018).[7] D. Bimberg and U. W. Pohl, Materials Today , 388 (2011).[8] N. S. Mohammad, Journal of Physics: Condensed Matter , 423202 (2014).[9] W. Liu and M. Asheghi, Applied Physics Letters , 3819 (2004).[10] J. E. Turney, A. J. H. McGaughey, and C. H. Amon, Journal of Applied Physics , 024317(2010).[11] A. Jain, Y.-J. Yu, and A. J. H. McGaughey, Physical Review B , 195301 (2013).[12] J. Cuffe, J. K. Eliason, A. A. Maznev, K. C. Collins, J. A. Johnson, A. Shchepetov, M. Prun-nila, J. Ahopelto, C. M. Sotomayor Torres, G. Chen, and K. A. Nelson, Physical Review B (2015).[13] X. Wang and B. Huang, Scientific Reports (2015), 10.1038/srep06399.[14] W. Kress and F. W. de Wette, Surface Phonons (Springer Berlin Heidelberg, 1991).[15] L. Lindsay, D. A. Broido, and N. Mingo, Physical Review B (2010).[16] A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau,Nano Letters , 902 (2008).[17] L. Lindsay, D. A. Broido, and N. Mingo, Physical Review B , 235428 (2011).[18] T. B. Boykin, N. Kharche, G. Klimeck, and M. Korkusinski, Journal of Physics: CondensedMatter , 036203 (2007).[19] T. B. Boykin, A. Ajoy, H. Ilatikhameneh, M. Povolotskyi, and G. Klimeck, Physical ReviewB , 205214 (2014).[20] Y. Ikeda, A. Carreras, A. Seko, A. Togo, and I. Tanaka, Physical Review B , 024305 (2017).[21] P. B. Allen, T. Berlijn, D. A. Casavant, and J. M. Soler, Physical Review B , 085322(2013).[22] A. J. C. Ladd, B. Moran, and W. G. Hoover, Physical Review B , 5058 (1986).[23] A. J. H. McGaughey and M. Kaviany, Physical Review B (2004).
24] N. W. Ashcroft and N. D. Mermin,
Solid State Physics (Saunders College Publishing, 1976).[25] A. J. H. McGaughey,
Phonon transport in molecular dynamics simulations: formulation andthermal conductivity prediction , Ph.D. thesis, University of Michigan (2004).[26] L. Lindsay and D. A. Broido, Physical Review B , 205441 (2010).[27] B. Fu, K. D. Parrish, H.-Y. Kim, G. Tang, and A. J. H. McGaughey, arXiv:1907.00969[cond-mat] (2019), arXiv: 1907.00969.[28] A. J. H. McGaughey, M. I. Hussein, E. S. Landry, M. Kaviany, and G. M. Hulbert, PhysicalReview B (2006).[29] W. Li, L. Lindsay, D. A. Broido, D. A. Stewart, and N. Mingo, Physical Review B , 174307(2012).[30] A. Togo, L. Chaput, I. Tanaka, and G. Hug, Physical Review B , 174301 (2010).[31] J. L. Feldman, M. D. Kluge, P. B. Allen, and F. Wooten, Physical Review B , 12589 (1993).[32] B. Fu, G. Tang, and Y. Li, Physical Chemistry Chemical Physics , 28517 (2017)., 28517 (2017).