Valley filtering effect of phonons in graphene with a grain boundary

TL;DR

This paper demonstrates that phonons in graphene can be selectively filtered by valley using grain boundaries, opening new avenues for phonon control and valley phononics in material science.

Contribution

It introduces the concept of valley filtering for phonons in graphene, a novel idea extending valleytronics to phononic systems.

Findings

Perfect valley filtering observed at specific energies for flexural phonon modes.

Valley filtering effectiveness is linked to phonon valley pocket anisotropy.

Fano-like resonance can enhance valley filtering performance.

Abstract

Due to their possibility to encode information and realize low-energy-consumption quantum devices, control and manipulation of the valley degree of freedom have been widely studied in electronic systems. In contrast, the phononic counterpart--valley phononics--has been largely unexplored, despite the importance in both fundamental science and practical applications. In this work, we demonstrate that the control of "valleys" is also applicable for phonons in graphene by using a grain boundary. In particular, perfect valley filtering effect is observed at certain energy windows for flexural modes and found to be closely related to the anisotropy of phonon valley pockets. Moreover, valley filtering may be further improved using Fano-like resonance. Our findings reveal the possibility of valley phononics, paving the road towards purposeful phonon engineering and future valley phononics.

Equations

$$\psi ={\mathbf{G}}^r{\left({\mathbf{G}}_0^r\right)}^{-1}{\psi }_0,$$

$${\psi }_l={\mathbf{G}}_{n,0}^r{\left({\mathbf{G}}_{L;0,0}^{(0)r}\right)}^{-1}{\mathbf{u}}_{Ln}(+)$$

$${\psi }_{N+1}=\sum _m{\tau }_{nm}{\mathbf{u}}_{Rm}\left(+\right).$$

$$t_{mn}$$

$${\Xi }_{mn}=|t_{mn}{|}^2.$$

$$\Xi (\omega )=\sum _{m,n}{\Xi }_{mn}(\omega ).$$

$$\begin{array}[]{l}{e^{i\left[{{k_{x}}\left({{x_{2}}-{x_{1}}}\right)+{k_{y}}\left({{y_{2}}-{y_{1}}}\right)}\right]}}={u_{n;2\alpha}}/{u_{n;1\alpha}},\\ {e^{i\left[{{k_{x}}\left({{x_{3}}-{x_{1}}}\right)+{k_{y}}\left({{y_{3}}-{y_{1}}}\right)}\right]}}={u_{n;3\alpha}}/{u_{n;1\alpha}},\end{array}$$

$$\left(\begin{array}[]{l}{k_{x}}\\ {k_{y}}\end{array}\right)={\left({\begin{array}[]{*{20}{c}}{{x_{2}}-{x_{1}}}&{{y_{2}}-{y_{1}}}\\ {{x_{3}}-{x_{1}}}&{{y_{3}}-{y_{1}}}\end{array}}\right)^{-1}}\left(\begin{array}[]{l}-i\log\left({{u_{2\alpha}}/{u_{1\alpha}}}\right)\\ -i\log\left({{u_{3\alpha}}/{u_{1\alpha}}}\right)\end{array}\right).$$

$$P_v=({\Xi }^K-{\Xi }^{K^{\prime }})/({\Xi }^K+{\Xi }^{K^{\prime }}),$$

$$\displaystyle{\bm{D}}_{C}=\left(\begin{array}[]{cc}\omega_{c}^{2}&t\\ t&\omega_{\textrm{defect}}^{2}\end{array}\right),$$

$${\bm{\Sigma}}^{r}_{L/R}=-\frac{i}{2}\left(\begin{array}[]{cc}\gamma(\omega)&0\\ 0&0\end{array}\right).$$

$${\mathbf{G}}_C^r=[(\omega +i\eta {)}^2\mathbf{I}-{\mathbf{D}}_C-{\mathbf{\Sigma }}_L^r-{\mathbf{\Sigma }}_R^r{]}^{-1}.$$

$$\Xi (\omega )$$

$$=\frac{{\left[\gamma \left(\omega \right)\right]}^2}{{[({\omega }^2-{\omega }_c^2)-t^2/({\omega }^2-{\omega }_{\mathrm{defect}}^2)]}^2+{\left[\gamma \left(\omega \right)\right]}^2}.$$

$$\psi \left(\mathbf{r}\right)$$

$$={\psi }_0\left(\mathbf{r}\right)+\int G^r\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)V\left({\mathbf{r}}^{\prime }\right){\psi }_0\left({\mathbf{r}}^{\prime }\right)d{\mathbf{r}}^{\prime },$$

$$\left[{\left(\omega +i\eta \right)}^2\mathbf{I}-{\mathbf{D}}_0\right]{{\mathbf{G}}_0}^r$$

$$\left[{\left(\omega +i\eta \right)}^2\mathbf{I}-{\mathbf{D}}_0-{\mathbf{D}}^{\prime }\right]{\mathbf{G}}^r$$

$$\left({\omega }^2\mathbf{I}-{\mathbf{D}}_0\right){\psi }_0=0,$$

$$\left({\omega }^2\mathbf{I}-{\mathbf{D}}_0-{\mathbf{D}}^{\prime }\right)\psi =0.$$

$$\left[{\left(\omega +i\eta \right)}^2\mathbf{I}-{\mathbf{D}}_0\right]\left({\psi }_0+{\mathbf{G}}_0^r{\mathbf{D}}^{\prime }\psi \right)$$

$$=$$

$$\left[{\left(\omega +i\eta \right)}^2\mathbf{I}-{\mathbf{D}}_0-{\mathbf{D}}^{\prime }\right]\left({\psi }_0+{\mathbf{G}}^r{\mathbf{D}}^{\prime }{\psi }_0\right)$$

$$=$$

$${\psi |}_{{\mathbf{D}}^{\prime }=0}={\psi }_0$$

$$\psi$$

$$={\psi }_0+{\mathbf{G}}^r{\mathbf{D}}^{\prime }{\psi }_0.$$

$${\mathbf{G}}^r={\mathbf{G}}_0^r+{\mathbf{G}}^r{\mathbf{D}}^{\prime }{\mathbf{G}}_0^r,$$

$$\psi ={\mathbf{G}}^r{\left({\mathbf{G}}_0^r\right)}^{-1}{\psi }_0.$$

$${\mathbf{u}}_{Ln}(+),n=1,...,\mathrm{M},$$

$${\psi }_l={\mathbf{G}}_{n,0}^r{\left({\mathbf{G}}_{L;0,0}^{(0)r}\right)}^{-1}{\mathbf{u}}_{Ln}(+)$$

$${\psi }_{N+1}=\sum _m{\tau }_{nm}{\mathbf{u}}_{Rm}\left(+\right).$$

$$t_{mn}$$

$$=\sqrt{\frac{v_m{a}_L}{v_n{a}_R}}{\tilde{\mathbf{u}}}_{Rm}^{†}\left(+\right){\mathbf{G}}_{N+1,0}^r{\left({\mathbf{G}}_{L;0,0}^{(0)r}\right)}^{-1}{\mathbf{u}}_{Ln}(+),$$

$${\tilde{\mathbf{u}}}_{Rm}^{†}\left(+\right){\mathbf{u}}_{Rn}\left(+\right)={\delta }_{mn}.$$

Full text

I Introduction

Owing to potential realization of low-dissipation devices, valleytronics has received considerable research interest. For electrons, degenerate but inequivalent valley states around the Fermi level form another degree of freedom besides the spin degree of freedom, providing a feasible way to design information-storage or -processing devices with low energy consumption. As the key factor for encoding information, valley degeneracy of electrons can be controlled by strain,Gunawan et al. (2006) magnetic fields,Zhu et al. (2012) optical pumping with circularly polarized light,Mak et al. (2012) substrates,Qi et al. (2015) $etc$. Also, valley-polarized current can be achieved by cyclic strain,Jiang et al. (2013) second-order nonlinear response,Yu et al. (2014) temperature gradient,Chen et al. (2015); Yu et al. (2015, 2016) and valley filters.Rycerz et al. (2007); Santos et al. (2009); Gunlycke and White (2011); Xu et al. (2016); Zheng et al. (2017)

For valleys of phonons, however, it remains largely unexplored. Particularly, phonons cannot be manipulated directly by electric fields or magnetic fields due to the lack of spin and charge degrees of freedom. To lift the degeneracy of phonon valleys, one needs to break the time-reversal symmetry and inversion symmetry, which is rather complicated.Liu et al. (2017a) Therefore, obtaining phonon valley polarization via valley filtering is an attractive and promising option.

For geometry-related valley filtering effects of electrons,Rycerz et al. (2007); Santos et al. (2009) the key ingredient is to make a single-valley region to filter out one valley of electrons. In these cases, no symmetry is explicitly required. For phonons, geometry-related valley filtering should be applicable similarly.

One particular type of valley filters, as learned from electrons, is grain boundaries (GBs), which are ubiquitous in chemical vapor deposition samples. GBs represent a class of interfaces exhibiting exotic electronic properties such as quasi-1D metallic statesLahiri et al. (2010); Zou and Yakobson (2015a, b) and selective valley filtering for electrons.Yazyev and Louie (2010); Gunlycke and White (2011) The influence of GBs on thermal conduction of phonons has been widely studied.Schelling et al. (2004); Kimmer et al. (2007); Lu and Guo (2012); Serov et al. (2013) However, valley filtering effect of phonons in GB systems lacks of investigations.

Motivated by recent development of both valleytronics and phononics,Schaibley et al. (2016); Li et al. (2004); Liu et al. (2017b) in this work, we explore the valley degree of freedom in an entirely different type of quasi-particles, $i.e.$, phonons. We show that phonon valley states can be filtered by using a grain boundary. In addition to valley-selective transport, Fano-like resonance appears in the transmission spectrum, enhancing the valley polarization of phonons further.

II Model and Methods

We illustrate our idea by investigating the transport properties of phonons through a 558 GB in graphene. The 558 GB is an extended line defect [Fig. 1(a)] containing carbon-atom pentagons and octagons.Serov et al. (2013); Huang et al. (2013) It is actually a 0∘ GB system,Phillips and Mele (2015) where selective valley filtering of electrons has been shown previously.Yazyev and Louie (2010); Gunlycke and White (2011) Graphene with a 558 GB was experimentally realizedLahiri et al. (2010) in 2010. Figure 1(a) shows the transport system, which consists of two semi-infinite thermal leads ($i.e.$, semi-infinite ideal graphene), and a central scattering region with a 558 GB. The structure is periodic along the transverse ($y$) direction.

In this work, we focus on the quasi-ballistic quantum transport regime, neglecting phonon-phonon and electron-phonon interactions. In realistic experimental systems, the graphene layer might be encapsulated with hBN layers or be placed on some substrates. In this work we shall not consider extrinsic effects due to phonon coupling between the different materials, i.e., the interlayer coupling of phonons is ignored such that the graphene layer behaves as a free-standing one.

II.1 Mode-Matching Method

To investigate transport behaviors of individual phonon modes of graphene, we utilize the mode-matching (MM) method,Ando (1991); Khomyakov et al. (2005); Rungger and Sanvito (2008) which has been widely used for investigating electronic transportZhang et al. (2014); Chen et al. (2015) and was recently applied to phonon transport.Ong and Zhang (2015); Sadasivam et al. (2017) In this part, we give a brief introduction of the mode-matching method. (More details can be found in the Appendix A.)

To apply the MM method, the transport model is artificially divided into“principal layers” along the transport direction and be labeled with layer number $\cdots-2,-1,0$,$1,2,\cdots,N$,$N+1,N+2,\cdots$ Here, $-\infty$ to [math] label the left thermal lead, $1$ to $N$ label the central scattering region, and $N+1$ to $+\infty$ label the right thermal lead. Principal layers are usually composed of several unit cells to ensure that only nearest-neighbor interaction exists between them. For one mode $\psi_{0}$ going through the device from the left thermal lead, it becomes $\psi$ at the right thermal lead. This evolution can be expressed in terms of phonon Green’s functions as:

[TABLE]

where ${{\bf G}^{r}}$ and ${{\bf G}_{0}^{r}}$ are the retarded Green’s of the central region with/without coupling between the thermal leads and the central region. More precisely, for a right-propagating eigen mode ${{\bf{u}}_{Ln}}(+)$ of the left thermal lead, it evolves into

[TABLE]

at layer $l$. To find out the transmission function, one may look at the wave function at layer $N+1$, which is the left-surface of the right thermal lead. Decomposing the wave function into eigen-modes ${{\bf{u}}_{Rm}}\left(+\right)$ of the right thermal lead, we have

[TABLE]

And the transmission amplitude of mode $n$ from the left thermal lead to mode $m$ of the right thermal lead is

[TABLE]

where $v_{Rm(Ln)}$ is the group velocity of mode $m(n)$ in R(L), $a_{L/R}$ is the principal-layer lattice constant of the $L/R$ thermal leads. Consequently, transmission of mode $n$ (in L) to mode $m$ (in R) is

[TABLE]

The total transmission function is obtained by adding up transmission coefficients of right-propagating modes:

[TABLE]

II.2 Computational Details

Force constants of the transport system were obtained using the finite-displacement methodParlinski et al. (1997); Togo et al. (2010) as implemented in the first-principles Realspace Electronic Structure CalcUlator (RESCU).Michaud-Rioux et al. (2016) Norm-conserving pseudo-potential for carbon atoms was used. Self-consistent calculations were performed within the local-density approximation (LDA) with energy criteria of $1\times 10^{-6}$ Hartree and charge variation criteria of $1\times 10^{-5}$.

For the calculation of phonon dispersion of graphene, the optimized carbon-carbon bond length of 1.40787 Å was used. In addition, a 4$\times$4$\times$1 supercell structure and a real-space resolution of 0.4 Bohr were employed.

For the calculation of the whole transport system, a $k$-mesh of 1$\times$3$\times$1 and a real-space resolution of 0.4 Bohr were chosen. The central region was relaxed until forces exerted on the central region were less than $1\times 10^{-3}$ Hartree/Bohr. Then, a $1\times 3\times 1$ supercell of the transport system was constructed with periodic boundary condition for the calculation of force constants. The constructed system contains 294 atoms and requires calculations of forces for 392 displaced configurations.

Transmission calculation based on the force constants was done using the PHonon MODEs (PHMODE) code. We adopted force constants of ideal graphene structure for the thermal leads to ensure that leads were in ideal crystalline structures. Then, for the utilization of the MM method, the transport system was divided into principal layers along the transport direction ($x$)Khomyakov et al. (2005) and a cutoff of 5 Å was applied to the force constants of the whole transport system to ensure that only nearest-neighbor interaction exists between principal layers.

II.3 Unfolding phonon band structures

By applying the mode-matching method, we can obtain mode-resolved transmission coefficients. However, the principal layers in leads are not primitive unit cells. What we solve from the mode-matching method are eigenvectors in the supercell and wave vectors $\bf K$ in the corresponding folded Brillouin Zone (BZ). To fully explore the mode dependence of transmission, it is vital for us to unfold the $\bf K$ to those $\bf k$-vectors within the $1^{\textrm{st}}$ BZ of the primitive cell.Huang et al. (2014); Zheng and Zhang (2017) Here, we propose a rather convenient way to realize this idea.

According to Bloch’s theorem, it is well known that wave function of a periodic system $\psi_{n{\bf k}}({\bf r})$ satisfies the relation: $\psi_{n{\bf k}}({\bf r+R})=\psi_{n{\bf k}}({\bf r})\exp(i{\bf k}\cdot{\bf R})$ for all $\bf R$ in the Bravais lattice, where $\bf k$ is the wave vector. As shown in Fig. 1(a) of the main text, we can pick up three equivalent atoms, for example, those labeled with 1, 2, 3, in a principal layer of the lead L. Then, we have

[TABLE]

where $x_{i},y_{i}(i=1,2,3)$ denote positions of the chosen atoms, $u_{n;i\alpha}$($\alpha=x$ or $y$ for in-plane modes, and $\alpha=z$ for out-of-plane modes) is the component of ${\bf u}_{n}$ on atom $i$ along the $\alpha$ direction. From the above equations, we can get the primitive wave vector as:

[TABLE]

Therefore, by combining the mode-matching method with Eq. (7), we can acquire the information of both wave vectors and transmission for an injection mode.

III Valley filtering effect

For completeness, phonon dispersion of graphene is shown in Fig. 1(b). In the left panel of Fig. 1(c), it is demonstrated that BZs are folded in our transport calculations. The unit cell used for the thermal lead L contains 4 primitive cells of graphene, and thus the corresponding BZ is folded. For example, a $\bf k$-vector in the primitive 1${}^{\textrm{st}}$ BZ ${\bf k}=(k_{x},k_{y})=(K_{x}+N_{1}{B}_{1},K_{y}+N_{2}{B}_{2})$ ($N_{1,2}=0,1,\cdots$; ${\bf B}_{1,2}$ are reciprocal lattice vectors of the unit cell used for thermal lead L; $|{\bf K}_{x(y)}|<|{\bf B}_{1(2)}|$) is reduced to be ${\bf K}=(K_{x},K_{y})$ in the folded BZ.

Figure 1(d,e) demonstrates mode transmission in both the folded BZ and primitive BZ at 200 and 400 cm*-1*. In the left panel of Fig. 1(d), the inner two circles in the folded BZ can be attributed to the TA and LA branches directly. Unfolding $(K_{x},K_{y})$ to primitive $(k_{x},k_{y})$ in the middle panel, one can observe that the other arcs actually originate from the ZA branch of graphene by comparison with the corresponding contour plot of graphene in the right panel.

For phonons at 200 cm*-1*, the TA and LA phonons almost fully transmit through the 558 GB; the ZA phonons have the largest transmission at perpendicular injection ($k_{y}=0$), and smaller transmission at a skew injection angle, which shows a similar trend as previous results in silicon with a high-energy $\Sigma$29 twist GB.Kimmer et al. (2007) The advantage of unfolding is more clearly seen at 400 cm*-1*. In this case, transmission of the LA and ZA branches is fairly high; by contrast, the TA branch has high transmission when the angle between $k_{y}$ and $k_{x}$ is roughly within 30 degree.

By examining the primitive modes, we are able to explore the valley-selective transmission in the GB system. Illustrated in Fig. 2(a), phonon “valleys” also exist in the ZA, ZO, TA, and LO branches, which emerge around the Dirac-like points [550 cm*-1* in Fig. 1(b)] and band extremum (1050 and 1227 cm*-1*) at $K$($K^{\prime}$). Similarly to the electron case, there are two inequivalent $K$ points, labeled as $K$ and $K^{\prime}$, at the corners of the $1^{\textrm{st}}$ BZ. Without losing generality, we choose $\omega=$500, 600, 800, 1300 cm*-1* to illustrate the valley-selective transmission spectrum in Fig. 2(b).

For phonons at 500 cm*-1*, TA modes have high transmission when the angle between wave vector $\bf k$ and the transport direction ($x$) is within 30 degree. The ZA valleys around the BZ corners, however, have different features. First, the shape of ZA valleys is triangular instead of circular or hexagonal. Second, transmission of ZA valleys at six corners is classified into two groups: three equivalent ones at $K$ valleys, and the other three equivalent ones at $K^{\prime}$ valleys. Third, high-transmission modes lie on the horizontal edge of the triangular valleys.

For phonons at 600 cm*-1*, TA phonons have much lower transmission. The contour shapes of the ZO valleys resemble those of ZA phonon valleys at 500 cm*-1*. However, high transmission modes lie in two different edges in ZO valleys from those edges in ZA valleys at 500 cm*-1*. Transmission of TA valleys at 800 cm*-1* has similar features to ZA valleys at 500 cm*-1*, and the LO valleys at 1300 cm*-1* to ZO at 600 cm*-1*.

To get a better idea of the valley filtering effect, we define $\theta=\arctan(v_{y}/v_{x})$ to be the injection angle, which ranges from $-90$ degree to $90$ degree for right-going modes. The angular dependence of valley transmission for those shown in Fig. 2(b) is depicted in Fig. 3(a). Angular dependence of valley transmission shows good symmetry, $\Xi^{K}(\theta)=\Xi^{K^{\prime}}(-\theta)$, at $\omega=500$, $600$ cm*-1*, and weak symmetry for $800$ and $1300$ cm*-1*. This deviation may be caused by the slight mirror-symmetry breaking of the optimized structure (the average deviation per atom is about 0.05 Å). The angles for achieving highest/lowest transmission for two valleys are roughly opposite to each other, $\theta_{m}(K)\approx-\theta_{m}(K^{\prime})$, which further leads to the valley filtering effect of the 558 GB system.

Defining valley polarization of transmitted modes as

[TABLE]

we obtain valley polarization at different $\omega$. Figure 3(b) depicts the valley polarization correspondingly to Fig. 3(a). All branches show significant valley polarization. Particularly, the ZA and ZO branches have highest valley polarization ($\pm$100% and $\pm$86%) at $|\theta|\to 90^{\circ}$, which is nearly parallel to the GB. Due to the triangular shape of valley pockets [Fig. 2(b)], transmission of right-going modes mostly occurs among $-90^{\circ}<\theta<30^{\circ}$ for one valley and $-30^{\circ}<\theta<90^{\circ}$ for the other valley. [Fig. 3(c)] As a consequence, valley polarization for ZA and ZO branches when $|\theta|>30^{\circ}$ is fairly high because one valley dominates transmission when $|\theta|>30^{\circ}$. Therefore, the filtering effect is closely related to the anisotropy of phonon valley pockets, implying that such a valley filtering effect is not limited in the 558 GB.

Also shown in Fig. 3(c), K and K’ valley pockets are nearly inverse symmetric. For a structure which has the mirror symmetry about the $x$ axis (the transport direction), the angular dependence of valley transmission functions around the K and K’ valleys should be $\Xi^{K}(\theta,\omega)=\Xi^{K^{\prime}}(-\theta,\omega)$ or $\Xi^{K}(\delta k_{x},\delta k_{y},\omega)=\Xi^{K^{\prime}}(\delta k_{x}^{\prime},-\delta k_{y}^{\prime},\omega)$, where $\delta{\bf{k}}={\bf{k}}-{\bf{K}},\delta{\bf{k^{\prime}}}={\bf{k^{\prime}}}-{\bf{K^{\prime}}}$ are measured from the bottom/top of the valley pockets. Therefore, when the transmission function of one phonon valley lacks the symmetry about $\theta=0$, we obtain $\Xi^{K}(\theta,\omega)\neq\Xi^{K^{\prime}}(\theta,\omega)$, which is valley-polarized, as illustrated in Fig. 4.

To further find out the filtering energy window for ZA and ZO modes, we plot the valley polarization for $\theta=30^{\circ}$ and $40^{\circ}$ in Fig. 3(d). From this plot, two perfect energy windows are found, [475, 510] and [555, 580], within which valley polarization of flexural modes (ZA/ZO) exceeds 95% when the injection angle reaches $40^{\circ}$.

IV Fano-like resonance

Besides transmission peaks, it is worth noting that there are transmission dips at $\omega=800$ and $1300$ cm*-1* in Fig. 3(a). The asymmetric line-shape around the transmission dip and peak implies a Fano-like resonance. To understand the physical origin of the Fano-like resonance, we illustrate the scattering wave functions of highly-transmitted modes [indicated by black arrows in Fig. 3(a)] in Fig. 5(a). In these plots, the atoms having the highest oscillation amplitude all locate at the GB region. However, the scattering wave function at $\omega=600$ cm*-1* is less local, which may explain the absence of dips. Based on the observation, we anticipate that the resonance is caused by the interplay between a continuum state from the left thermal lead, $i.e.$, graphene, and a localized state at the grain boundary. Localized quasi-1D metallic statesLahiri et al. (2010); Zou and Yakobson (2015a, b) for electrons at grain boundaries have been found previously. Therefore, it is a reasonable conjecture that local phonon states may also form at the GB.

The transmission as a function of phonon energy for a fixed $K_{y}$ mode is further demonstrated in Fig. 5(b). The asymmetry line-shape maintains in the energy domain. Previous theoretical investigations have demonstrated Fano resonances due to the coupling of optical phonon modes,Lee et al. (2006); Li et al. (2012) bound excitons to LO phonons,Xu et al. (2005) and coupling of phonon modes in a 1D nonlinear chain model.Kim and Kim (2001) Inspired by the role of local states in Fano resonance, we build a simple model to interpret the GB system: as shown in Fig. 5(c), there is one local vibration mode $\omega_{\textrm{defect}}$ that interacts with a conducting mode $\omega_{c}$, which couples directly to two thermal leads.Kim et al. (2012) The dynamical matrix of the central region without thermal baths can be written as

[TABLE]

where $t$ represents the interaction of these two modes. The influence of two thermal baths can be counted in using self-energy matrices as

[TABLE]

The corresponding bandwidth functions are ${\bm{\Gamma}}_{L/R}=-2\textrm{Im}({\bm{\Sigma}}^{r}_{L/R})$. The retarded Green’s function of the central region is

[TABLE]

And the transmission spectrum can be calculated asChen et al. (2018)

[TABLE]

Because ${\bf\Sigma}^{r}(-\omega)={\bf\Sigma}^{r*}(\omega)$,Wang et al. (2007); Xu et al. (2008) $\gamma(\omega)$ should be an odd function of $\omega$. Then, the simplest expression of $\gamma(\omega)$ is $\gamma(\omega)=\omega\gamma_{0}$, which is actually the wide-band limit model for phonons. The expression of $\gamma(\omega)$ is not important, because odd functions such as $\omega^{3}$, $1/\omega$, $\sin\omega$, and $\omega/(\omega^{2}+1)$ also lead to similar Fano-like line-shapes.

Based on Eq. (15), Fig. 5(d) shows that when the interaction $t$ is rather weak compared to $\gamma_{0}^{2}$, Fano-like resonance manifests around the local mode $\omega_{\textrm{defect}}$. The line-shape shown is partially consistent with the one shown in Fig. 5(b). This simple model demonstrates the existence of resonance and anti-resonance points of transmission as a function of energy. However, it does not explicitly contain angular dependence, which may be further considered by introducing $\bf{k}$-dependence in the model. Our calculations enrich the physics by showing that Fano-like resonance from pure phonon modes also occur in 2D materials. Since dips and peaks emerge at such resonance points, valley polarization can be significantly enhanced around the resonant points.

V Conclusion

In summary, we investigated the phonon transport properties in a 2D graphene sheet which contains an extended zero-angle grain boundary (558 GB). Combining real-space first-principles calculations with the extended mode-matching method to recover the primitive wave vectors of injection modes, we revealed that the 558 GB selectively scatters different phonon modes, leading to valley-selective transport. The valley filtering effect is closely related to the anisotropy of valley pockets of phonons. In particular, the 558 GB can be used as a perfect phonon valley filter for ZA and ZO valleys within [475, 510] and [555, 580] cm*-1*, respectively. In addition, Fano-like resonance was revealed, which can further enhance the valley polarization of phonons.

Due to crucial computational load, only one type of grain boundary was investigated. However, since the valley filtering effect is closely related to the anisotropy of valley pockets, we anticipate that phonon valley filtering effect can be generally observed in various GB systems and expect that our method will be beneficial in revealing novel properties of phonons in emerging 2D materials.Zou and Yakobson (2015b); Xu et al. (2013); Zou et al. (2013)

Acknowledgements.

The authors want to thank Prof. Xiaolong Zou, Prof. Jian-Sheng Wang, Prof. Xiaohong Zheng, and Dr. Lei Zhang for helpful discussions. We gratefully acknowledge financial support by NSF-China [Grant Nos. 11704257 (X.C.), 11874035 (Y.X.), and 11704238 (J.W.)], the General Research Fund (Grant No. 17311116), the University Grant Council (Contract No. AoE/P-04/08) of the Government of HKSAR (J.W.), and the NSERC of Canada (H.G.). We thank Calcul Québec and Compute Canada for the computation facilities.

Appendix A Mode-matching method

In this section, we briefly introduce the basic ideas of modified mode-matching method for studying quasi-ballistic transport properties of phonons.

A.1 Lippmann-Schwinger Equation

Suppose we have an incoming wave function $\psi(\bf r)$ going through the scattering region. To get the wave function across the scattering region, we can use either quantum mechanics methods plus some connection conditions to solve the problem, or the Lippmann-Schwinger (LS) equation. The LS equation states that the scattering wave function $\psi$ is composed of the original incoming wave function ${\psi_{0}}$ and the scattered wave function as:Wang and Guo (2009)

[TABLE]

where $V(\bf r^{\prime})$ is the scattering potential, $G^{r}\left({{\bf{r}},{\bf{r}}^{\prime}}\right)$ and ${G_{0}^{r}}\left({{\bf{r}},{\bf{r}}^{\prime}}\right)$ are the propagators of the mode with and without scattering from the central region. Using discrete lattice model, the retarded phonon Green’s functions for the system without and with scattering ($G_{0}^{r}$ and ${G^{r}}$) can be defined as

[TABLE]

where ${\bf D}_{0}$ is the dynamical matrix of the central part and ${\bf D}^{\prime}$ is the perturbed dynamical matrix inducing scattering. If ${\psi_{0}}$ is an eigen mode with the energy $E=\hbar\omega$ ($\hbar=1$ hereafter) from the left thermal lead, $\psi$ should also be an eigen state of the two-probe system with the same energy due to the law of energy conservation in quasi-ballistic regime, i.e.,

[TABLE]

Combining Eqs. (18-21), we have

[TABLE]

Therefore, under boundary conditions that

[TABLE]

we can express the scattering wave function $\psi$ in terms of phonon Green’s functions as

[TABLE]

A.2 Transmission

With the information of scattering states, it is straight forward to calculate the transmission function. Combining Eq. (28) and the Dyson Equation,

[TABLE]

we can reformulate $\psi$ in terms of phonon Green’s functions as

[TABLE]

Further by decomposing $\psi$ into eigen channels of the right thermal lead, we can obtain mode-resolved transmission. In detail, the right-going modes from the left thermal lead are denoted as

[TABLE]

where $M$ represents the degree of freedom in a principal layer. Going through the central scattering region, ${{\bf{u}}_{Ln}}(+)$ becomes

[TABLE]

at the $l^{\textrm{th}}$ principal layer. Here, the $0^{\textrm{th}}$ layer is the surface principal layer of the left thermal lead and ${{\bf G}^{(0),r}_{L;0,0}}$ is the retarded Green’s function at the $0^{\textrm{th}}$ principal layer of the semi-infinite left thermal lead, which hosts eigenmodes $\{{\bf{u}}_{Ln}(+)\}_{n=1}^{M}$. Decomposing the scattering wave function at the surface layer of the right thermal lead, $\psi_{N+1}$, into right-going eigenmodes of the right thermal lead, we have

[TABLE]

Then, we get elements of the transmission amplitude matrix as

[TABLE]

where $v_{m(n)}$ is the group velocity of mode $m(n)$, $a_{L/R}$ is the principal-layer lattice constant of the $L/R$ thermal leads, and dual functions $\tilde{\bf{u}}_{Rm}\left(+\right)$ satisfy

[TABLE]

From Eq. (LABEL:eq:tmn), it is clear that transmission coefficients of phonon modes going from the left to the right thermal lead can be obtained if we have the information of eigen states of both thermal leads, phonon Green’s function of the central scattering region, and surface phonon Green’s function of the left thermal lead.

A.3 Obtaining incoming modes at a given energy

In this section, we shall discuss how to get eigen phonon modes at a given energy for a periodic system. Eigen modes of the system can be solved according to the following equation as:

[TABLE]

where ${\bf D}$ is the dynamical matrix of the system. This equation leads to a complete set of orthogonal eigenstates with the eigenvalue $\omega_{0}^{2}$. We may write the wave function in real space in terms of local wave functions at each principal layers (or unit cells) as

[TABLE]

Substituting the above equation into Eq. (42), we have a chain-like relationship between adjacent principal layersAndo (1991):

[TABLE]

Because this is a periodic system, the Bloch’s theorem guarantees that wave functions at adjacent principal layers differ by only a phase factor:

[TABLE]

with $\lambda_{n}={e^{iK_{n}a}}$. ($a$ is the lattice constant, which is $A_{1}$ in our calculation.) Replacing ${\bf{c}}_{i\pm 1}$ in Eq. (43) by ${\bf{c}}_{i}$ according to Eq. (44) and using an eigenmode ${\bf u}_{n}$ instead, we get

[TABLE]

which is a quadratic equation for $\lambda_{n}$. To facilitate the solution, we can rewrite the above equation into a linear one:

[TABLE]

or in the form of matrix:

[TABLE]

Note that the matrix form is not unique. Then, we can solve the equivalent equation

[TABLE]

where

[TABLE]

and get 2$M$ eigenmodes and the corresponding eigenvalues. Among the 2$M$ eigenmodes, there are $M$ right-going and $M$ left-going modes. And both of them may contain propagating and evanescent modes, judging by the value of $\lambda_{n}$:Khomyakov et al. (2005) (1) if $|\lambda_{n}|<1$ or $>1$, ${{\bf{u}}_{n}}$ is an evanescent mode; (2) if $|\lambda_{n}|=1$, ${{\bf{u}}_{n}}$ is a Bloch propagating mode.

A.4 Velocity of incoming modes

Differentiating Eq. (45) over $K_{n}$, we have

[TABLE]

Therefore, the group velocity of a normalized propagating eigenmode is ($\hbar=1$)

[TABLE]

There is an extra factor of $1/\omega_{0}$ compared to the electron case.Khomyakov et al. (2005); Rungger and Sanvito (2008) Considering that $\omega_{0}$ is a given constant when obtaining mode-resolved transmission spectrum in Eq. (LABEL:eq:tmn), we may use

[TABLE]

instead of $v_{n}$ in Eq. (LABEL:eq:tmn) for calculating transmission.

A.5 Effective Dynamical Matrix

To get mode-resolved transmission, we constructed the effective dynamical matrix at each $K_{y}$ point as:Ong and Zhang (2015); Sadasivam et al. (2017)

[TABLE]

where $A_{y}$ is the lattice constant of the transport system along the $y$-direction, and $n_{x}$($n_{y}$) marks the position of (adjacent) unit cells. In our model, $A_{y}=A_{2}\approx 4.88$ Å . Using the effective 1D dynamical matrix, we obtained $K_{x}$-resolved transmission at each $K_{y}$ point based on the mode-matching method.Wang and Wang (2009); Ong and Zhang (2015); Chen et al. (2015)