Adiabatic dynamics of coupled spins and phonons in magnetic insulators

In conventional \textit{ab initio} methodologies, phonons are calculated by solving equations of motion involving static interatomic force constants and atomic masses. The Born-Oppenheimer approximation, where all electronic degrees of freedom are assumed to adiabatically follow the nuclear dynamics, is also adopted. This approach does not fully account for the effects of broken time-reversal symmetry in systems with magnetic order. Recent attempts to rectify this involve the inclusion of the velocity dependence of the interatomic forces in the equations of motion, which accounts for time-reversal symmetry breaking, and can result in chiral phonon modes with non-zero angular momentum even at the zone center. However, since the energy ranges of phonons and magnons typically overlap, the spins cannot be treated as adiabatically following the lattice degrees of freedom. Instead, phonon and spins must be treated on a similar footing. Focusing on zone-center modes, we propose a method involving Hessian matrices and Berry curvature tensors in terms of both phonon and spin degrees of freedom, and describe a first-principles methodology for calculating these. We then solve Lagrange's equations of motion to determine the energies and characters of the mixed excitations, allowing us to quantify, for example, the energy splittings between chiral pairs of phonons in some cases, and the degree of magnetically induced mixing between infrared and Raman modes in others. The approach is general, and can be applied to determine the adiabatic dynamics of any mixed set of slow variables.

Phonons play a crucial role in determining various thermodynamic and electronic properties of materials, including heat capacity, heat transport, electronic conductivity, and superconductivity.Conventionally, phonons are treated within the Born-Oppenheimer approximation [24], i.e., assuming that the electronic degrees of freedom (DOF) can adiabatically follow the motion of the atoms.In these calculations, the potential energy in the phonon Hamiltonian is computed as a function of atomic displacements.The usual harmonic approximation involves keeping only the leading quadratic dependence of the energy on displacement, encoded in the interatomic force constant (IFC) matrix, i.e., the Hes-sian matrix of the energy.Anharmonic treatments go further by taking account of higher-order tensors which describe third and higher derivatives with respect to displacements.These harmonic and anharmonic tensors are all invariant under time-reversal symmetry (TRS), so that at this level of description the phonons are assumed to preserve TRS and to possess the symmetries of the nonmagnetic group.
However, this assumption is incorrect for materials with magnetic order or in the presence of an external magnetic field.Recent efforts have been made to address this issue by incorporating the nuclear Berry potential into the effective Hamiltonian [25][26][27][28][29][30][31][32][33], which arises naturally when the electronic DOF are integrated out under the Born-Oppenheimer approximation, as first pointed out by Mead and Truhlar [25].The nuclear Berry potential introduces a velocity-dependent force into the equations of motion (EOM), which is determined by the nuclear Berry curvature multiplied by the nuclear velocity.
To date, much of the theoretical literature treats models in which the nuclear Berry curvature is an adjustable parameter.The first-principles calculation of the nuclear Berry curvature is still in its infancy, with only a few calculations for molecular [30] and crystalline [33] systems.
Spin-wave excitations, or magnons, constitute additional DOF in magnetic materials.These excitations are chiral from the outset; for example, in the presence of easy-axis anisotropy, spins always precess clock-wise when viewed end-on.While it is not so widely appreciated, magnon dynamics can also be formulated and computed in the context of a geometric-phase framework [50,51].There is no inertial term in the spin dynamics, but the Berry curvature tensor enters the EOM by describing the spin precession in response to a torque [52].This raises the possibility of a uniform treatment of both nuclear and spin degrees of freedom in a common theoretical framework.
Crucially, the frequency (or energy) range of the magnons strongly overlaps that of the phonons.In metals these also overlap with electron-hole excitations, but we shall restrict our attention here to insulators in which there is a clear energy separation between both phonons and magnons on the one hand, and cross-gap electronic excitations on the other.In this case, both types of bosonic excitations can be treated as "slow" DOF, with the remainder of the electronic system following adiabatically.In such cases it is crucial to treat both phonon and spin DOF on the same footing.
A first step was taken in this direction in Ref. [33] for the case of bulk ferromagentic CrI 3 .In that work, a minimal model was proposed in which the nuclear Berry curvature arose entirely from the canting of Cr spins in response to atomic displacement.Though this model captured the essential physics of CrI 3 , a general and quantitatively accurate theory should include other contributions to the Berry curvature including "phonononly" Berry curvature arising from atomic displacements at fixed spin, and mixed "spin-phonon" Berry curvature.Such contributions can play an important role, particularly in predicting energy splittings of high-energy chiral phonons.Additionally, Ref. [33] relied on the input of experimental magnon energies, which may not always be available.
Motivated by the need for better fundamental and quantitative understanding of phonons in TRS-broken systems, we present a generalized adiabatic treatment incorporating all relevant Hessian matrices and Berry curvatures for both phonons and spins.It should be noted that the methodology presented in this work provides a general theoretical framework for treating dynamics beyond the specific case of coupled spins and phonons.I.e., it allows for efficient ab initio calculation of the adiabatic dynamics of any mixed set of slow variables.This goes significantly beyond the more approximate approach of Ref. [33], which we shall refer to as the minimal spin-phonon model below.We also demonstrate an ab initio methodology to calculate these matrices.We conduct four case studies covering ferromagnetic (FM) and antiferromagnetic (AFM) materials in both threedimensional (3D) and two-dimensional (2D) form.We investigate the energy splittings and symmetries of chiral phonons in these systems.
We show that, interestingly, circularly polarized chiral phonons do not always exhibit energy splittings when TRS is broken; whether this splitting occurs depends on the point group symmetry in the presence of magnetic order.Also, terms beyond the IFC are needed to correctly capture the symmetry of phonons in cases where magnetic order breaks a spatial symmetry g, but the combination of g and time reversal T remains a symmetry.This limitation arises because conventional phonon calculations rely on real-symmetric IFC and atomic mass matrices that preserve TRS. 1 We show that correct accounting for phonon symmetries is crucial for determining Raman/IR activities.
This paper is organized as follows.In Sec.II, we establish the theoretical framework via the Lagrangian formalism of adiabatic dynamics, and derive EOM that treat spins and phonons on an equal footing.We also revisit phonon angular momentum and introduce the concept of atom-resolved phonon angular momentum.Section III details the proposed method for calculating all the matrices involved in the present approach, along with the computational details.In Sec.IV, we present our results, which include the analysis of phonons and magnons in 3D FM CrI 3 bulk (Sec.IV A), 3D AFM Cr 2 O 3 bulk (Sec.IV B), and 2D systems (Sec.IV C) such as FM CrI 3 monolayer (Sec.IV C 1) and AFM VPSe 3 monolayer (Sec.IV C 2).Section V provides a discussion of the role of spin-orbit coupling (Sec.V A) and future experimental investigations (Sec.V B).Finally, we summarize our findings and present concluding remarks in Sec.VI.

II. THEORETICAL BACKGROUND
We restrict our considerations to the case of insulating magnetic materials in which there is a clear separation between the energy scales of the phonons and the crossgap electronic excitations.The magnetic order insures that there will also be dynamics associated with spin fluctuations, i.e., the magnons.If the magnon frequencies would be higher than, and robustly gapped from, those of the phonons, it would be possible to treat all electronic excitations, including the magnons, as adiabatically following the phonon DOF.However, this is almost never the case.In the present work, we therefore treat both the lattice and spin DOF on a similar footing, while assuming that there is still a large energy gap between the top of the phonon or magnon spectrum and the onset of cross-gap electronic excitations.
We formulate our theory in the context of a firstprinciples mean-field theory such as DFT, where the dynamics of the nuclei is treated classically while the electronic system is evolved according to the time-dependent Shrödinger equation.This is essentially the domain of time-dependent DFT (TDDFT), but here we aim to treat the spin DOF as slow semiclassical variables alongside the nuclear displacements.This requires a separation of the electronic DOF into a small number of spin DOF and the remaining large number of electronic excitations on the scale of the band gap or above.
To do so, we define a "spin" unit vector on a magnetic ion to be the direction of the average spin density inside a Wigner-Seitz sphere centered at that site, and subsequent calculations of electronic ground states and energies are always computed with these spin variables constrained.The essential requirement is that the remaining electronic system, so constrained, should be free of any remaining slow DOF, i.e., any below-band-gap excitations.The implementation of the Wigner-Seitz sphere constraint is not a serious obstacle in practice, as most DFT code packages have features for carrying out electronic minimizations under the constraint of fixed spin orientations defined in this way.We emphasize that we treat not only the lattice displacements, but also the spin cantings, in a harmonic approximation about the ground-state reference structure.We treat only collinear easy-axis systems here, and assume that the spin cantings with respect to this axis are small.
To finish a discussion of the approximations of our theory, we note that we use "ordinary" adiabatic dynamics, in which the adiabatic perturbation theory is carried only to first order in the rate of change of nuclear or spin variables.This is well justified as long as the gap separating phonons and magnons from interband electronic excitations is large.And finally, we shall shortly make a harmonic approximation, in which the atomic displacements and spin cantings are expanded to leading order around a ground-state reference configuration.
With these understandings, we turn now to a detailed presentation of our methodology.

A. Lagrangian formulation of adiabatic dynamics
While the Hamiltonian formalism is commonly used in the literature to analyze chiral phonons [25,26,30,31,33], we will instead start with the Lagrangian formalism.We show in Sec.II C that this approach is well-suited to developing a comprehensive model of the coupled spinphonon dynamics, while avoiding the difficulties associ-ated with defining canonical momentum for spins.The Lagrangian takes the form where the configuration Q can represent any slow variables.While the formulation is general, we shall focus on the case that Q represents both the nuclear coordinates and the spin variables, where the latter act as constraints on the spin moments as explained above.For clarity of presentation, we assume a finite number of nuclear and spin DOF, as for a molecule or the Γ-point modes of a periodic crystal.The first term in Eq. ( 1) is the kinetic energy associated with the i-th degree of freedom, where M i is the nuclear mass for the phonon variables or zero for the spin-canting variables.The second term is the potential energy, and the third represents the coupling between the time derivative of the adiabatic variable Q i and the Berry potential A i [50].The latter is defined as where |ψ(Q)⟩ represents the electronic wave function at constrained nuclear coordinates and spin orientations Q.
By solving the Euler-Lagrangian equation for the Lagrangian defined in Eq. (1), we obtain where ∂ i denotes ∂/∂Q i .By using Ȧi = j ∂ j A i Qj , we can simplify Eq. (3) to Here where Ω ij is the Berry curvature, and is therefore gaugeinvariant.Although the gauge-dependent quantity A i (Q) appears in the Lagrangian, the EOM are gauge-invariant since only G ij (Q) appears.
In this work, we focus on small oscillations near equilibrium.To simplify the analysis, we introduce the generalized displacement vector q defined via Q is the equilibrium value of the i-th degree of freedom.Introducing the harmonic approximation, we expand the potential energy ϵ(Q) in terms of q as where K ij is the Hessian matrix in q i and q j , i.e., The EOM for q is then given by where G = G(Q)| Q=Q0 is also computed at the reference configuration Q 0 .Conventional treatments of phonons in isolation typically use only mass and force-constant matrices M and K, while spin dynamics in isolation is described by the anisotropy tensor K and Berry curvature G.Note that M is real diagonal, K is real symmetric, and G is real antisymmetric.
To determine the frequencies, we substitute q i (t) = e −iωt q i into Eq.( 7), yielding where |q n ⟩ is a column vector with the i-th component q n,i corresponding to the n-th mode associated with DOF i. Equation ( 8) is easily solved using, e.g., the methods of Sec.III C. The above treatment provides a semiclassical theory of the adiabatic dynamics of the system.
B. Mead-Truhlar approach without explicit spin degrees of freedom If we limit the slow variables Q to include only atomic coordinates, i.e., allowing all electronic degrees of freedom (including spins) to be in their instantaneous ground state for a given Q, we restore the treatment of Mead and Truhlar in Ref. [25].To indicate the specialization to atomic coordinates, in this section we replace Q and q by R and u, where R denotes the equilibrium position for atomic coordinates, and u denotes the atomic displacement from the equilibrium.Now the EOM for u in the harmonic approximation is where l and m are composite indices for Iα, I runs over atoms and α represents a Cartesian direction.Equation (9) shows that G MT mn un corresponds to a force acting on coordinate m that is proportional to the velocity of coordinate n.An alternative derivation of Eq. ( 9) using the quantum theory in the Hamiltonian framework is given in Appendix A.
The conventional treatment of phonons [1,2] is recovered by discarding the term involving the nuclear Berry curvature G MT in Eq. ( 9).This is justified in TR-invariant systems, where G MT vanishes by symmetry.However, G MT is often neglected even when the system is not TR symmetric, with the consequences that Eq. ( 9) obeys TRS and the phonons will not have the correct symmetry of the TRS-broken system.This will force modes with opposite chirality to be degenerate at the zone center, which is not necessarily the case in a magnetic material.
Recent works have used Eq. ( 9) to demonstrate splitting of chiral modes as a result of TRS breaking [31,33].However, it was demonstrated in Ref. [33] that for CrI 3 , the main contribution to G MT comes from canting and precession of spins.As mentioned earlier, the energy scale for spin rotations corresponds to the frequency of the magnons, which is close to that of the phonons in most systems.Consequently, the assumption that atomic displacements are the only slow DOF in the system, which led to Eq. ( 9), is not valid.Reference [33] developed a Hamiltonian formalism for coupled spin-phonon dynamics that we refer to here as the "minimal spinphonon model."In the next section we will develop a more general Lagrangian-based approach that is well suited to treating phonon and spin dynamics together.

C. Treatment of spins and phonons on the same footing
We now return to the framework of Eq. ( 7) in which q i includes both nuclear and spin DOF.The Euler-Lagrange EOM derived from Eq. ( 7) where i runs over both phonon and spin DOF.In the phonon sector, q i = u Iα is a shorthand for a small displacement of atom I in Cartesian direction α = {x, y, z}.
In the spin sector, q i denotes a small spin canting q i = s Jβ , where J runs only over magnetic ions and β indexes the spin canting in the two directions orthogonal to the ground-state spin orientation.Specializing to easy-axis systems with spin axis along ±ẑ, we let β run over only the two in-plane Cartesian directions.That is, q i = s Jβ = S Jβ /|S J | describes the component β = {x, y} of the unit vector of spin S J located on the J-th magnetic ion.
The matrix M ij in Eq. ( 10) is the diagonal mass matrix introduced in Eq. (1), with zero entries for the spin DOF; K ij = ∂ i ∂ j ϵ(q) is a generalized Hessian matrix; and embodies the Berry curvature as in Eq. (5).In these last two expressions, ∂ j denotes ∂/∂q j , which can be a derivative with respect to either nuclear or spin-canting coordinates, evaluated at the reference ground-state configuration.It should be emphasized that in our formulation, the Hessian matrices are expanded to quadratic order in phonon or magnon amplitudes, resulting in a harmonic theory with infinite lifetimes for both phonons and magnons.This is also evident in the hermiticity of Eq. (10).While including anharmonic interactions beyond quadratic order would allow decay of a phonon or magnon excitation into two or more lower-energy excitations, thereby rendering their lifetimes finite, such effects are not considered in this work and remain an area for future investigation.
In this context, derivatives with respect to nuclear coordinates must be taken at fixed spin.That is, the ma-trices K and G are now computed in terms of the electronic quantum state |ψ(R, s)⟩, rather than |ψ(R)⟩ as in Sec.II B. As explained earlier, this requires a calculation of the electronic ground state subject to constrained spin orientations.While there is some freedom in the definition of the spin unit vector, we follow the established approach of defining it in terms of the integrated spin density inside a Wigner-Seitz sphere, as discussed in Sec.III B. This choice encodes the distinction between "spin" and "other electronic" DOF in our theory.
To simplify the analysis, we can partition all DOF i into phonon DOF (labeled as p) and spin DOF (labeled as s).The matrices in Eq. ( 10) can then be represented using a block structure as We use the term "bare phonons" to refer to phonons that are calculated without considering G (pp) , G (ps) (G (sp) ), or K (ps) (K (sp) ), while phonons calculated with the inclusion of those terms are referred to as "perturbed phonons."The term "perturbed" is used in recognition of the fact that the influence of these terms is generally small, although we solve Eq. ( 10) exactly.Nevertheless, we also provide a perturbation analysis of K (sp) , G (sp) , and G (pp) in Appendix D, illuminating the physical implications of each term, and thereby offering valuable insight.
It is worth noting that G (pp) , G (ps) (G (sp) ), and K (ps) (K (sp) ) are zero in the absence of spin-orbit coupling (SOC) in collinear systems, a point elaborated upon in Sec.V A below.In that case, the broken TRS in the spin sector is never communicated to the orbital electronic sector or, in turn, to the phonon sector.
Importantly, in the case that K (ps) (K (sp) ) and G (ps) (G (sp) ) vanish, phonons and spins decouple, and Eq.(10) reduces to the EOM for phonons and magnons separately.In the phonon sector, it reduces to Eq. ( 9) which corresponds to the approach proposed by Mead and Truhlar [25].Meanwhile, in the magnon sector, it reduces to which aligns with the EOM for magnons presented in Ref. [50].Additionally, when Ω (ss) is simplified to consider only isolated spinors, rendering inter-spin elements negligible, Eq. ( 12) reduces to the well-known Landau-Lifshitz equation [52,53].Further, if the energies of phonons are considerably lower than those of magnons, Eq. (10) in the phonon sector also reduces to the Mead-Truhlar approach, as discussed in Ref. [33].The present approach defined by Eqs.(10) and (11) uses Hessian matrices and Berry curvature tensors in terms of all DOF.In contrast, the minimal spin-phonon model presented in Ref. [33] included only the spin-spin component of the Berry curvature, so that G (pp) , G (ps) , and G (sp) were assumed to vanish.Furthermore, the minimal spin-phonon model in Ref. [33] only includes one bare phonon doublet and one bare magnon.In Appendix B, we will discuss a more comprehensive version of the spin-phonon model which incorporates all bare phonons and magnons but still neglects G (pp) , G (ps) , and G (sp) , and discuss its connection to the minimal spinphonon model and its equivalence with the widely employed Landau-Lifshitz equation [52].We demonstrate that our approach, through the inclusion of non-trivial Ω (pp) , K (sp) , and Ω (sp) terms, provides a fruitful generalization of the Landau-Lifshitz equation.

D. Phonon angular momentum
Before discussing the first-principles methodology to calculate the terms in Eq. (10) and Eq. ( 11), we review the definition of phonon angular momentum and introduce the concept of atom-resolved phonon angular momentum, as this will be important for characterizing the chiral modes in Sec.IV.The definition of phonon angular momentum can be found in the literature [26], and we briefly revisit the relevant definitions here.First, we note that by solving Eq. (10), one can obtain an energy eigenvalue ω n and the corresponding mode with both phonon and spin components, i.e., |q n ⟩ = |u n ⟩⊕|s n ⟩, where n runs over different solutions of Eq. (10).Since we are mainly interested in the phonon sector in the present work, we define the atom-resolved phonon angular momentum using the phonon part |u n ⟩ of the mode.
The solutions to the EOM described by Eq. ( 10) yield energies that differ from those of the bare phonons.However, in all of the systems examined in this work, the differences between these energies and those of the bare phonons are relatively small.As such, we can identify the mode |q n ⟩ as"phonon-like" if its energy ω n is in close proximity to that of a bare phonon.Furthermore, we note that for a phonon-like mode |q n ⟩, the |s n ⟩ contributions are very small compared to |u n ⟩.In systems where the frequencies of the phonons and magnons coincide, the aforementioned conditions may not be met; nonetheless, in all cases we consider in this work, the zone-center phonons and magnons of relevance exhibit distinct energies, thus ensuring well-defined phonon-like and magnon-like modes.
For phonon-like modes, we continue to adopt the normalization convention even though ⟨u m |M (pp) |u n ⟩ is no longer exactly zero for m ̸ = n.It is also possible to have a "perturbed magnon" solution, which is a magnon-like solution with tiny phonon components.
The definition of phonon angular momentum was originally proposed in Ref. [26].For a phonon-like mode |u n ⟩, the atom-resolved phonon angular momentum (ARPAM) L n,Iz , for atom I of mass M I in the z direction is with L n,Ix and L n,Iy defined similarly by cyclic permutation of Cartesian indices.Note that if the projection of a mode vector on a given atom is of the form x + iŷ, the atom undergoes a counterclockwise rotation when viewed from above and contributes a positive L z .The total phonon angular momentum (PAM) 2 L n,z is defined as which is the sum of the angular momenta L n,Iz over all atoms I.The bare phonons solve the secular equation involving only M (pp) and K (pp) , which is equivalent to the conventional phonon treatment of Eq. ( 9) with G = 0.These bare phonons can always be chosen real, in which case it follows from Eq. ( 14) that the full PAM and individual ARPAM always vanish.In the case of degenerate modes it may be possible to choose chiral linear combinations, but the trace over the degenerate set of modes always results in zero PAM and ARPAM.This is a consequence of the fact that TRS has not yet been broken at the bare level of description.

A. Finite difference method
This section demonstrates how to compute all matrices in Eq. ( 11) using finite-difference methods in the context of first-principles calculations.The nuclear masses M (pp) are trivially known.All calculations are carried out using the DFT methodology described in Sec.III B at fixed atomic coordinates and fixed spin orientations, where the latter are defined in terms of an integration of the spin density over a Wigner-Seitz sphere as mentioned earlier. 2 We are aware that in some literature, the PAM is referred to as "pseudo angular momentum" [45], which is a non-zero integer multiple of ℏ if the phonon respects Cn symmetry with a nontrivial eigenvalue.However, it is important to note that the phonon angular momentum discussed in this paper is the kinetic angular momentum and is not conserved in the absence of infinitesimal rotational symmetry in the lattice.In contrast, the pseudo angular momentum is conserved up to nℏ if the system possesses Cn symmetry.In Sec.IV, we refer to the C 3 eigenvalue as χ(C 3 ).
The force-constant matrix K (pp) , which is also the Hessian matrix of the energy with respect to nuclear DOF, is defined as Here l and m run over nuclear DOF and ϵ is the total energy ϵ(s, R) of the configuration s, R. The matrix element lm is computed by taking the finite difference of the Hellmann-Feynman forces F l [1,2] while constraining the spin directions to lie along ẑ.
The Hessian matrix of the energy with respect to spin DOF is denoted by K (ss) , with elements defined as where a and b run over spin DOF.To calculate K (ss) , we compute the second derivatives of the total energy with respect to small canting of the spins.Specifically, letting s a = s Iα and s b = s Jβ , for each pair (I, J) we compute the energies of (α, β) = ±(0.02,0), ±(0, 0.02), and ±(0.02, 0.02) relative to the ground state, while constraining all other spin moments to remain along ẑ.
The spin-phonon Hessian matrix K (sp) is defined as where a and l respectively run over the spin and nuclear DOF.K (sp) al is computed by taking the finite-difference derivative of the force F l with respect to the spin coordinate s a .
An alternative approach to calculating K (sp) was described in Ref. [33].Near the ground state, the energy ϵ of the entire system can be expanded as We define a spin response matrix χ (sp) as where the second equality holds if both s a and u l are small.Then we can obtain where we have restored the matrix form for simplicity.In practice, one can perturb the ground state structure with u l by manually moving the atoms from the equilibrium position and calculate the spin canting s a with respect to u l .We used Eq. ( 22) to calculate K (sp) for bulk CrI 3 , and the result is consistent with Eq. ( 18).However, due to the slow convergence of spin relaxation, we recommend using Eq. ( 18) to calculate K (sp) .In Eq. ( 11), G is simply ℏ times Ω, where Ω is the Berry curvature.The latter is computed using Stokes' theorem as expressed by where i and j run over all DOF, and Φ ij is the Berry phase around a diamond-shaped region of parameter space whose area appears in the denominator.Specifically for, e.g., a finite system with electronic ground state wave function |ψ⟩, In an extended crystal with a single occupied band, one must sum over the Bloch wave vector k in the Brillouin zone to obtain ij is defined as in Eq. ( 24) but with ψ replaced by the Bloch function ψ k .Φ ij now has the interpretation of a Berry phase per unit cell, consistent with the interpretation of the Hessian K as an energy per unit cell.To extend Φ (k) ij to the multiband case, we can replace the inner product of two Bloch states with the overlap matrix in the usual way [54] as Here, the overlap matrices are defined as where m and n are band indices.When computing the matrices, a finite difference of 0.015 Å is used for the phonon DOF, and 0.02 is used for the spin DOF.We take care to choose these finite differences to ensure they remain within the linear regime.The numerical values of K (ss) and G (ss) for all four studied materials can be found in Appendix C.
All structures are relaxed using the local spin-density approximation (LSDA), and the convergence criteria for forces and energies are 10 −3 eV/ Å and 10 −8 eV, respectively.After relaxation, the wave functions are calculated using static calculations with a convergence criterion of 10 −10 eV for energies.Spin-orbit coupling, which is essential to the physics described in this work, is included in all static calculations except structural relaxations.

C. Solution of the equations of motion
In this section, we introduce a practical method for solving the EOM in Eq. (10).Given that the generalized mass matrix is non-invertible, conventional approaches are not directly applicable.We proceed as follows, although other potential methods may be suitable.
First, we rewrite Eq. ( 10) as which can be shown to reduce to Eq. ( 10) with the definitions The eigenvalues (multiplied by i) of the matrix in Eq. ( 27), denoted as ω n , correspond to the solutions of Eq. (10).We restrict the index n to run only over solutions for which ω n > 0, which we take to be the physically meaningful ones.
In practice we find that small numerical errors remain in the eigenvectors, which arise from the fact that a Hermitian eigensolver cannot be applied in the context of Eq. ( 27).We have found that these errors can easily be removed by following with a second step in which we substitute the computed energies ω n back into Eq.( 10) and then apply a Hermitian eigensolver to recalculate the eigenvectors.

IV. RESULTS
We conduct four case studies covering both FM and AFM systems, in both 3D and 2D.In Sec.IV A, we present the results obtained from our present approach for bulk CrI 3 , focusing on chiral phonons.The significance of the Berry curvatures neglected in Ref. [33], is discussed in Sec.IV A 1, while the relevant solutions for the magnons are discussed in Sec.IV A 2. Subsequently, in Sec.IV B, we investigate the phonons and magnons in bulk Cr 2 O 3 using the present approach.Shifting our focus to monolayer systems in Sec.IV C, we first consider 2D FM CrI 3 in Sec.IV C 1, and then 2D AFM VPSe 3 in Sec.IV C 2.

A. Chiral phonons in bulk ferromagnetic CrI3
CrI 3 is a hexagonal van der Waals material that exhibits FM order in both bulk and monolayer phases [68].It is an insulator with inversion and three-fold rotational symmetry around the z-axis in both phases.The magnetic moments on the Cr ions do not break the inversion symmetry.The symmetries will be discussed in detail below.
The crystal structure of the bulk CrI 3 unit cell is depicted in Fig. 1.For bulk CrI 3 , the structural and magnetic symmetries are identical (space group R 3, point group S 6 ), regardless of the presence of FM order.That is, the magnetic space group is of Type I ("colorless"), in which no symmetry operations involve TR [69].Therefore, the zone-center phonons can be categorized into the irreducible representations (irreps) of the S 6 point group as 4E g ⊕ 4E u ⊕ 4A g ⊕ 4A u .Among these, one A u mode and one pair of E u modes are acoustic modes.In this work, we will focus on the optical modes.
The E g and E u irreps are complex-conjugate irreps, which are more properly decomposed further into onedimensional irreps of opposite chirality.In the presence of TRS, or when neglecting the term involving the Berry curvature G tensor in Eq. ( 9), the two modes making up one of these E irreps are degenerate.However, this degeneracy is broken by the FM order, which induces a nonzero Berry curvature via the SOC as discussed Sec.II.
We first present our results for the case of the E g and E u phonons in Table I.The energy shifts ∆E represent the differences between the bare phonon frequencies E 0 computed using only the force-constant term K (pp)  in Eq. ( 9), and the modified energy obtained from the present approach of Eq. ( 10) or from the spin-phonon model of Appendix B. 3 The PAM L z of each phonon is also reported for the present approach.
From Table I we observe that all E g and E u bare phonons are doubly degenerate, while the degeneracy is broken for the perturbed phonons.Although our numerical solution for the bare modes initially yields a pair of real phonons |u 1 ⟩ and |u 2 ⟩, we resolve these into eigenstates of the C 3 symmetry operator.We first ensure that ⟨u 2 |C 3 |u 1 ⟩ is positive, flipping the sign of |u 2 ⟩ if not, and then construct the bare chiral modes |u ± ⟩ = (|u 1 ⟩∓i|u 2 ⟩)/ √ 2 belonging to eigenvalues ε = exp(i2π/3) and ε * = exp(−i2π/3) respectively.We shall designate these as '+' and '−' modes, and refer to them as belonging to the ε and ε * symmetry sectors, respectively.The former (latter) are characterized by a clockwise (counterclockwise) rotation of the Cr atoms when viewed from above.
The numerical solutions for the perturbed phonons automatically generate chiral C 3 eigenstates, and we find that each of these is almost identical to the bare chiral mode of the same symmetry that is closest in energy.There is only a small admixture of other bare modes belonging to the same sector.This allows a straightforward association of bare and perturbed modes as shown in Table I.
Table I demonstrates that the E u phonons with ε chirality have the largest energy corrections.This follows from their close energy proximity to the E u magnons, which also possess ε chirality.A detailed discussion of magnons in bulk CrI 3 is provided in Sec.IV A 2, and a perturbative treatment of phonon-magnon interactions is presented in Appendix D. The pronounced effect on the E u phonons is further evident in Eq. (D20), where the phonons with ε chirality have the smallest energy denominators.
We briefly address the expected uncertainties in our first-principles calculations of the K (sp) , G (sp) , and G (pp) matrices, as these are the matrices that contribute to the energy shifts (see Appendix D).First, we analyze the fitting errors for K (sp) associated with our finite-difference methodology, which is done via linear regression of the forces against spin canting angles.We find that the norm of the fitting error is roughly 1.5% of the norm of K (sp) itself.A second useful metric is the norm of the residual resulting from symmetrization of a matrix relative to the norm of the matrix itself, which we find to be 0.4%, 2.0%, and 2.7% for K (sp) , G (sp) , and G (pp) respectively.The perturbation analysis of Eqs.(D20) and (D32) indicates that the energy shifts and splittings are quadratic in K (sp) and G (sp) and linear in G (pp) .This implies overall errors on the order of 3-4% in the values reported above for these shifts and splittings.
In Fig. 2, we visualize the real and imaginary parts of two E g chiral phonons around 7 meV, denoted as E (1) g for the mode with lower energy and E (2) g for the mode with higher energy.It is apparent that those two chiral phonons are approximately complex conjugates of each other.
Note that the argument about energy splitting based on point group irreps can be generalized to any magnetic material with a point group that includes an E-type irrep (E, E g , or E u ) consisting of two 1D irreps when complex representations are allowed, indicating that similar phonon energy splittings can be expected in such materials.
The phonon angular momentum is calculated using Eqs.( 14) and (15), and the results are listed in Table I.Due to the C 3 symmetry, the total angular momentum of a given chiral phonon can only have a z-component.We find values of |L z | ranging from 0.13 to 0.95 in units of ℏ, indicative of a substantial angular momentum for many of the chiral phonons, especially the higher-frequency ones.
We find that the bare circular modes, constructed as 2 from the real bare modes, can be used to evaluate the PAM to a good approximation.The results are close to those shown in Table I, with a typical error of 10 −3 ℏ in |L z |.However, when calculating the PAM in this way, the angular momentum vanishes exactly when summed over any pair of chiral phonons, which have equal and opposite L z values.In contrast, solutions using the full matrices of Eq. (10) reveal that the cancellation is no longer perfect, and the total angular momentum of a pair is small but non-zero.The slight discrepancy between the two |L z | values arises from the magnon-mediated phonon-phonon mixing, which appears at second-order in the phonon-magnon coupling and therefore makes only a minor contribution.We provide a more detailed discussion on this topic based on a perturbation approach in Appendix D 1.The methodology of Ref. [33] for CrI 3 differs from the approach of this work (i.e., Eq. ( 10) and Eq. ( 11)) in the following ways.First, the only Berry curvature considered in Eq. ( 11) was G (ss) .Second, the minimal model in Ref. [33] neglects the mixing between bare phonon doublets, considering the interaction which each doublet and one magnon branch separately.In Appendix B, we discuss a more comprehensive version of that model which still neglects G (pp) , G (sp) , and G (ps) , but incorporates all bare phonons and magnons.To show the importance of G (pp) and G (sp) (G (ps) ) in accurately determining energy splittings, we also apply the model present in Appendix B to calculate chiral phonon energies for E g and E u modes, and the results are present in Table I.
As shown in Table I, both methods exhibit an energy splitting of chiral modes compared to bare phonons.The energy splitting of low-energy E g and E u phonons is similar to that obtained using the present approach, suggesting that G (pp) and G (ps) have little effect on those modes.However, for high-energy E g and E u phonons, the spinphonon model results differ significantly, indicating the importance of G (pp) and G (ps) for these phonons.Specifically, the absolute values of the G (pp) matrix elements for these four E g modes are {0.995,0.254, 0.972, 5.623} in units of µeV.The highest-energy E g modes exhibit substantially greater values, leading to a pronounced energy splitting.Appendix D 3 details the perturbative analysis of G (pp) , with Eq. (D32) specifying the energy corrections attributed to the G (pp) term.As indicated by Eq. (D32), G (pp) primarily accounts for the substantial splitting observed in the highest-energy E g mode, with both direct computation and perturbative treatment yielding a correction of ±0.0028 meV.
Interestingly, we find that the PAM calculated using the spin-phonon model (not shown) is in good agreement with the full method with an error of the order of 10 −4 ℏ, suggesting that the phonon vectors are not very different regardless of whether G (pp) and G (ps) are included or not.
We now turn to a consideration of the A g and A u modes.The spin-phonon model presented in Appendix B has no effect on the A g and A u phonons, since K (sp) FIG. 3. Visualization of the real and imaginary components of the perturbed phonon A gives the coupling to magnons, and there are no zonecenter magnons in these symmetry sectors.Considering the present approach, although all (ss) and (sp) terms are absent in the A g and A u sectors, the G (pp) matrix does not vanish, and causes a slight mixing between two different A g (or A u ) bare phonons.That is, we can write {n,m} ⟩ represents the bare phonons from A g (or A u ) irrep.Here, iδ m is purely imaginary and expected to be small, attributed to the minor size of G (pp) matrix element, as discussed in Appendix D 3.
This results in a nonzero PAM and nonzero energy shift of each phonon relative to its bare energy.Although these phonons exhibit non-zero PAM, we refrain from referring to them as "chiral phonons" because they belong to the χ(C 3 ) = 1 sector of modes that are invariant under the C 3 operation.Instead, we refer to them as (1 + iδ)-type perturbed phonons, or simply perturbed phonons.
The energy shifts and PAM for A g and A u perturbed phonons are presented in Table II, and the real and imaginary parts of the A g perturbed phonon near 10 meV (A (1) g ) are visualized in Fig. 3.

Bare and perturbed magnons
So far we focused solely on the phonon-like solutions of the equation of motion given by Eq. (10) in the present approach (or Eq. (B1) for the spin-phonon model).Nevertheless, these equations also admit magnon-like solutions, which we call "perturbed magnons."In this sec- tion, we investigate the impact of phonons on the magnon spectrum, using bulk CrI 3 as a case study.The numerical values for the matrices K (ss) and G (ss) are provided in Appendix C.
The EOM for a bare magnon without any phononmagnon interaction is given by Eq. ( 12).The energies of bare magnons are listed in Table III.The E g magnon corresponds to the acoustic mode, where the two Cr spins have the same canting, while the E u magnon corresponds to the optical mode, where the two Cr spins have opposite canting.Both magnons have positive energies and belong to the '+' sector with χ(C 3 ) = ε, indicating that the two Cr spins are rotating clockwise. 4he perturbed magnons have slightly modified energies (Table III).The energy change from the bare magnon for the E u mode is more significant compared to the E g magnon, since the E u phonons are closer to the E u magnon in the energy spectrum.By far the most significant contribution to the energy shifts is the inclusion of K (sp) ; neglecting the G (pp) , G (sp) , and G (ps) terms changes the energies by less than 0.4 µeV.A perturbation treatment of magnon energies is provided in Appendix D 1.
Note that the equations of motion have two additional negative-energy magnon solutions belonging to the '−' sector with χ(C 3 ) = ε * .These correspond to the counterclockwise precession of spins and are not physically observable.Nonetheless, they still influence the system dynamics through their interaction with physical '−' phonons, which acquire some magnon dressing in which the spin vectors are forced to precess in the unnatural counterclockwise sense.This becomes clear from the perturbation analysis presented in Appendix D, where the summation index µ in Eq. (D16) and Eqs.(D19-D20) runs over all solutions of Eq. (D3), including those of negative energy.However, since the energy denominators in Eq. (D20) are larger when coupling to negative-energy solutions, the magnon dressing is typically smaller.This explains why the '+' phonons, which couple to positiveenergy magnon solutions, are more strongly perturbed than the '−' ones, which do not.This can be seen in Table I, where it is especially noticeable for the E u modes.

B. Phonons in antiferromagnetic Cr2O3
In this section, we present our results for the phonons in bulk Cr 2 O 3 .As shown in Fig. 4, Cr 2 O 3 has an AFM insulating ground state with antiparallel magnetic moments aligned along the threefold-rotational z axis.
The magnetic symmetry of Cr 2 O 3 is richer than that of CrI 3 , in that it belongs to a Type III ("black-white") magnetic group [69].Since we are interested in zonecenter phonons, we frame our discussion in terms of the magnetic point group G.The black-white character means that half the the elements of G come with the TR operator T and half come without.The latter form a "unitary subgroup" H, and the others are of the form gH where g is one of the antiunitary operators in G.The "structural point group" Ḡ contains the same list of operators as in G, except that T is removed from all the antiunitary operators.The important point in what follows is that while the Hessian matrices are even under all elements of Ḡ, the Berry curvature matrices are even under the elements of H and odd under the elements of Ḡ − H.As a consequence, when the Berry curvature is included, the irreps of unitary subgroup H, not the structural group Ḡ, should be used to label the perturbed modes of this system.
In the case of Cr 2 O 3 , the magnetic space group is R 3′ c ′ , and the structural point group is D 3d .The magnetic ordering on the Cr sublattice breaks inversion (i), dihedral mirror (σ d ), and rotoinverion (S 6 ) symmetries, reducing the unitary point group to D 3 , whose irreps will be used to label our perturbed phonons and magnons.The corresponding time-reversed operators iT , σ d T , and S 6 T are present in the magnetic point group, but being antiuni- tary, do not induce any additional irreps.
While the experimental magnetic moment of Cr 2 O 3 is found to be along the z direction, DFT calculations predict an in-plane magnetic moment.This discrepancy can be resolved by applying a 2% epitaxial strain, which results in an easy-axis ground state [62].To ensure consistency with the experimental ground state, we apply a 2% epitaxial strain in our calculations, as our model assumes spins to be oriented along the z-direction.
As previously discussed, the irreps of point group D 3 should be used to characterize zone-center perturbed phonons, which can be decomposed into 4A 1 ⊕6A 2 ⊕10E.Among them, one A 2 and two E modes correspond to acoustic phonons.According to the character table of D 3 , the E irreps have to be 2D irreps, indicating that there is no degeneracy breaking or energy splitting compared to bare phonons.This is confirmed by the numerical results of perturbed phonon energies, which are presented in Table IV.It is worth noting that the perturbed phonon energies differ from those of the bare phonons, which are also included in Table IV for reference purposes.

Mixing of phonon irreps
As mentioned earlier, although i, σ d , and S 6 are no longer symmetries, iT , σ d T , and S 6 T remain symmetries of the system.On the other hand, the matrices K (pp) and M (pp) are both real.Applying the symmetry operations iT , σ d T , and S 6 T to these matrices is equivalent to applying i, σ d , and S 6 , respectively.Therefore, K (pp) and M (pp) are symmetric not only under the D 3 symmetry operations, but also under i, σ d , and S 6 , resulting in a point group symmetry of D 3d .
This has an important impact on the symmetries of zone-center bare phonons, which are calculated using K (pp) and M (pp) : bare phonons at the Γ point must have D 3d symmetry, and thus can be labeled by the irreps of D 3d .The Γ point bare phonons decompose into , where one E u doublet and one A 2u mode are acoustic modes.The irreps of the bare phonons are included in Table IV.
The presence of magnetic moments in bulk Cr 2 O 3 results in the breaking of i, σ d , and S 6 symmetries associated with the bare phonons.Consequently, the original D 3d point group symmetry is reduced to D 3 , leading to inter-irrep mixing between bare phonons of different irreps.Referring to the correlation table of the D 3d point group, we find that the E g and E u irreps combine to form the E irrep, the A 1g and A 1u combine to form the A 1 irrep, and the A 2g and A 2u irreps combine to form the A 2 irrep of D 3 .An important question arises: How significant is this inter-irrep mixing, and to what extent does the spatial symmetry breaking due to the magnetic order impact the behavior of the phonons?
To address this question, we project the bare phonons belonging to irreps of the D 3d group onto the perturbed phonons.This allows us to determine the components of each irrep in D 3d that contribute to the given perturbed phonon.We introduce the concept of irrep decomposition component, denoted as ρ n (irrep), which is defined as where n is the label for the perturbed phonon mode and m runs over bare phonons from a particular irrep.The last column of Table IV reports our results for ρ n (E u ) for n ∈ E g and vice versa, and similarly for A 1g -A 1u and A 2g -A 2u components.In each case, the majority component (e.g., ρ n (E g ) for n ∈ E g ) is almost unity, so we list only the minority components.These are non-zero as expected, but we find them to be quite small, typically between 10 −3 and 10 −4 .This implies that the phonon sector exhibits only weak breaking of inversion symmetry, allowing us to continue to refer to the perturbed phonons as E g -like or E u -like.The experimental implication of the inter-irrep mixing is that the perturbed phonons will exhibit distinct Raman and IR activities compared to bare phonons, as summarized in Table V. Above the Néel temperature, where TRS is preserved, the symmetry of the phonons corresponds to that of the bare phonons, so that E u phonons are IR-active but Raman-inactive.However, upon cooling the sample below the Néel temperature, E u phonons undergo mixing with E g phonons, resulting in the emergence of E u -like perturbed phonons.These E u -like phonons possess Raman activity (in addition to IR activity), as they technically belong to the E irrep of the D 3 point group.However, since the inter-irrep mixing is small, the Raman activity of the E u -like phonons is relatively weak.Nonetheless, recent Raman measurements have confirmed the presence of these features [70].In a similar way, the Raman-active E g -like phonons acquire some small IR activity.Clearly, an approach such as ours, which treats the coupling of phonons and spins in a realistic and symmetry-consistent manner, is needed to describe these effects.

Atom-resolved phonon angular momentum
As mentioned earlier, inversion is no longer a symmetry in Cr 2 O 3 , but inversion times time reversal (iT ) remains a symmetry.This means that iT maps total angular momentum from ⃗ L to − ⃗ L, resulting in ⃗ L = ⃗ 0 for any non-degenerate single mode or for the sum over two degenerate modes.However, each atom can still possess a nonzero atom-resolved phonon angular momentum (ARPAM), defined by Eqs. ( 14) and ( 15) above.ARPAM is a pseudovector assigned to each atom and has the same symmetry as a local magnetic moment.The configuration of the local magnetic moment is determined by the magnetic space group, and this is also true for the ARPAM.For Cr 2 O 3 , Cr and O atoms occupy Wyckoff positions 12c and 18e respectively.The possible configurations of ARPAM (and local magnetic moments) are listed in Table VI 5 , indicating that Cr atoms can only have out-ofplane angular momentum L z or −L z , while O atoms can only have in-plane angular momentum in a way that respects C 3 symmetry.This is further supported by the numerical results presented in Table IV.In Fig. 5, we visualize the ARPAM for the E It is important to note that there exists a gauge freedom for each E doublet, where the two degenerate modes can be unitarily mixed with each other through a U (2) matrix.Therefore, discussing ARPAM for each individual E mode is meaningless, as it is gauge-dependent.However, the sum of ARPAM for the two degenerate E modes is gauge-independent.To demonstrate this, let us consider an E doublet labeled by n = {1, 2}, where L n,Iz denotes the ARPAM for atom I along the z-direction.Then, the sum of L 1,Iz and L 2,Iz is which is the trace of a product of three matrices.While different gauges correspond to different bases for these matrices, they do not affect the trace, making L 1,Iz + L 2,Iz a gauge-invariant quantity.This is also true for the x and y directions.The ARPAM for E-irrep perturbed phonons shown in Table IV is traced over the two degenerate modes in the doublet.For E modes with nonzero ARPAM, there is no gauge choice in the degenerate subspace where the phonon mode displacements can be expressed as purely real eigenvectors.

Chiral decomposition of E modes
Despite the fact that perturbed E phonons are always doubly degenerate, it is still worthwhile to decompose the doublet into two single E chiral phonons that respect C 3 symmetry individually and have complex eigenvalues.There are several reasons for this.Firstly, these chiral phonons are excited by circular polarized photons [45].Secondly, they exhibit different energies in the presence of an external magnetic field along the z-direction.Finally, these modes correspond to a special gauge choice in which each individual mode has the largest L z magnitude.To decompose the E doublets, we diagonalize the C 3 matrix using the bases formed by two degenerate E modes.The resulting eigenvectors correspond to two E chiral phonons that respect C 3 symmetry individually.Since the two E modes together respect the C 3 symmetry, the C 3 matrix has to be a 2 × 2 unitary matrix, and is therefore diagonalizable.
We have decomposed the E g -like perturbed phonons around 36 meV into two chiral phonons, denoted as E (1)+ g and E (1)− g , and have visualized their real parts, imaginary parts, and ARPAM in Fig. 6.The real parts of the two modes are almost identical to each other, while the imaginary parts are almost opposite.As a consequence, the ARPAMs of the two modes are also nearly opposite to each other.The PAMs for the '+' and '−' modes are ±0.8593ℏ,respectively.However, the ARPAMs of the two modes do not cancel out perfectly, and their sum is shown in the left panel of Fig. 5.Other E g -like and E ulike perturbed phonons can also be decomposed in the same manner.

Magnons
We now turn to a discussion of the magnons.As a reminder, bare magnons are solved as described in Sec.IV A 2, while perturbed magnons are magnon-like solutions of the EOM of the present approach.The energies of bare and perturbed magnons are presented in Table VII, and the numerical values for the matrices K (ss)  and G (ss) describing the bare magnons are given in Appendix C.
We find that bare and perturbed magnons have almost identical energies, which results from the weak SOC for both the Cr and the O atoms.Since magnons belong to the E irrep, they must be doubly degenerate.These degenerate magnons can be decomposed into two magnons, each with χ(C 3 ) equal to ε or ε * , following the same method we used with E phonons.Additionally, these magnons can also be excited by circular polarized photons.
We next analyze the symmetries of magnons by considering the E g and E u components of each magnon mode, in a manner analogous to Eq. ( 29).Unlike the bare phonons, however, the bare magnons do not have well- defined parity, because the magnetic order strongly violates inversion symmetry.Therefore, we base our analysis on the eigenvectors of the anisotropy matrix K (ss) instead.Note that K (ss) , being quadratic in the spin DOF, is real and symmetric, so that the symmetry operator iT behaves like i for K (ss) .Thus, it has well defined E g and E u eigenvectors that we denote as |t µ ⟩.Then the E g and E u components of a general solution |s µ ⟩ are defined via where both |t ν ⟩ and |s µ ⟩ are normalized.These ρ values can be used to quantify the extent to which the inversion symmetry is broken for magnons, and the numerical values are listed in Table VII.We find that the acoustic magnons, which have lower energies, are nearly E u modes.However, the optical magnons have comparable ρ(E g ) and ρ(E u ) values, indicating that the inversion symmetry is strongly broken for these modes.This is not surprising, and is true also for the bare magnons, since the equations of motion for spin strongly violate TRS.Note that we apply a 2% epitaxial strain to fix the sign of the magnetic anisotropic energy, which is related to the energy of acoustic magnons.Although the variation of the magnon energies is sensitive to strain, the impact of such variations on the phonon energy corrections is negligible.This is clarified by the perturbative analysis in Appendix D, particularly evident in Eq. (D20).Given that the energy difference between the E phonons and acoustic magnons is over 30 meV, a 1 meV variation in the magnon energy leads to a change of only about 3% in the phonon energy correction.Thus, changes in acoustic magnon energy, or the magnetic anisotropic energy, do not significantly influence the phonon energy corrections.

Summary of Cr2O3
In summary, our analysis shows that perturbed phonons possess D 3 symmetry, while bare phonons exhibit D 3d symmetry.The E-type irreps of both D 3 and D 3d are two-dimensional, implying that there is no energy splitting despite the breaking of TRS.However, an E doublet can still be decomposed into two chiral phonons with different chiralities that respect C 3 symmetry.These chiral phonons can be excited by circularpolarized photons with different handedness, and their degeneracy will be lifted if an external magnetic field is present.Each perturbed phonon consists of components from two irreps of bare phonons, yet the level of inter-irrep mixing remains minimal.This characteristic implies that the inversion symmetry within the phonon sector is not strongly broken.Despite this minimal mixing, it is necessary to label perturbed phonons according to the irreps of the unitary subgroup of the magnetic point group.As a result, the infrared and Raman activity properties of these perturbed phonons differ from those of the bare phonons -a phenomenon that is confirmed by recent experimental observations.Magnons belong to the E irrep, and are therefore doubly degenerate.The magnons are mixtures of E g and E u sectors, indicating a strong inversion symmetry breaking for magnons.

C. Perturbed phonons and magnons in 2D systems
In Secs.IV A and IV B, we have investigated chiral phonons in 3D systems, considering both FM and AFM cases.In this section, we shift our focus to 2D systems, specifically a monolayer of FM CrI 3 and a monolayer of AFM VPSe 3 .

Chiral phonons in the monolayer CrI3
Bulk CrI 3 has van der Waals gaps between layers, which allows it to be exfoliated to a 2D single layer while maintaining FM order [68].In this section, we report our calculations of chiral phonons in monolayer CrI 3 , with a focus on comparing the results with those in the bulk cases.
The crystal structure of monolayer CrI 3 is depicted in Fig. 7.The magnetic group P 31m ′ is a Type-III blackwhite group, as was the case for bulk Cr 2 O 3 .The CrI 3 monolayer has higher structural symmetry than that of bulk CrI 3 due to the presence of dihedral mirrors (σ d ) and two-fold rotations about an in-plane axis (C ′ 2 ) that are absent in the bulk.This results in a structural D 3d point group, which is twice the size of the structural group S 6 of bulk CrI 3 .However, the newly added symmetries are all antiunitary in the presence of the FM spin ordering, so that the unitary group is again just S 6 .
Thus, we again label the perturbed phonons at Γ using irreps of S 6 , under which they decompose to 4A g ⊕ 4E g ⊕ 4A u ⊕ 4E u , where one A u mode and two E u modes are acoustic modes.Similar to bulk CrI 3 , the E g and E u irreps in the S 6 point group correspond to two 1D complex irreps rather than a single 2D irrep.As a result, the E g and E u chiral phonons are no longer doubly degenerate and have energy splittings when compared to bare phonons.We refer to the E g and E u phonons as chiral phonons because they individually respect C 3 symmetry, but instead of having identity C 3 eigenvalues, they have complex eigenvalues (ε or ε * ).In Table VIII, we present the energies, angular momenta, and C 3 eigenvalues of the E g and E u chiral phonons, along with the energies of the bare phonons as a reference.Additionally, we include the results for magnons.
Comparing Tables VIII and I, we observe that E g chiral phonons in monolayer CrI 3 have a similar energy splitting to those in bulk CrI 3 , but with larger PAM because the in-plane motion has been enhanced relative to out-ofplane motion.In monolayer CrI 3 , the E u chiral phonons around 10 meV have the largest energy splitting because they are closest in energy to the E u magnon.Overall, the mechanisms for E g and E u chiral phonons in bulk and monolayer CrI 3 are similar despite differences in energy splittings and magnon energies.
In order to discuss the phonons belonging to the A irreps, it is essential to first determine the symmetries of the bare phonons.Similar to our argument in bulk Cr 2 O 3 , since σ d T and C ′ 2 T remain as symmetries, the bare phonons must possess D 3d symmetries instead of S 6 .We have found that the Γ-point bare phonons can be decomposed into 2A 1g ⊕2A 2g ⊕4E g ⊕1A 1u ⊕3A 2u ⊕4E u , wherein one A 2u mode and two E u modes are acoustic.The presence of magnetic moments breaks the σ d and C ′ 2 symmetries of the bare phonons, which leads to interirrep mixing between the bare phonons of different irreps in the perturbed phonons.According to the correlation table of the D 3d point group, the A 1g and A 2g irreps combine to form the A g irrep of S 6 , while the A 1u and A 2u irreps combine to form A u .However, as in bulk Cr 2 O 3 , the inter-irrep mixing in monolayer CrI 3 is found to be relatively weak, as demonstrated by the results presented in Table IX.This allows us to continue to refer to the perturbed phonons as, e.g., A 1g -like based on their pre- dominant character.Nevertheless, it is important to emphasize that the correct labeling of zone-center phonons should consider the irreps of the S 6 group rather than the D 3d group.Table IX contains the energies, angular momenta, and inter-irrep mixing of A g and A u perturbed phonons.Bare phonon energies are also included as a reference for the energy shifts.By comparing Tables II and IX, we observe that the perturbed phonons from the A irreps in both bulk and monolayer CrI 3 exhibit small energy shifts and have a small PAM, indicating a weak effect of the G (pp) matrices in both cases.
Concerning the magnons, we find that the energy of the optical magnon in monolayer CrI 3 is lower than in the bulk.This is attributable to the absence of ferromagnetic inter-layer exchange in the monolayer [75].In contrast, the energy for the acoustic magnon is greater than that in the bulk.As the acoustic magnon energy is related to the magnetic anisotropy energy (MAE), this suggests a greater MAE for monolayer CrI 3 .This observation agrees with findings from another study [76], where the generalized-gradient approximation to the exchangecorrelation functional was used [77].The numerical values for the matrices K (ss) and G (ss) are provided in Appendix C.
In summary, as is the case for the chiral phonons in bulk CrI 3 , the E g and E u chiral phonons in monolayer CrI 3 are no longer doubly degenerate and possess significant angular momentum.The A g and A u perturbed phonons exhibit small energy shifts with respect to the bare phonons and acquire non-zero angular momentum.However, unlike the bulk case, A g and A u perturbed phonons in monolayer CrI 3 exhibit inter-irrep mixing due to the symmetry reduction caused by the presence of magentic order.From the perspective of optical activity, the A 2g -like perturbed phonons are weakly Ramanactive, since they contain some components from the A 1g sector, whereas the A 1u -like perturbed phonons acquire some IR activity due to components from A 2u .

Phonons in monolayer VPSe3
In this section, we investigate the phonons of a VPSe3 monolayer, which is a 2D AFM insulator.It is predicted to have a Néel-type AFM structure [63], with the magnetic moment oriented along the z-direction.The crystal structure of monolayer VPSe 3 is depicted in Fig. 8.The magnetic space group is P 3′ 1m, which is again a Type-III black-white group.The structural point group is D 3d , but the magnetic moments from V atoms break both inversion (i) and two-fold rotational (C ′ 2 ) symmetries, reducing the unitary point group to C 3v .The three-fold rotation (C 3 ) and dihedral mirror (σ d ) remain symmetries of the magnetic group, together with operations iT , C ′ 2 T , and S 6 T .According to the character table for the C 3v point group, the E irrep is a true 2D irrep, i.e., not a complex conjugate pair of 1D irreps.Therefore, perturbed phonons belonging to E irreps must be doubly degenerate, which is confirmed by the numerical results in Table X.The energies of the perturbed phonons are described by energy shifts ∆E with respect to the bare phonon energies (E 0 ).
The VPSe 3 monolayer maintains iT symmetry, similar to Cr 2 O 3 .This symmetry forbids the existence of a nonzero total PAM but permits the presence of nonzero ARPAMs.The magnetic space group of VPSe 3 is P 3′ 1m, with V, P, and Se atoms occupying the 2c, 2e, and 6k Wyckoff positions, respectively.We find that V atoms possess only nonzero L z , Se atoms can only have nonzero in-plane L, and P atoms cannot carry angular momentum.The ARPAMs for the perturbed phonons, traced over the subspace of degenerate doublet for phonons with E irreps, are provided in Table X.It is possible to further decompose the doublets from E irreps into two chi- ral phonons, which individually respect the C 3 symmetry but have different complex eigenvalues.
Building on our earlier discussion, similar to monolayer CrI 3 and Cr 2 O 3 , the D 3d symmetry is preserved by K (pp) and M (pp) in monolayer VPSe 3 , allowing us to use irreps of D 3d to label zone-center bare phonons.These bare phonons decompose to 5E g ⊕ 5E u ⊕ 3A 1g ⊕ 4A 2u ⊕ 1A 1u ⊕ 2A 2g , with one A 2u and two E u modes being acoustic modes.However, due to the presence of magnetic moments on V atoms, i and C ′ 2 symmetries are broken, reducing the symmetry to C 3v and requiring the use of irreps from this point group to label perturbed phonons instead of D 3d .The inter-irrep mixing resulting from the symmetry reduction can be determined from the correlation table, which shows that E g and E u bare phonons form E perturbed phonons, A 1g and A 2u bare phonons form A 1 perturbed phonons, and A 2g and A 1u bare phonons form A 2 perturbed phonons.As in the cases of monolayer CrI 3 and Cr 2 O 3 , the inter-irrep mixing caused by the magnetic order is small.Thus, we still refer to the perturbed phonons as E g -like, etc., but again the correct labeling of those perturbed phonons should involve the irrep of C 3v rather than those of D 3d .
The inter-irrep mixing ρ is defined in the same manner as with Cr 2 O 3 using Eq. ( 29), and the results are presented in Table X.The inter-irrep mixing leads to anomalous Raman/IR activities as well.For example, the A 1g -like and E g -like perturbed phonons, which are strongly Raman-active, now acquire small IR activities, In addition, we have studied the energies of bare and perturbed magnons.Both bare and perturbed magnons are doubly degenerate as they belong to the E irrep.The bare magnon energy is found to be 8.58 meV, while the energy shift of the perturbed magnon is −0.81 µeV, which is more significant than in Cr 2 O 3 .This can be attributed to the larger SOC in VPSe 3 .We further analyzed the symmetry of the perturbed magnon by decomposing it into E g and E u eigenmodes of K (ss) .Our findings reveal that ρ(E u ) = 0.9992, while ρ(E g ) = 0.0405, indicating that the magnon is almost entirely of E u symmetry, similar to the acoustic magnon in Cr 2 O 3 .The numerical values for the matrices K (ss) and G (ss) are provided in Appendix C.
To summarize, we have computed the perturbed phonons of monolayer VPSe3.While bare phonons retain D 3d symmetry, perturbed phonons only possess C 3v symmetry due to the presence of magnetic moments.As in the case of Cr 2 O 3 , the E-irrep in C 3v is a 2D irrep, preserving degeneracy despite the broken TRS.The total PAM of each perturbed phonon is zero due to iT symmetry, but perturbed phonons can possess nonzero ARPAM.Moreover, we have also investigated the energies of bare and perturbed magnons, finding that the energy shift of magnons is greater than in Cr 2 O 3 .Although the magnon is almost entirely of E u symmetry, the E g -E u mixing in the magnon sector is more pronounced than in the phonons.
V. DISCUSSION A. Role of spin-orbit coupling Spin-orbit coupling (SOC) is essential to the physics described above.We identify three matrices, namely K (sp) , G (sp) , and G (pp) , all of which either violate TRS or mediate the interaction between phonons and magnons.
Here we show that in collinear systems such as those considered here, all three of these matrices vanish in the absence of SOC.Since these matrices control the splittings of degenerate phonon modes, it follows that these splittings also vanish without SOC.This is demonstrated here using formal arguments and then confirmed via numerical calculations on bulk CrI 3 .
In the absence of SOC, the presence of global spin rotational symmetry implies that the exchange interactions are of pure Heisenberg form, H = ⟨ij⟩ J ij S i • S j , where the J ij depend on the atomic coordinates.Spin-phonon coupling in noncollinear systems is often described in terms of exchange striction, i.e., the first-order changes of J ij with atomic displacements (see, e.g., Ref. [78]).However, such variations of J ij do not induce any spin canting in an SOC-free collinear system.This follows because the energy is stationary with respect to canting of any spin, since δH = ⟨ij⟩ J ij [S i • δS j + δS i • S j ], and S i and δS j are orthogonal in a collinear spin system.
A complementary point of view comes from noting that the spinor wavefunctions are separable into real spatial wavefunctions with pure spin-up or spin-down character in an SOC-free collinear magnet.This remains true as atoms are displaced, so that there is no induced spin canting.This also explains why G (pp) vanishes.For example, consider three structural configurations, a reference '0' and configurations with displacements δq i and δq j .The Berry phase clearly vanishes, since all inner products are real.
The argument for the vanishing of G (sp) is slightly more subtle.This time one of the displacements, say δq j , is replaced by a spin canting δs j of one spin.The spin system is no longer collinear, but it is still coplanar.In this case the spinors can be represented using Pauli matrices σ 3 and σ 1 to span the plane in which the spins lie.Both of these Pauli matrices are real, so the overall spatial-spinor wavefunctions remain real, and the Berry phase around the loop remains zero.
To confirm these conclusions, we repeated the calculation of the spin splittings in bulk CrI 3 while varying the strength of the SOC used in the first-principles calculations.This adjustment required the introduction of a multiplicative factor in the computational code, followed by recompilation. Figure 9 shows the variation of the phonon energy splittings ∆E with the scaling factor λ for the SOC, where λ = 1 is the physical strength.Clearly the splittings get smaller as λ is reduced and vanish at λ = 0, consistent with our expectations.
We find that K (sp) , G (sp) and G (pp) have linear as well as quadratic terms in λ (with the exception of the K (sp) coupling to acoustic magnons).The quadratic terms are dominant in most cases in Fig. 9, although the linear behavior dominates for the highest-energy E g mode.
We emphasize that these conclusions hold for collinear spin systems.While we have not studied noncollinear systems in this work, we fully expect the present approach to apply in that case.In particular, exchange striction will induce bilinear couplings between atomic displacements and spin cantings as captured by our spinphonon Hession K (sp) , even in the absence of SOC.The role of SOC in the spin-phonon coupling of noncollinear systems may be an interesting avenue for future research.

B. Connection to experiment
Before concluding this paper, we would like to provide some comments on experiments.In materials where chiral phonons exhibit an energy splitting, it is possible to measure this splitting using Raman or IR spectroscopy.Since the splittings for E g chiral phonons are usually small, it may be more feasible to measure splittings of E u chiral phonons using IR spectroscopy.Additionally, circular-polarized photons with different handedness can excite chiral phonons with different chiralities, and observing the peak shift on the spectrum measured using different circular-polarized photons is possible.Another approach is to measure the magnetic moment of the excited chiral phonon using circular-polarized photons, which can be applicable regardless of the presence of energy splitting.Light-induced demagnetization has been observed in many materials, including bulk CrI 3 [79], due to electronic excitation of crystal field levels.However, by adjusting the energy of photons to resonate with a chiral phonon but not to excite electrons, one can selectively probe the chiral phonons because the energy differences between crystal field levels or semicore levels are typically larger than phonon energies.Therefore, electronic excitations can be avoided when probing chiral phonons selectively.

VI. CONCLUSION
Based on a Lagrangian formulation, we developed a theoretical formalism and computational methodology to determine adiabatic dynamics in systems with multiple slow degrees of freedom.Our computational methodology is based on static constrained DFT calculations of Hessians and Berry curvatures with respect to the slow parameters in order to extract the semiclassical dynamics.This constitutes a general and computationally efficient approach that does not require explicit timedependent ab initio calculations as in, e.g., TDDFT, or other specialized capabilities beyond those commonly found in, or easily added to, widely available DFT codes.We demonstrate the utility of this methodology by applying it to phonons and spins (magnons) in magnetic insulators, as the dynamics of each is generally on the same energy scale.This represents a more systematic way of treating spin-phonon dynamics compared to conventional approaches which rely on building spin models and parameterizing their dependence on atomic displacements.
Specifically, we showed that the inclusion of the Hessians and Berry curvatures involving spin and phonon DOF are required to accurately describe the oftenneglected effects of TRS breaking on phonon modes at the zone center.Results for four case-study materials were presented, covering both FM and AFM ordering, as well as 2D and 3D materials.FM CrI 3 (both 3D bulk and 2D monolayer) exhibits energy splittings of its Etype modes (which are doubly degenerate under TRS) at the zone center.This results in chiral phonons exhibiting circular atomic motion leading to a significant phonon angular momentum.In AFM Cr 2 O 3 (3D bulk) and VPSe 3 (2D monolayer), the E-type modes remain doubly degenerate, but the inversion symmetry breaking from the magnetic order results in modes with mixed E g and E u character.The angular momentum of a phonon in these materials is zero in total, but can exhibit welldefined nonzero atom-resolved contributions.In addition, the mixing of E g and E u modes implies that the Raman modes acquire some small IR character and vice versa, a feature that is open to experimental confirmation.
This work opens up various directions for future study.The general methodology will be useful in any material with interesting dynamics involving multiple slow DOF.Specifically in the area of spin-phonon dynam-ics, it allows the search for materials with large splittings of chiral modes, due for example to large SOC or approximate degeneracies between phonon and magnon branches.It also motivates the exploration of dynamics of materials with more complicated magnetic structures, such as noncollinear or so-called altermagnetic [80,81] orders.Finally, we focus here on zone-center modes, but the methodology shows promise for being generalized for computing the coupled spin-phonon dynamics across the Brillouin zone, which is relevant for the thermal Hall effect.
Overall, we expect that the developments in this work will allow efficient and accurate calculations of generalized adiabatic dynamics, and thus the exploration of much resulting novel physical phenomena, in spinphonon coupled systems and beyond.in order to make γ real.Therefore, we assume γ to be real, resulting in ⟨s ± |K (sp) We can expand |v⟩ and |s⟩ in the bases of |v ± ⟩ and |s ± ⟩ as where x ± and s ± are coefficients representing the magnitudes of the contributions from the corresponding basis vectors.In the minimal spin-phonon model, it is assumed that the mixing between different E g doublets is negligible.However, in reality, this mixing, which is mediated by magnons, plays a crucial role in explaining why the angular momenta of the two circular polarized modes do not cancel out exactly.A more detailed analysis of this mixing based on the perturbation approach can be found in Appendix D 1.
By multiplying ⟨v ± | to Eq. (B2) and substituting Eqs.(B4), (B10), and (B8), we obtain Similarly, by multiplying ⟨s ± | to Eq. (B3) and substituting Eqs.(B5), (B7), and (B8), we get Thus, we have successfully derived Eq. ( 6) in Ref. [33].Although Eq. (B1) is derived from an adiabatic Lagrangian formalism, we aim to demonstrate that it can also emerge from the undamped Landau-Lifshitz equation [52].In its original form, the Landau-Lifshitz equation is presented as where H is the Hamiltonian, dH/d ⃗ M is the effective magnetic field, and γ represents the gyromagnetic ratio, which is the ratio of the magnetic moment ⃗ M to its corresponding angular momentum ⃗ S. As the dynamic variable in our study is the unit vector of spin, ⃗ s = ⃗ S/S, we can recast Eq. (B14)) in terms of ⃗ s as Referring back to Eq. (B6), it is noteworthy that −Sϵ is equivalent to the G (ss) .Since the Hamiltonian has been expanded up to the quadratic order, the derivative dH/d⃗ s can be expressed as where we have chosen to use the bra-ket notation for ease of representation.This gives rise to which aligns with the form presented in Eq. (B1).Thus, we have successfully shown the equivalence of the EOM for the spin-phonon model presented in this Appendix with the undamped Landau-Lifshitz equation.
Appendix C: Spin-spin Hessians (K (ss) ) and Berry curvatures (Ω (ss) ) for all four materials In this appendix, we present the spin-spin Hessians and Berry curvatures for all four materials investigated in our study.
For bulk CrI 3 , both the spin-spin Hessian (K (ss) ) and Berry curvature (Ω (ss) ) matrices have the form where both a and b are 2 × 2 matrices that are restricted to be of the form and similarly for b, as a result of the three-fold symmetry.
In Table XI we provide the symmetric and antisymmetric parts of a and b for both K (ss) and Ω (ss) of bulk CrI 3 .It is noteworthy that the eigenvalues of Ω (ss) are 1.499 and 1.566 for the E g and E u modes, respectively.These values represent the effective spins of the respective magnon modes.They are are close to 3/2 as expected from the nominal spin of the Cr 3+ ion, but they do differ slightly, especially for the optical E u magnon.Instead, the deviation is very small for the acoustic E g magnon.These deviations were not captured by the minimal spin-phonon model proposed in Ref. [33].
The symmetry-allowed matrix elements for monolayer CrI 3 exhibit a similar structure to those of bulk CrI 3 , and they can also be decomposed following the same procedure as described in Eqs.(C1) and (C2).The numerical values of these matrix elements are provided in Table XI.  ss) and Ω (ss) for materials understudy (ML = monolayer).Symmetric ('sym') and antisymmetric ('asym') parts are defined in Eq. (C2).For Ω (ss) of Cr2O3, the values of csym = −6.211× 10 −5 and dsym = −2.409× 10 −4 appear rounded to zero.For bulk Cr 2 O 3 , the non-zero matrix elements of K (ss) are given by Here, a, b, c, and d can also be decomposed in the same manner as in Eq. (C2).The numerical results for these matrix elements are provided in Table XI.
For monolayer VPSe3, the matrix elements of both K (ss) and Ω (ss) have the form where a and b are 2 × 2 matrices.The values of these matrix elements can be found in Table XI.
From Table XI we can find that, for the systems with at most two spins, the asymmetric (symmetric) parts of K (ss) (G (ss) ) all vanish.However, the symmetry is more complicated in Cr 2 O 3 , where some of the b, c, and d components in Eqs.(C3) and (C4) have both symmetric and asymmetric components.In the main text, we have directly solved the equation of motion (EOM) of the present approach (Eq.( 10)).However, employing perturbation theory to solve the EOM can offer valuable physical insights into the influence of Hessian matrices and Berry curvatures on the results.In this section, we will present a perturbation treatment for K (sp) , G (ps) , and G (pp) individually.Furthermore, we will provide numerical results from the perturbation treatment of K (sp) and compare them with the predictions of the spin-phonon model.
We begin by reformulating Eq. ( 10) in matrix form as where the matrices M , K, and G are defined in Eq. (11).At variance with the main text, in this Appendix we let i, j label mixed modes, while m, n and µ, ν label bare phonon and bare magnon modes respectively.We utilize the bare phonon and magnon states |u µ ⟩ as the basis functions for the perturbation theory, where the superscript '(0)' denotes the zeroth order in the perturbation expansion.
The EOM for unperturbed phonons can be expressed in matrix form as (K (pp) − ξ n M (pp) )|u (0)  n ⟩ = 0 (D2) where for convenience we have introduced the squared frequency ξ n = ω 2 n , and for the magnons as where ω µ may be positive or negative even though solutions with negative energy are not physically observable.The fact that ⟨u As discussed in the main text, phonons and magnons belonging to the E irreducible representations can be doubly degenerate.To address this, we need to employ degenerate perturbation theory.However, we can simplify the analysis by considering the '+' and '−' sectors separately.Phonons and magnons from different sectors do not mix, allowing us to work with these bases throughout this section.
1. Perturbation treatment of K (sp)   We begin by examining the perturbation of K (sp) alone, neglecting G (pp) and G (sp) .We replace K (sp) with λK (sp) , where λ serves as our perturbation parameter, and expand Eq. (D1) up to second order in λ as K = K (0) + λK (1) , K (0) = K (pp) 0 0 K (ss) , K (1) = 0 K (ps) K (sp)  0 , i + λ 2 ω (2) i + . . ., ξ i = ξ (D20) This is a second major result.The summation over µ in Eqs.(D16) and (D20) runs over all solutions of Eq. (D3), including those of negative energy.Thus, the unphysical negative-energy states can still contribute to the perturbation of the phonons.It is worth noting that Eq. (D20) explains why phonons from the '+' sector exhibit larger energy splitting.Since positive-energy magnons also belong to the '+' sector, the interactions between magnons and phonons are stronger in the '+' sector due to the smaller energy denominator.Furthermore, we can obtain |u If i labels a phonon and j labels a magnon, both the second and third terms in Eq. (D21) are zero, implying that the first term is also zero.It follows that |q m ⟩ , (D22) which represents the effect of the magnon-mediated phonon-phonon interaction at second order.It provides an explanation for the non-canceling angular momentum observed in chiral phonon pairs, as shown in Table I and VIII.

b. Perturbation of magnons
We now turn to the corresponding perturbation treatment of the magnons.First, we replace |q i ⟩ and |q j ⟩ in Eq. (D11) with |s µ ⟩ and |u n ⟩, respectively.Then, we multiply Eq. (D11) by |u    From Eq. (D25), we can determine the coefficients d nµ .Therefore, the first-order perturbation of the magnons is given by |s (1)  µ ⟩ = (D27) The second order perturbation of the magnon states can be obtained in a similar way as was done for the phonons.We replace |q i ⟩ and |q j ⟩ with |s µ ⟩ and |s ν ⟩ and substitute Eq. (D26) into Eq.(D21).This substitution yields |s (2)  µ ⟩ = ν̸ =µ n σ ν ⟨s (D28) We are now prepared to present numerical results obtained using the perturbation approach.The neglect of G (sp) and G (pp) corresponds to the spin-phonon model presented in Appendix B, which can be solved exactly.

FIG. 1 .
FIG. 1. Visualization of crystal structure and magnetic moments of bulk CrI3 unit cell.Cr atoms are depicted in blue, while top-layer and bottom-layer I atoms are in bright and dark magenta, respectively.The black vectors indicate the magnetic moments, which are oriented along the z-direction.

FIG. 2 .
FIG. 2. Visualization of the real and imaginary components of the chiral phonons E(1) g and E (2) g .Cr atoms are depicted in blue, while top-layer and bottom-layer I atoms are in bright and dark magenta, respectively.The green vectors represent the real part of the phonon displacement, and the blue vectors represent the imaginary part.The real part of E (1) g and E(2) g are nearly identical, while the imaginary parts are in opposite directions.

( 1 )
g .Cr atoms are depicted in blue, while top-layer and bottom-layer I atoms are in bright and dark magenta, respectively.The green vectors represent the real part of the phonon displacement, and the blue vectors represent the imaginary part, which are amplified 200 times.

FIG. 4 .
FIG. 4. A visualization of the crystal structure and magnetic moments of the unit cell of bulk Cr2O3.The Cr atoms are depicted in blue, while the top-layer and bottom-layer O atoms are shown in bright and dark red, respectively.The black vectors indicate the magnetic moments, which are oriented along the z-direction.

( 1 )
g perturbed phonon around 36 meV and the A (2) 2g perturbed phonon near 57 meV, where the L z for Cr atoms and the in-plane angular momentum for O atoms are clearly visible.

FIG. 5 .
FIG. 5.Atom-resolved phonon angular momentum (ARPAM) for two selected modes in bulk Cr2O3: E (1) g -like mode near 36 meV and A (2) 2g -like mode near 57 meV.The ARPAM for the Eg-like mode is traced over two degenerate modes.The amplitude of ARPAM in this figure is amplified by a factor of 1000 compared to Fig. 6.

FIG. 6 .
FIG. 6.The real and imaginary components, as well as atomresolved phonon angular momentum (ARPAM), of the chiral phonon E (1)+ g around 36 meV in Cr2O3.The Cr atoms are shown in blue, while the top-layer and bottom-layer O atoms are depicted in bright and dark red, respectively.The green vectors indicate the real part of the phonon displacement, the blue vectors indicate the imaginary part, and the brown vectors represent the ARPAM.The real parts of E (1)− g are almost identical to E (1)+ g , while its imaginary parts and ARPAM are almost opposite to those of E (1)+ g .

FIG. 7 .
FIG. 7. Side (a) and top (b) views of the crystal structure of monolayer CrI3, with a dashed line indicating the unit cell.Cr atoms are depicted in blue, while top-layer and bottomlayer I atoms are in bright and dark magenta, respectively.The black vectors indicate the magnetic moments, which are oriented along the z-direction.

FIG. 8 .
FIG. 8. Side (a) and top (b) views of the crystal structure of monolayer VPSe3, with a dashed line indicating the unit cell.Bright and dark shades denote V atoms with spin up and down, respectively.Purple atoms represent P atoms, while top-layer and bottom-layer Se atoms are in bright and dark green, respectively.The magnetic moments, oriented along the z-direction, are denoted by black vectors.

FIG. 9 .
FIG.9.Energy splittings (∆E) for the Eg chiral phonon modes (a) and Eu chiral phonon modes (b) in bulk CrI3, plotted as a function of the spin-orbit coupling strength scaling factor (λ).The ∆E represents the energy difference between the '+' mode and the '−' mode near that energy.
a b c d b T a T d c T c T d T a b T d T c b a T non-zero matrix elements of G (ss) are G (ss) b T a T d c T −c T −d T a b T −d T −c −b a T Appendix D: Perturbation treatment of phonon-magnon dynamics n ⟩ = 0 for ξ m ̸ = ξ n and that ⟨s (0) µ |G (ss) |s (0) ν ⟩ = 0 for ω µ ̸ = ω ν allows us to adopt the normalization conditions⟨u (0) m |M (pp) |u (0) n ⟩ = δ mn , (D4)⟨s (0) µ |iG (ss) |s (0) ν ⟩ = δ µν σ µ ,(D5)where σ µ takes the values ∓1 for ω µ > 0 and ω µ < 0 respectively.The solution to Eq. (D1) is a composite vector that includes both the phonon and magnon sectors, living in the space spanned by basis vectors |u equations into Eq.(D1) and expand in orders of λ, making use of the relations ξ

( 2 )
i ⟩ can only possess phonon character.Replacing |q i ⟩ and |q j ⟩ by |u n ⟩ and |u m ⟩ and substituting Eq. (D16) into Eq.

TABLE II .
Computed properties of zone-center Ag and Au phonons in bulk CrI3.Energy shifts ∆E are defined relative to the bare phonon energies E0.Lz is the phonon angular momentum of Eq. (15).
1. Importance of G(pp)and G (ps) in calculating energies

TABLE III .
Energies of the bare magnons E0, energy shifts of perturbed magnon modes comparing to bare magnons ∆E, and the corresponding C3 eigenvalues χ(C3) for both Eg and Eu modes in bulk CrI3.

TABLE IV .
Bare phonon irreps, energies, atom-resolved phonon angular momentum (ARPAM) of Cr atoms in zdirection (Cr Lz), ARPAM of O atoms in x-direction (O Lx), and inter-irrep mixing (ρ) of perturbed phonons in bulk Cr2O3.The bare phonon energies E0 are included as a reference for the energy shifts ∆E.E-irrep (Eg and Eu) phonons are doubly degenerate.Note that '0' denotes an entry that is precisely zero by symmetry, while '0.000' signifies a non-zero entry that has been rounded to zero.The energy shifts ∆E for A-irrep perturbed phonons are on the order of 10 −4 µeV, and Lx for the O atom of the A

TABLE V .
Raman and infrared (IR) activities for bare and perturbed phonons.For perturbed phonons, the Irrep * column indicates the irrep of perturbed phonons labeled by irreps of D 3d , which is made possible by the weak inter-irrep mixing.In cases where a perturbed phonon exhibits both Raman and IR activity, the minor activity is presented within parentheses and is subordinate to the major activity.

TABLE VI .
Possible configurations of atom-resolved phonon angular momentum (ARPAM) for bulk Cr2O3 according to the magnetic space group R 3′ c ′ .Cr atoms can only have outof-plane angular momentum Lz or −Lz, while O atoms can only have in-plane angular momentum and must respect C3 symmetry.

TABLE VII .
Energies of the bare magnons (E0) and perturbed magnons (E) in bulk Cr2O3, as well as their Eg and Eu components defined in Eq. (29).

TABLE VIII .
Energies, z-direction total PAM (Lz), and C3 eigenvalue (χ(C3)) of Eg and Eu chiral phonons and magnons in monolayer CrI3 calculated using the full model.The bare phonon and magnon energies E0 are included as a reference for the energy shift ∆E.

TABLE X .
Bare phonon irreps, energies, atom-resolved phonon angular momentum (ARPAM) of V atoms in the z-direction (V Lz), ARPAM of Se atoms in the y-direction (Se Ly), and inter-irrep mixing (ρ) of perturbed phonons in monolayer VPSe3.Bare phonon energies E0 are included as a reference for energy shifts ∆E.The E-irrep phonons are doubly degenerate.The energy shifts ∆E for A-irrep perturbed phonons are on the order of 10 −5 µeV.
α and γ are confined to x and y Cartesian directions.By reintroducing the matrix form and transferring the ϵ tensor to the left-hand side, we arrive at

TABLE XI .
Matrix elements K

TABLE XII .
Comparison of perturbation approach and exact solution in the spin-phonon model.The perturbation approach provides second-order perturbations to phonon and magnon energies E νn .The exact solutions for phonon and magnon energy shifts ∆E = ℏ∆ω and the spin component of phonon-like solutions cνn are included as benchmarks.The bare phonon and magnon energies E