Frequency splitting of chiral phonons from broken time reversal symmetry in CrI₃

The atomic vibrations that are present in molecules and solids at zero and finite temperature play a crucial role in their thermodynamic and transport properties. First-principles calculations based on density-functional theory (DFT) have been established as a powerful tool for understanding and predicting lattice-dynamical properties, including phonon dispersion [1, 2] and electron-phonon coupling [3, 4]. The key quantity underlying the calculation of lattice dynamics is the interatomic force constant (IFC) matrix, which is constructed by finding derivatives of the nuclear forces with respect to nuclear positions, either directly via finite displacements or though density functional perturbation theory [1, 2].

In presence of magnetic ordering, the change in the electronic ground state compared to the nonmagnetic case propagates to the IFCs [5, 6, 7, 8, 9]. However, since it is defined and calculated as a static response function, the IFC matrix is invariant under time reversal by construction. Thus, a description of the nuclear dynamics based solely on the IFCs will not correctly reflect the vibration mode degeneracies in a magnetic system; instead, the phonon frequency spectrum will be determined by the nonmagnetic symmetry group.

There has been significant recent work on the explicit inclusion of time-reversal symmetry (TRS) breaking in the nuclear equations of motion via the velocity dependence of the interatomic forces, applied to models [10, 11, 12, 13] and magnetic molecules [14]. This “velocity-force” coupling can be obtained from the nuclear Berry curvature, which describes the evolution of the phase of the electronic wavefunction with changes in nuclear coordinates [10, 11, 14]. A key result of including this coupling is that degenerate vibrational modes may split into nondegenerate chiral modes [15] with a well-defined finite angular momentum [16, 14], even at the Brillouin-zone center. Thus, the correct treatment of magnetic symmetry is crucial for elucidating the role of atomic vibrations in thermal Hall [17, 11, 18, 19, 20] and other effects involving TRS-broken lattice dynamics [21, 22, 23, 24, 25, 26, 27, 28].

A key assumption underlying the velocity-force approach in previous works [10, 11, 12, 13, 14] is that the time scale for electronic dynamics is fast compared to nuclear dynamics. However, this may completely break down in some systems, e.g., when the nuclear Berry curvature results from nuclei coupling to spins, whose dynamics are not necessarily faster than the atomic vibrations. The breakdown of the adiabatic picture could have a profound effect on the predicted splitting of chiral phonon modes.

In this work, we illustrate such a situation using the bulk-layered magnetic insulator CrI₃ as an example system. We first develop a DFT methodology amenable to both molecules and solids for computing phonons in the presence of velocity-force coupling. We apply this method to calculate the zone-center phonons in CrI₃, demonstrating the splitting of otherwise degenerate phonons into chiral modes. Next, we show that the velocity-force response is dominated by the canting of the spins on the Cr sites caused by atomic displacements. The dynamics of such spin canting are characterized by the zone-center magnon frequencies. The fact that magnon frequencies in CrI₃ are on the same order or smaller than the optical phonon frequencies [29, 30, 31, 32, 33, 34, 35, 36] means that we must treat spins and atomic displacements on the same footing. We develop a minimal model of this kind, and show that the adiabatic velocity-force approach can greatly overestimate the frequency splitting of the chiral modes.

We begin by reviewing the formalism of the velocity-force coupling from previous works [10, 11, 12, 13, 14], which we will refer to as the “Mead-Truhlar” (MT) approach. In order to simplify the discussion, we assume a finite system, generalizing to an infinite crystal below. The starting point of the derivation is the Born-Oppenheimer approximation, where the system wavefunction is factored into nuclear and electronic parts such that the ground-state electronic wavefunction $|\psi(\mathbf{R})\rangle$ depends parametrically on the nuclear coordinates $\mathbf{R}$ [37]. Once the electronic degrees of freedom are integrated out, the effective Hamiltonian for the nuclear wavefunction becomes [13]

where Roman indices (here $i$) run over nuclei, Greek indices (here $\alpha$) run over Cartesian directions, $p_{i\alpha}$ is the momentum operator for nucleus $i$ along direction $\alpha$, and $m_{i}$ is the mass of nucleus $i$. Using the notation $\partial_{i\alpha}=\partial/\partial R_{i\alpha}$, $A_{i\alpha}(\mathbf{R})=i\braket{\psi(\mathbf{R})\mid\partial_{i\alpha}\psi(\mathbf{R})}$ is a nuclear Berry potential, and $V_{\mathrm{eff}}(\mathbf{R})=\epsilon(\mathbf{R})+\sum_{i\alpha}\frac{\hbar^{2}}{2m_{i}}\bigl{(}\braket{\partial_{i\alpha}\psi(\mathbf{R})\mid\partial_{i\alpha}\psi(\mathbf{R})}-A_{i\alpha}(\mathbf{R})^{2}\bigr{)}$ is an effective scalar potential, where $\epsilon(\mathbf{R})$ is the ground-state energy for a given fixed nuclear configuration. As first pointed out by Mead and Truhlar [13], the nuclear Berry potential $A_{i\alpha}$ cannot always be made to vanish by changing the gauge of $\ket{\psi(\mathbf{R})}$ via the choice of an $\mathbf{R}$-dependent phase factor.

The resulting vibrational modes are then found by solving the equations of motion which, to harmonic order in nuclear displacements, are given by [10, 13, 11, 12]

Here $\mathbf{M}$ is a diagonal nuclear mass matrix $M_{i\alpha,j\beta}=m_{i}\delta_{i,j}\delta_{\alpha,\beta}$, $\omega_{n}$ is the frequency of mode $n$, and $\eta_{n}(\tau\alpha)$ is the component of the eigendisplacement of nucleus $\tau$ along direction $\alpha$ normalized so that $\eta_{n}^{\dagger}\mathbf{M}\eta_{m}=\delta_{nm}$. $\mathbf{G}$ is the “velocity-force matrix,” whose elements $G_{i\alpha,j\beta}$ relate the force on nucleus $i$ along direction $\alpha$ to the velocity of nucleus $j$ along direction $\beta$, and is expressed as

which is just $\hbar$ times the nuclear Berry curvature. The matrix $\mathbf{C}$ is the IFC matrix, which we define to be $C_{i\alpha,j\beta}=\partial_{i\alpha}\partial_{j\beta}\epsilon(\mathbf{R})$. Note that following Eq. (1) one could alternatively define $\mathbf{C}$ in terms of the Hessian of $V_{\text{eff}}$, which includes additional terms compared to the conventional IFC. However, the additional terms are higher order in the inverse nuclear mass and do not involve breaking of time-reversal symmetry, so that we neglect them in this work.

This formalism can be extended to the calculation of phonons in infinite crystals within DFT. We restrict ourselves to phonons at $\Gamma$, the Brillouin zone center, and work from a discrete set of DFT calculations on primitive-cell structures in which each sublattice displacement in turn is considered. The C matrix is constructed in the usual way from finite differences of the forces, while G is built from Berry phases computed over triangular configuration paths involving a given pair of sublattice-displaced structures and the undistorted structure. Berry phases are computed on a per-unit-cell basis for the Bloch manifold at each wave vector $\mathbf{k}$ and then averaged over the Brillouin zone to get the desired G matrix elements. The details are given in Sec. S1 of the supplemental material (SM) [38].

We perform calculations on CrI₃ in the ferromagnetic ground state using the vasp code [39, 40, 41], the local density approximation exchange-correlation functional [42], and projector-augmented wave potentials [43]. Semicore ($3s$, $3p$) and ($5s$, $5p$) states are included in the valence for Cr and I respectively. A $5\times 5\times 5$ Monkhorst-Pack grid [44] is used to sample the Brillouin zone, and the energy cutoff for the plane-wave basis set is 520 eV. Spin-orbit coupling, which is essential to the physics described here, is included in all calculations.

As mentioned above, the conventional calculation of phonons neglects G in Eq. (2), resulting in an equation of motion with TRS ($\tilde{\omega}_{n}^{2}\,\mathbf{M}\,\tilde{\eta}_{n}=\mathbf{C}\tilde{\eta}_{n}$). The resulting zone-center phonon modes with frequencies $\tilde{\omega}_{n}$ can be represented with real eigendisplacements $\tilde{\eta}_{n}$. In CrI₃ the representation of these modes at the zone center can be decomposed into the real irreducible representations (irreps) of $\bar{3}$ as $4A_{g}\oplus 4A_{u}\oplus 4E_{g}\oplus 4E_{u}$. The frequencies $\tilde{\omega}_{n}$ range up to 32 meV (see the SM [38] Table SIII). We can consider the velocity-force coupling between phonon modes obtained from the IFC alone by defining $\bar{G}_{nm}=\tilde{\eta_{n}}^{{\dagger}}\mathbf{G}\tilde{\eta}_{m}$; $\bar{\mathbf{G}}$ is block diagonal, with each block corresponding to a real irrep. Once G is included in Eq. (2), TRS is broken and the twofold-degenerate $E_{u}$ and $E_{g}$ modes split, as the corresponding irreps further break up into complex one-dimensional representations.

In Fig. 1(b), we show a selection of zone-center phonon frequencies of bulk ferromagnetic CrI₃ (see the SM [38] Table SIII for a complete list of frequencies). On the left side we plot $\tilde{\omega}_{n}$, i.e., the frequencies neglecting the velocity-force contribution, and on the right the frequencies including the velocity-force contribution via the MT approach (i.e., $\omega_{n}$). In Fig. 1(a) we plot the sum of the magnitudes of the velocity-force coupling terms relevant to each phonon eigenvector, as a measure of the strength of the coupling. We see that the $E_{g}$ modes have the strongest coupling. We shall see why shortly.

The splitting of each of the two-fold degeneracies shown in Fig. 1(b) results in chiral phonons with well-defined angular momentum in the $z$ direction (due to the symmetry of CrI₃). For mode $n$, the angular momentum is found via $l_{z}^{(n)}=2\hbar\sum_{\tau}m_{\tau}\mathrm{Im}[\eta_{n}^{*}(\tau x)\eta_{n}(\tau y)]$[21, 45], where $\tau$ runs over the atomic sublattices and $m_{\tau}$ is the mass of nucleus $\tau$. If we neglect G in Eq. (2), then the angular momentum of modes in a degenerate subspace spanned by ($\tilde{\eta}_{1},\tilde{\eta}_{2}$) depend on the basis we choose. Clearly, $l_{z}$ vanishes in the real basis $\tilde{\eta}_{n}$ introduced above, while a complex “circularly polarized” combination of degenerate modes of the form $\eta^{\prime}_{\pm}=\tilde{\eta}_{1}\pm i\tilde{\eta}_{2}$ will have equal and opposite $l_{z}$ (since we are at the $\Gamma$ point [16]), with the magnitude determined by the mode displacement patterns. Neglecting terms in the $\bar{\mathbf{G}}$ matrix that mix different degenerate subspaces, the chiral modes after splitting will have exactly the eigendisplacements $\eta^{\prime}_{\pm}$. Fig 1 (c) shows the angular momentum ($l_{z}$) for each mode after the inclusion of the G matrix. The fact that some of the split modes do not have equal and opposite $l_{z}$ indicates mixing between the degenerate subspaces in $\bar{\mathbf{G}}$. See Sec. S4 of the SM for a full analysis of the angular momentum of the chiral modes[38].

We now analyze the mechanisms responsible for the velocity-force coupling in CrI₃. In this material, the magnetic moments reside on the Cr atoms and are initially oriented out of plane along $z$. We note that the modes with large $\bar{\mathbf{G}}$ matrix elements in Fig. 1 are mostly those involving displacements of the I sublattices, which carry strong spin-orbit coupling. Under such displacements, it is natural that the Cr moments may cant, reflecting a local change in magnetic easy axis. This canting will result in a spin Berry curvature, and thus a contribution to G.

In fact, we find that for CrI₃, the spin Berry curvature is the dominant mechanism of velocity-force coupling. We demonstrate this by calculating the matrix elements of G under the assumption that only the spin-Berry-curvature mechanism is present. More specifically, we first approximate Eq. (3) as $\bar{G}_{mn}\simeq B_{Ia,m}G_{Ia,Jb}B_{Jb,n}$, where $G_{Ia,Jb}=-2\hbar\,\textrm{Im}\langle\partial_{Ia}\psi|\partial_{Jb}\psi\rangle$ is the Berry curvature of the wavefunctions in the parameter space spanned by the Cr spins; $B_{Ia,n}=\partial s_{Ia}/\partial\tilde{u}_{n}$ is a “spin canting matrix” describing the static change in the equilibrium spin unit vector on magnetic Cr site $I$ in direction $a$ resulting from phonon perturbation $n$; and $\tilde{u}_{n}$ is the amplitude of mode $n$ such that the set of atomic displacements are given by $\tilde{u}_{n}\tilde{\eta}_{n}$ ¹¹1Here $a$ runs only over $x,y$ directions corresponding to the two possible tilt angles. We have assumed the total magnitude of the spin is unchanged and that tilt angles are small so that $S_{Ix}=S\sin(\theta_{Ix})\approx Ss_{Ix}$. Under the assumption that the spin Berry curvature dominates, we can further approximate $G_{Ia,Jb}=-S\,\delta_{IJ}\epsilon_{ab}$, where spin $S=3\hbar/2$ for Cr and $\epsilon_{ab}$ is the antisymmetric tensor; then

A comparison between this spin-Berry approximation and the full Berry-phase calculation of the $\bf G$ matrix is presented in Table SII of the SM [38]. For some modes, the spin-Berry approximation accounts for only $\sim 60$% of the full velocity-force contribution; this means that other contributions, e.g., “phonon-only Berry curvature” (the Berry curvature from atomic displacements at fixed spin) are significant. However, these modes exhibit relatively small splittings from the full velocity-force contribution, so the errors are not significant in absolute terms. For the two $E_{g}$ modes that exhibit the largest splitting in the MT calculation [green, orange, and blue curves in Fig. 1(b)], the errors arising from this approximation are less than one percent.

This is a remarkable result, and is one of the main findings of the present work. By adopting the spin-Berry approximation, one can bypass the cumbersome computation of Berry phases for paths connecting distorted structures. Instead, all we need to take from the DFT calculations is the information about the spin canting in response to phonon distortions. In particular, it is now clear why the $\mathbf{G}$ tensor elements are so much smaller for the E_u modes; these are the ones that couple to the optical magnons, whose much larger stiffness strongly suppresses the spin canting.

A critical implication of the fact that the velocity-force coupling in CrI₃ results from spin canting is that the assumption underlying the MT approach of Eqs. (1-3), namely, that all electronic dynamics are fast compared to that of the phonons, is clearly unfounded. This is because the relevant time scale for spin dynamics is that of the magnon frequencies in the system. The experimentally measured zone-center magnons of CrI₃ have frequencies of 0.3 meV (i.e., the magnetocrystalline anisotropy) for the acoustic branch, and 17 meV for the optical branch [36], while the relevant phonons with the largest velocity-force coupling have frequencies in the range of 6-14 meV [see Fig. 1(b)].

Thus, an appropriate description of the low-energy dynamics must treat spins and phonons in this system on the same footing. To illustrate how this can be done, we focus on a single pair of $E_{g}$ or $E_{u}$ modes, which couple respectively either to an effective acoustic spin unit vector ${\bf s}=({\bf s}_{1}+{\bf s}_{2})/\sqrt{2}$ or its optical counterpart ${\bf s}=({\bf s}_{1}-{\bf s}_{2})/\sqrt{2}$. Denoting the phonon mode amplitudes and momenta as $(x,y)$ and $(p_{x},p_{y})$, the coupled spin-phonon Hamiltonian takes the form

Here $\tilde{\omega}$ is the bare phonon frequency, $\alpha=\partial^{2}E/\partial^{2}s_{x}$ is the spin anisotropy energy, and $\gamma=\partial^{2}E/\partial x\partial s_{x}$ is the coupling between the spin and the pair of phonons (which have been chosen such that the $xs_{y}$ and $ys_{x}$ terms vanish). Note that the unperturbed magnon frequency is related to the anisotropy by $\omega_{\rm m}=\alpha/S$, where $S=3\hbar/2$ is the Cr spin, and that the $\bf B$ matrix introduced above reduces in this minimal model to $B=\gamma/\alpha$. Going over to circularly polarized coordinates via $x_{\pm}=(x\pm iy)/\sqrt{2}$ and $s_{\pm}=(s_{x}\pm is_{y})/\sqrt{2}$, the equations of motion become

which are easily solved numerically.

In Table 1, we compare the frequency splittings of doubly-degenerate modes determined by the adiabatic Mead-Truhlar (MT) approach and spin-phonon (SP) model. The last column in Table 1 presents the results for the spin-phonon model; the $E_{g}$ modes couple to the acoustic magnon branch and the $E_{u}$ modes to the optical branch (we use experimental [36] magnon frequencies of 0.3 and 17 meV, respectively) ²²2A more general description of spin-phonon coupling, treating all relevant nuclear and spin degrees of freedom simultaneously, will be presented in a forthcoming communication.. A more detailed comparison of Mead-Truhlar approach and spin-phonon model is presented in Sec. S2 of the SM[38].

In order to understand the general features of the spin-phonon mixing embodied in Eq. (6), we plot in Fig. 2 the solutions of Eq. (6) as a function of phonon frequency $\tilde{\omega}$, with the magnon frequency set to $\omega_{\rm m}=17$ meV$/\hbar$ (i.e., the optical branch [36]), and $\gamma$ fixed at 2 meV^3/2/$\hbar$, a typical value for the $E_{u}$ modes in CrI₃ (since those are the phonons that couple to the optical magnon). The solid blue curve in the inset shows the splitting of the modes with dominant phonon character also as a function of $\tilde{\omega}$. Outside of the small “resonant” regime $\tilde{\omega}\simeq\omega_{\rm m}$, where significant magnon-phonon hybridization occurs, the mode splitting is well described by $SB^{2}F$ (see orange dashed curve in inset of Fig. 2), where $SB^{2}$ is the relevant force-velocity term in the adiabatic spin-Berry approximation [i.e., Eq. (4)] and $F=|1-(\tilde{\omega}/\omega_{\rm m})^{2}|^{-1}$. At small $\tilde{\omega}$, the magnon can be treated as a high-energy degree of freedom which renormalizes the phonons, and the splitting from the MT approach is recovered. Increasing $\tilde{\omega}$ toward $\omega_{\rm m}$ enhances the splitting of modes over the value at the adiabatic (MT) limit, peaking at the point where the magnon and phonon frequencies coincide and the modes have maximum hybridization. Above $\omega_{\rm m}$, the splitting decreases from its peak value, but is still enhanced over that determined by the MT approach for a range of $\tilde{\omega}$. At $\tilde{\omega}\gg\omega_{\rm m}$ the splitting is suppressed, and vanishes as $\tilde{\omega}\rightarrow\infty$.

This behavior is reflected in the splittings given in Table 1. Two of the three $E_{u}$ modes which couple to the optical magnon have frequencies below $\omega_{\rm m}=17$ meV$/\hbar$, but still in ranges where the splitting is enhanced above the adiabatic MT limit. The physical interpretation here is that the phonons are becoming hybridized with the magnon instead of just having their frequencies renormalized by it. However, since the splittings in the adiabatic limit (and thus the corresponding G matrix elements) are proportional to $\omega_{\rm m}^{-2}$, the splittings remain quite small. The largest frequency $E_{u}$ mode has $\omega_{\rm m}\ll\tilde{\omega}$, so that the splitting is reduced compared to values obtained from the velocity-force approach.

The $E_{g}$ modes couple to the acoustic magnon ($\omega_{\rm m}=0.3$ meV$/\hbar$), so for all $E_{g}$ modes $\tilde{\omega}\gg\omega_{\rm m}$, which is precisely the opposite of the limit in which the adiabatic velocity-force theory is applicable. This results in the drastic reduction ($F\ll 1$) in the splitting of the $E_{g}$ modes in Table 1 compared to the MT description in Eq. (2). The physical interpretation of this regime is that the Cr spins cannot keep up with the phonons, and thus the area swept out from the spin canting is greatly reduced compared to the assumption that they follow the nuclear motion adiabatically.

Clearly, these results have significant implications for experimental measurements of chiral phonons in CrI₃. Optical techniques, possibly utilizing circular polarization, constitute a powerful tool for studying such properties [48, 22, 49, 50, 25]. The strongly suppressed frequency splitting (SP column of Table 1) of the Raman-active $E_{g}$ modes is likely to be difficult to detect. This is consistent with recent work on CrBr₃, a related system, which found signatures of the chiral phonons but did not report a splitting of these modes [25]. The larger, though still quite modest, splitting of the $E_{u}$ modes could in principle be measured by peak shifts in infrared absorption [51], while direct detection of chirality would require circularly polarized infrared spectroscopy as in Refs. [48, 49].

More generally, our results also have several implications for finding other systems with a large splitting of chiral phonon modes at $\Gamma$. In materials like CrI₃, where the spin-Berry mechanism is responsible for the majority of the velocity-force coupling, it is most promising to look for (or engineer via, e.g., strain or magnetic field) cases where the relevant phonon and magnon frequencies coincide. This avoids the suppression of the splitting in the $\tilde{\omega}\gg\omega_{\rm m}$ regime, and the small spin canting likely when $\omega_{\rm m}\gg\tilde{\omega}$. Systems with lower-frequency optical magnons that maintain similar or larger spin-phonon coupling $\gamma$ are also strong candidates for observing larger effects. In any case, we can see from Table 1 that correctly accounting for the relative dynamics of spins versus phonons is necessary to avoid significantly overestimating the splitting of certain chiral modes.

In conclusion, we have developed and implemented a first-principles methodology for capturing time-reversal-symmetry-broken lattice dynamics in magnetic solids, and applied it to the case of bulk ferromagnetic CrI₃. We show that in this system, the previously-made assumption of fast electron dynamics compared to the lattice breaks down, since the relevant coupling is between Cr spins and atomic displacements. With a minimal model, we demonstrate that spins and phonons must be treated on the same footing to avoid large qualitative errors in the splitting of chiral modes.