Mapping Phonon Modes from Reduced-Dimensional to Bulk Systems

TL;DR

This paper introduces an algorithm to map phonon modes from 2D film systems to their 3D bulk counterparts, analyzing its effectiveness across different materials and film thicknesses.

Contribution

The paper presents a novel normal mode decomposition algorithm for mapping phonon modes from reduced-dimensional films to bulk systems, accounting for anisotropy and thickness effects.

Findings

The algorithm performs well for thick LJ argon films but less so for thinner ones.

For graphene, the mapping remains accurate across all film thicknesses due to inherent anisotropy.

Phonon density of states shifts with film thickness, affecting mode correspondence.

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 are defined by a three-dimensional (3D) Brillouin zone, is proposed. The algorithm is based on normal mode decomposition and is inspired by the observation that the atomic motions generated by the 2D eigenvectors lead to standing-wave-like behaviors in the cross-plane direction. It is applied to films between two and ten unit cells thick built from Lennard-Jones (LJ) argon, whose bulk is isotropic, and graphene, whose bulk (graphite) is anisotropic. For LJ argon, the density of states deviates from that of the bulk as the film gets thinner due to phonon frequencies that shift to lower values. This shift is a result of transverse branch splitting due to the film's anisotropy and the emergence of a quadratic acoustic branch. As such, while the mapping algorithm works well for the thicker LJ argon films, it does not perform as well for the thinner films as there is a weaker correspondence between the 2D and 3D modes. For graphene, the density of states of even the thinnest films closely matches that of graphite due to the inherent anisotropy, except for a small shift at low frequency. As a result, the mapping algorithm works well for all thicknesses of the graphene films, indicating a strong correspondence between the 2D and 3D modes.

Tables

Table 1. Table 1 : Crystal structures for LJ argon (conventional unit cell) and graphite.

	LJ Argon	Graphene and Graphite
Lattice vectors	${𝐯}_1=[a,0,0]$, ${𝐯}_2=[0,a,0]$,	${𝐯}_1=[a,0,0]$, ${𝐯}_2=[\frac{a}{2},\frac{a\sqrt{3}}{2},0]$,
	${𝐯}_3=[0,0,a]$	${𝐯}_3=[0,0,c]$
Basis vectors	${𝐛}_1=[0,0,0]$, ${𝐛}_2=[0,\frac{a}{2},\frac{a}{2}]$	${𝐛}_1=[0,0,0]$, ${𝐛}_2=[\frac{a}{2},\frac{a\sqrt{3}}{6},0]$
	${𝐛}_3=[\frac{a}{2},\frac{a}{2},0]$, ${𝐛}_4=[\frac{a}{2},0,\frac{a}{2}]$	${𝐛}_3=[0,\frac{a\sqrt{3}}{3},\frac{c}{2}]$, ${𝐛}_4=[\frac{a}{2},\frac{a\sqrt{3}}{6},\frac{c}{2}]$
Lattice constant(s)	$a=5.269\AA$ McGaughey et al. (2006)	$a=2.50\AA$ Lindsay and Broido (2010), $c=6.70\AA$ Lindsay et al. (2011)
Unit cell
${\kappa }_z=0$ BZ Slice

Equations

$$\varphi (r_{ij})=4{ϵ}_{\mathrm{LJ}}\left[{\left(\frac{{\sigma }_{\mathrm{LJ}}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{\mathrm{LJ}}}{r_{ij}}\right)}^6\right],$$

$${\omega }^2(\mathbf{\kappa },\nu )\mathbf{e}(\mathbf{\kappa },\nu )=\mathbf{D}(\mathbf{\kappa })\mathbf{e}(\mathbf{\kappa },\nu ),$$

$$\begin{array}{rl}{D}_{3(j-1)+\alpha ,3(j^{\prime }-1)+\beta }(j{j}^{\prime },\mathbf{\kappa })& =\\ \frac{1}{(m_j{m}_{j^{\prime }}{)}^{1/2}}& \sum _{k^{\prime }}{\Phi }_{\alpha \beta }(j{k}_0;j^{\prime }{k}^{\prime })\exp (i\mathbf{\kappa }\cdot [\mathbf{r}(j^{\prime }{k}^{\prime })-\mathbf{r}(j{k}_0)]).\end{array}$$

$${\Phi }_{\alpha \beta }(jk;j^{\prime }{k}^{\prime })=\frac{\partial U}{\partial u_{jk}^{\alpha }\partial u_{j^{\prime }{k}^{\prime }}^{\beta }}=-\frac{\partial F_{jk}^{\alpha }}{\partial u_{j^{\prime }{k}^{\prime }}^{\beta }}.$$

$${\mathbf{u}}_{2D}({\mathbf{\kappa }}_{2D},{\nu }_{2D},t)=\mathrm{Re}\{{\mathbf{e}}_{2D}({\mathbf{\kappa }}_{2D},{\nu }_{2D})\exp [i{\omega }_{2D}({\mathbf{\kappa }}_{2D},{\nu }_{2D})t]\}.$$

$$q({\mathbf{\kappa }}_{3D},{\nu }_{3D},{\nu }_{2D},t)=\sum _{\alpha ,j,l}\sqrt{\frac{m_j}{N}}{u}_{2D,\alpha ,jl}([{\kappa }_{3D,x},{\kappa }_{3D,y}],{\nu }_{2D},t)e_{3D,\alpha ,jl}^{†}({\mathbf{\kappa }}_{3D},{\nu }_{3D})\exp (-i{\mathbf{\kappa }}_{3D}\cdot [\mathbf{r}(jl)-\mathbf{r}(j{l}_0)]),$$

$${\kappa }_z=\frac{n\pi }{N{a}_{\mathrm{cross}}},n=\pm 0,1,2,...,N$$

$$\overline{q}({\mathbf{\kappa }}_{3D},{\nu }_{3D},{\nu }_{2D})=\sum _t\left|q({\mathbf{\kappa }}_{3D},{\nu }_{3D},{\nu }_{2D},t)\right|.$$

$$\Vert \mathrm{DO}{\mathrm{S}}_{2D-3D}\Vert =\sqrt{\sum _i(h_{2D,i}-h_{3D,i}{)}^2},$$

$$C_{\mathrm{ph}}=\frac{1}{V}\sum _{\mathbf{\kappa },\nu }\frac{k_{\mathrm{B}}{x}^2{e}^x}{(e^x-1{)}^2}.$$

$$P({\mathbf{\kappa }}_{2D},{\nu }_{2D})={\left(\sum _{j,k}{\left[\sum _{\alpha }{e}_{2D,\alpha ,jk}({\mathbf{\kappa }}_{2D},{\nu }_{2D})e_{2D,\alpha ,jk}^{†}({\mathbf{\kappa }}_{2D},{\nu }_{2D})\right]}^2\right)}^{-1}.$$

Full text

I 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 a variety of phenomena, including thermal energy storage, thermal expansion, and thermal transport Chen (2005); Dove (2010). Every phonon mode is defined by a wave vector, which specifies its location in reciprocal space [i.e., the Brillouin zone (BZ)], polarization vector, and frequency. The phonon modes represent a Fourier transform of the atomic positions to a coordinate system composed of independent harmonic oscillators Ziman (2001).

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)] Kaiser (2002); Xia et al. (2003); Pu et al. (2018). Changes in the material dimensionality also affect the dimensionality of the BZ, thereby limiting the allowed wave vectors. This reduction in BZ dimensionality has been exploited to manipulate the properties of electrons in applications such as field-effect transistors, solar cells, and lasers Bimberg and Pohl (2011); Mohammad (2014).

Our objective is to probe how the dimensionality change impacts phonons in films, which are periodic in their two in-plane directions. For example, it is known that as a film’s thickness is decreased, phonons begin to scatter with its boundaries in addition to intrinsic scattering mechanisms (i.e., with other phonons), which reduces their mean free paths Liu and Asheghi (2004); Turney et al. (2010); Jain et al. (2013); Cuffe et al. (2015). Films also limit the allowed cross-plane wave vectors Turney et al. (2010); Wang and Huang (2015) and generate surface phonons Kress and de Wette (1991). In contrast, flexural modes in single-layer graphene have a dramatically reduced scattering phase space compared to multi-layer graphene due to the emergence of new symmetries that limit their interactions Lindsay et al. (2010); Balandin et al. (2008).

The phonon properties of a film can be predicted using the slab method, which uses a 2D BZ with a unit cell that spans the film thickness, as shown in Fig. 1 Kress and de Wette (1991). The slab method is commonly used to model surface phonons Kress and de Wette (1991) and to predict the phonon properties of 2D materials such as multi-layer graphene Lindsay et al. (2011). The slab method stands in contrast to a common practice of predicting film phonon properties by using the modes corresponding to the 3D bulk unit cell (see Fig. 1) and assumptions related to the allowed cross-plane wave vectors and/or boundary scattering Turney et al. (2010); Jain et al. (2013).

While the 2D slab method is rigorous, modeling phonon scattering in this framework is computationally challenging. For example, if three-phonon scattering events are to be included, the required calculations scale with the fourth power of the number of atoms in the unit cell. Furthermore, any required ab initio calculations with density functional theory 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 sufficient to describe the phonon modes in a film.

We will investigate this question by developing and applying a mapping algorithm that links 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 virtual crystal representation Boykin et al. (2007, 2014); Ikeda et al. (2017), both of which are 3D systems. Allen et al. developed a general unfolding algorithm that they applied to (i) phonons in a diatomic 1D chain described by an exact two-atom basis and an approximate one-atom basis (i.e., the corresponding virtual crystal) and (ii) electrons in structures that are finite in one dimension (i.e., slabs) Allen et al. (2013). Their slab calculation, however, requires the assumption of an artificial periodicity in the cross-plane direction to accommodate their unfolding formulation, which is periodic in nature. Our proposed mapping algorithm eliminates the cross-plane periodicity requirement by searching for signatures of 3D phonon modes in the 2D phonon mode eigenvectors using normal mode decomposition Ladd et al. (1986); McGaughey and Kaviany (2004).

The proposed algorithm is applied to films of Lennard-Jones (LJ) argon (whose bulk form is isotropic) and graphene (whose bulk form, graphite, is anisotropic). The material details, underlying harmonic lattice dynamics theory, and mapping algorithm are described in Sec. II. The 2D and 3D density of states (DOS) and heat capacities are compared in Sec. III.1. In Sec. III.2, the 2D film dispersions are mapped onto the 3D bulk dispersions to assess the algorithm and to identify similarities and differences between the 2D and 3D modes. The differences for the LJ argon system reveal a transition from an isotropic bulk material to an anisotropic film. This transition results in the separation of pairs of degenerate transverse branches into a cross-plane flexural branch and an in-plane transverse branch. This separation is naturally present in graphene, which is inherently anisotropic.

II Methodology

II.1 Materials

We model argon using the 12-6 LJ potential

[TABLE]

where $\phi(r_{ij})$ is the potential energy between atoms $i$ and $j$, $r_{ij}$ is the distance between them, and $\sigma_{\mathrm{LJ}}$ and $\epsilon_{\mathrm{LJ}}$ are the LJ length and energy scales ($\sigma_{\mathrm{LJ}}=3.4\times 10^{-10}\,\mathrm{m}$ and $\epsilon_{\mathrm{LJ}}=1.67\times 10^{-21}\,\mathrm{J}$ Ashcroft and Mermin (1976)). A cutoff of $2.5\sigma_{\mathrm{LJ}}$ is applied to limit the interaction range. Bulk argon is a face-centered cubic crystal with a one-atom primitive unit cell and a zero-temperature, zero-pressure lattice constant of $5.269\,\mathrm{\AA}$ McGaughey (2004). We use the four-atom conventional unit cell due to its simple-cubic structure, which has three orthogonal lattice vectors. This choice results in the cross-plane direction of the films aligning with one of the lattice vectors and the other two lattice vectors lying in the same plane as the in-plane directions. The alignment of the lattice vectors and film directions allows for straightforward visualizations and matches the setup of the graphene system. The in-plane directions are labeled with the Cartesian coordinates $x$ and $y$, and the cross-plane direction is labeled as $z$. The lattice vectors, basis vectors, and unit cell are shown in Table 1. In subsequent sections, the LJ argon films are identified by their thickness in bulk conventional unit cells. An $N$-unit cell film therefore contains $4N$ atoms in its unit cell.

The atomic interactions within a graphene layer are modeled using the optimized Tersoff potential developed by Lindsay and Broido Lindsay and Broido (2010). The in-plane carbon atoms take on a honeycomb structure with lattice constant of $2.50\,\mathrm{\AA}$ Lindsay and Broido (2010). The interlayer coupling in multi-layer graphene and graphite is modeled by a 12-6 LJ potential with $\sigma_{\mathrm{LJ}}=3.276\times 10^{-10}\,\mathrm{m}$ and $\epsilon_{\mathrm{LJ}}=7.37\times 10^{-22}\,\mathrm{J}$ Lindsay et al. (2011). The layer separation is $3.35\,\mathrm{\AA}$ Lindsay et al. (2011) and a cutoff of $6.5\,\mathrm{\AA}$ was applied to ensure that a layer only interacts with its neighboring layers. The layers are stacked in an AB structure that repeats every two layers to form a four-atom unit cell for graphite. The structure is shown in Table 1. The thicknesses of the multi-layer graphene films are defined by the number of bulk unit cells. An $N$-unit cell graphene film therefore contains $4N$ atoms in its unit cell. This choice of unit cell allows the same labeling scheme for the in-plane and cross-plane directions as for LJ argon. For both the LJ argon and multi-layer graphene films, a vacuum region with a size greater than the interaction cutoff is placed above and below the top and bottom surfaces.

The mapping algorithm (Sec. II.3) was designed with the assumption of a uniform cross-plane layer separation. Relaxing a film results in a non-uniform layer separation due to the non-symmetrical forces near its surfaces Fu et al. (2019). We determined the variation in layer separation 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 the five unit cell graphene film, the separation of the two outermost layers is 1% smaller than that at the film center. We believe that these variations are small enough that the impact on the phonon modes will not be significant. As such, the LJ argon and graphene films were not relaxed.

II.2 Harmonic Lattice Dynamics

The frequencies $\omega$ and polarization vectors $\bf{e}$ of the phonon modes with the wave vector $\boldsymbol{\kappa}$ can be found by solving the eigenvalue problem

[TABLE]

where $\nu$ labels the polarization Dove (2010). $\bf{D}$ is the dynamical matrix, whose entries are

[TABLE]

Here, $m$ is the atomic mass and $\mathbf{r}$ is the equilibrium position of basis atom $j$ in unit cell $k$, with $k_{0}$ being the central unit cell. The $\Phi_{\alpha\beta}$ terms are the second-order (i.e., harmonic) force constants, where $\alpha$ and $\beta$ 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$ with respect to displacements $u$, or equivalently as the negative of the first derivative of the force $F^{\alpha}_{jk}$, as

[TABLE]

We obtain the harmonic force constants with a four-point central difference formula on the force. Translational invariance is enforced using a Lagrangian approach Li et al. (2012). The harmonic lattice dynamics calculation can be performed either on a set of wave vectors evenly distributed in the BZ to obtain the DOS or along high-symmetry directions to obtain the dispersion.

II.3 Mapping Algorithm

Our mapping algorithm makes use of the normal mode decomposition technique, which maps atomic trajectories onto phonon (i.e., normal) modes Ladd et al. (1986); McGaughey and Kaviany (2004). While the atomic motions in a 2D film cannot form a traveling wave in the cross-plane direction due to a lack of periodicity, they may form standing waves. These standing waves can then be interpreted as a superposition of two traveling waves with opposing cross-plane wave vectors that are associated with 3D phonons.

Each 2D phonon mode has a polarization vector that describes the displacements from equilibrium of the atoms in the unit cell. We calculate the atomic motions of all atoms in the unit cell, $\mathbf{u}_{\mathrm{2D}}$, at a temporal phase $t$ as

[TABLE]

The temporal phase is varied in the range of $0\leq t<\frac{2\pi}{\omega_{\mathrm{2D}}}$ in intervals of $\frac{\pi}{50\omega_{\mathrm{2D}}}$ to cover a full period. Here, time is not a dynamic variable. Shifting the polarization vector by a temporal phase allows for the full range of atomic motions in that phonon mode to be explored, which helps 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 ($\boldsymbol{\kappa}_{\mathrm{3D}}$) and their polarization vectors ($\mathbf{e}_{\mathrm{3D}}$). As we are working in orthogonal lattices and only mapping to determine the cross-plane wave vector component $\kappa_{z}$, the in-plane components $\kappa_{x}$ and $\kappa_{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

[TABLE]

where $u_{\mathrm{2D},\alpha,jl}$ is the $\alpha$-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-plane direction in the 2D unit cell, with $l_{0}$ referring to the bottom of the 2D unit cell. Similarly, $e_{\mathrm{3D},\alpha,jl}$ refers to the component of the polarization vector corresponding to movement of basis atom $j$ and unit cell layer $l$ in the $\alpha$-direction in the 3D description, with the $\dagger$ symbol indicating its complex conjugate.

The sampling of $\kappa_{z}$ is limited to

[TABLE]

to only consider modes that can be described by the discrete number of bulk unit cells in the film thickness. The cross-plane lattice constant $a_{\mathrm{cross}}$ is $a$ for argon and $c$ for graphene (Table 1). The wavelength, $2\pi/\kappa_{z}$, can be at most twice the film thickness (i.e., at least half of the wavelength must fit into the film). Eq. (7) thus results in $N+1$ mapping destinations for the $N$ polarizations for a given $\boldsymbol{\kappa}_{\mathrm{2D}}$. As such, there will not be a complete one-to-one correspondence, as will be discussed in Sec. III.2.

The normal mode coordinate magnitudes are then summed over the entire period for each 2D mode for every possible 3D cross-plane wave vector and polarization, generating quantities of the form

[TABLE]

For each 2D mode, the summed normal-mode coordinate is plotted versus the cross-plane wave vector. The $\kappa_{z}$ with the maximum $\overline{q}$ is used as the mapping destination.

This application of normal mode decomposition to the 2D atomic trajectories is similar to the Bloch-wave unfolding algorithm described by Allen et al. Allen et al. (2013). Their unfolding approach, 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 translated 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 a film, with the polarization vectors reflecting across it with some corrections. This strategy allowed 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 3D polarization vectors, such that the midline approximation is unnecessary.

III Results

III.1 Phonon Density of States

Before applying the mapping algorithm, we first compare the 2D and 3D DOS for the LJ argon [Fig. 2(a)] and graphene [Fig. 2(b)] systems. The in-plane wave vector sampling was the same for the 2D and 3D systems, at $100\times 100$ for LJ argon and $44\times 44$ for graphene. The 3D 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 area to 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 DOS with a shift to lower frequencies. The ten unit cell film DOS, in comparison, is similar to the 3D DOS. We quantify the difference between the 2D and 3D DOS by calculating the Euclidean distance $\lVert\mathrm{DOS_{2D-3D}}\rVert$ between them as

[TABLE]

where $i$ indicates the frequency bin and $h_{i}$ is the DOS for that bin. A hundred uniform bins were used to divide the frequency range. This quantity was calculated for two to ten unit cell LJ argon films and the results are plotted in Fig. 2(c). As the film thickness increases, the difference between the 2D and 3D DOS steadily decreases. This result suggests that while the phonon modes of the two unit cell film may differ from those of the 3D system, at larger thicknesses 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 DOS even at two unit cells thick, with small differences only at low frequencies. The cross-plane bonding is much weaker than the in-plane bonding. As such, the in-plane forces are not affected by the removal of layers and the majority of the frequencies are unchanged. The acoustic flexural phonon modes are the most affected, which is reflected in the deviation of the low frequency region of the DOS where they reside. In contrast, the LJ argon films introduce anisotropy that is not present in the bulk system, which leads to differences between the 2D and 3D phonon modes.

The differences between the 2D and 3D DOS can be further examined by calculating the specific heat, $C_{\mathrm{ph}}$, from Togo et al. (2010)

[TABLE]

Here, $x=\frac{\hbar\omega(\boldsymbol{\kappa},\nu)}{k_{\mathrm{B}}T}$, $T$ is the temperature, $k_{\mathrm{B}}$ is the Boltzmann constant, $\hbar$ is the reduced Planck constant, and $V$ is the crystal volume. The results for both materials at a temperature of 10 K are plotted as a function of film thickness in Fig. 2(d), where they are normalized by the corresponding bulk value. At this low temperature, the specific heat is mainly determined by the low-frequency modes. Since the films of both materials have DOS that deviate from bulk at low frequencies, their specific heats should also deviate. This result is seen in both LJ argon (with a specific heat up to 1.45 times the bulk value) and graphene (with a specific up to 1.75 times the bulk value).

III.2 Phonon Dispersion

III.2.1 LJ Argon

The mapping algorithm was applied to LJ argon films with thicknesses of two to ten unit cells. Examples of the $\overline{q}$ sweep through $\kappa_{z}$ [Eq. (8)] are plotted in Figs. 3(a)-3(d) for a ten unit cell LJ argon film at the 2D $\mathrm{\mathbf{\Gamma}}$-point (i.e., the center of the BZ, where $\kappa_{x}=\kappa_{y}=0$) for four polarizations. For simplicity, we define $\kappa_{z}^{*}=\frac{\kappa_{z}}{2\pi/a_{cross}}$. The $\overline{q}$ sweeps are symmetrical about $\kappa_{z}^{*}=0$ for each frequency. This result is expected, as propagating cross-plane modes cannot exist in a 2D film. Any cross-plane periodicity must be captured by a standing wave built from two propagating waves that are traveling in opposite directions. There are clear maxima in $\overline{q}$ at each polarization, which are shown with dashed vertical lines and red markers. 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 film with varying $\kappa_{x}$ and $\kappa_{y}=0$ (i.e., along the $\mathbf{\Gamma}-\mathbf{X}$ direction) is plotted in Fig.. Plotted below the 2D dispersion are corresponding bulk 3D dispersions at $\kappa_{z}^{*}$ values of 0.15, 0.30, and 0.40. Portions of the 2D dispersion are highlighted in red, green, and blue. Applying the mapping algorithm matches these modes to $\kappa_{z}^{*}$ values of $0.15$ (red), $0.30$ (green), and $0.40$ (blue). The red and blue $\boldsymbol{\Gamma}$-point modes correspond to the data from Figs. 3(b) and 3(d).

If there are strong similarities between the 2D and 3D phonon modes, then the mapped 2D dispersions should trace out the 3D branches. The results for the ten unit cell LJ argon film are shown in Figs. 5(a)-5(c) from $[0,0,\kappa_{z}]$ to $[\pi/a,0,\kappa_{z}]$ at $\kappa_{z}^{*}$ values of 0, 0.05, and 0.30. With the exception of $\kappa_{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-plane wave vectors. There are some deviations in the acoustic branches for $\kappa_{z}^{*}=0$ and $\kappa_{z}^{*}=0.05$. For example: (i) A 3D branch traced with small pieces of different 2D branches. This piece-wise mapping happens because the shapes of the 2D branches differ from those of the 3D branches, such that multiple branches are required to map across the BZ. (ii) A mapped 2D branch has frequencies lower than any 3D branch.

The second discrepancy can be examined by considering the mapped dispersion along $\kappa_{x}$ for a two unit cell film for $\kappa_{y}=\kappa_{z}=0$, as shown in Fig. 6. The downward frequency shift is present in all branches and can be attributed to two phenomena. First, there are additional branches (e.g., the optical branch with a $\boldsymbol{\Gamma}$-point frequency near $6$ Trad/s). Such branches are a result of the film anisotropy. As the film thickness decreases, the degenerate bulk transverse 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 at the top and bottom of the film, the restoring forces in the in-plane and cross-plane directions will be different, which leads to the branch splitting. This effect will occur for all degenerate acoustic and optical branches. Second, one of the acoustic branches is quadratic, which is expected in a two-dimensional system. A similar behavior will be present in the ten unit cell film and leads to a reduction in frequency.

The success of the mapping algorithm is also compromised by the existence of surface phonon 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 participate in a given mode and is defined as Feldman et al. (1993)

[TABLE]

The participation ratios for the 10 unit cell LJ argon film at the $\boldsymbol{\Gamma}$-point are plotted in Fig. 7(a), where four modes with $P<0.1$ are observed (the two plotted points each correspond to a degenerate mode). Upon a closer examination of their polarization vectors, the atomic motions in these modes are mainly at the film surfaces, with little to no motion deep within the film. Due to the spatial localization of surface modes, the mapping algorithm cannot map them to any 3D mode. The application of the mapping algorithm to the surface mode with a frequency of 5.72 Trad/s from Fig. 7(a) (colored red) results in the $\overline{q}$ sweep plotted in Fig. 7(b), which shows no distinct peaks.

The presence of the above effects at low thicknesses indicates an inherent difference between the 2D phonon modes from their 3D counterparts in LJ argon, going further than a simple limitation on the possible cross-plane wave vectors. In such cases, the rigorous 2D calculation is necessary to obtain accurate phonon properties.

There is also a source of numerical error that is separate from the anisotropy and the surface modes. For mapping destinations with the smallest cross-plane wave vector component for a given thickness [i.e., $n=\pm 1$ in Eq. (7)], more mis-mapped modes are present compared to other values of $n$. This behavior is evident in Figs. 5(b) and 5(c), which show the mapped 2D dispersions for a ten unit cell film at $\kappa_{z}^{*}=0.05$ ($n=1$) and at $\kappa_{z}^{*}=0.30$ ($n=6$). The $\kappa_{z}^{*}=0.30$ dispersion has a larger number of 2D mapped modes tracing the 3D branches, but the $\kappa_{z}^{*}=0.05$ dispersion has some 3D branches with almost no 2D mapped modes. These missing modes are a result of the difficulty of fitting the largest possible wavelength into the thickness of the 2D system. There are thus mis-mapped modes at other $\kappa_{z}^{*}$ because the modes that should have mapped to $\kappa_{z}^{*}=0.05$ must end up somewhere. This effect is present for all film thicknesses at the smallest allowed wave vector, but the increasing number of mapping destinations as the film gets thicker means that the overall effect for the entire 2D dispersion is reduced.

III.2.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 dispersions for the two and five unit cell graphene films are plotted in Figs. 8(a) and 8(b). Only the low frequency branches are shown, as higher frequency branches do not vary greatly with respect to $\kappa_{z}^{*}$ and are difficult to differentiate. Also plotted is the bulk graphite $\boldsymbol{\Gamma}-$M dispersion at $\kappa_{z}^{*}=0$. The mapped 2D modes in both the two and five unit cell films closely trace the 3D branches. The transverse branches are split in both 2D and 3D due to the anisotropy. A small downward shift in frequency is only observed in the two unit cell film, consistent with the DOS [Fig. 2(b)]. The 2D branches tend to completely map the 3D branches, rather than in the piece-wise manner seen for LJ argon. This result indicates that the shapes of the 2D branches are very similar to those of the 3D branches, which can be attributed to the weak interactions between the graphene layers.

IV Summary

We applied lattice dynamics calculations to study the relationship between phonon modes in corresponding 2D and 3D systems. Specifically, we developed a mapping algorithm based on normal mode decomposition to link 2D film and 3D bulk dispersions (Sec. II.3). LJ argon and graphene systems were analyzed due to their respective isotropy and anisotropy in bulk.

As shown in Figs. 2(a) and 2(c), the DOS of LJ argon films converges to the bulk DOS as the thickness increases. At low thicknesses, however, the 2D and 3D DOS diverge, showing the effect of anisotropy in the films that is not present in bulk. In contrast, the graphene DOS 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-layer interactions.

Upon application of the mapping algorithm, differences between the 2D and 3D phonon modes were observed in LJ argon. At larger thicknesses, the differences are small and the 2D modes 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 shifts and branch splitting (Fig. 6). These effects can be attributed to the emergence of anisotropy and surface modes in the LJ argon films. The graphene systems, with their weak inter-layer interactions, do not show the same deviations. Instead, their 2D phonon modes are strongly matched to the 3D phonon modes at specific cross-plane wave-vectors, with only small downward frequency shifts at low frequencies [Figs. 8(a) and 8(b)].

The observed differences between mapped 2D modes and their 3D counterparts are relevant to predictions of film thermal conductivity. The majority of previous calculations have assumed the existence of 3D modes whose mean free paths are modified with a boundary scattering model Jain et al. (2013); Cuffe et al. (2015); Fu et al. (2017); Wang and Huang (2015). In some studies, the number of cross-plane wave vectors was limited by the film thickness Turney et al. (2010); Wang and Huang (2015). While a treatment based on 3D modes may be appropriate 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 surface modes. In such cases, an accurate thermal conductivity prediction requires use of the 2D slab 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 is published, it will be found at Link.