Giant effective magnetic fields from optically driven chiral phonons in 4$f$ paramagnets

We present a mechanism by which optically driven chiral phonon modes in rare-earth trihalides generate giant effective magnetic fields acting on the paramagnetic $4f$ spins. With cerium trichloride (CeCl$_3$) as our example system, we calculate the coherent phonon dynamics in response to the excitation by an ultrashort terahertz pulse using a combination of phenomenological modeling and first-principles calculations. We find that effective magnetic fields of over 100 tesla can possibly be generated that polarize the spins for experimentally accessible pulse energies. The direction of induced magnetization can be reversed by changing the handedness of circular polarization of the laser pulse. The underlying process is a phonon analog of the inverse Faraday effect in optics that has been described recently, and which enables novel ways of achieving control over and switching of magnetic order at terahertz frequencies.

We present a mechanism by which optically driven chiral phonon modes in rare-earth trihalides generate giant effective magnetic fields acting on the paramagnetic 4f spins. With cerium trichloride (CeCl3) as our example system, we calculate the coherent phonon dynamics in response to the excitation by an ultrashort terahertz pulse using a combination of phenomenological modeling and first-principles calculations. We find that effective magnetic fields of over 100 tesla can possibly be generated that polarize the spins for experimentally accessible pulse energies. The direction of induced magnetization can be reversed by changing the handedness of circular polarization of the laser pulse. The underlying process is a phonon analog of the inverse Faraday effect in optics that has been described recently, and which enables novel ways of achieving control over and switching of magnetic order at terahertz frequencies.
Particularly interesting are circularly polarized, or chiral, phonons, where the ions in a solid move on closed elliptical or circular orbits, generating angular momentum that can be transferred to the spins [29][30][31][32][33][34]. Previously, nondegenerate chiral phonons have been described at the Brillouin-zone edges of materials with hexagonal symmetries [35][36][37][38][39][40], but their direct excitation with an ultrashort laser pulse is prohibited due to the large momentum mismatch between photons and phonons. In contrast, degenerate chiral phonons consist of superpositions of two orthogonal components of doubly or triply degenerate phonon modes that can be found at the Brillouin-zone center of materials with uniaxial or cubic symmetries and can therefore be resonantly excited with light. In recent years, a number of studies have shown or predicted effective magnetic fields arising from coherent chiral phonon driving that reach the milli tesla range [8, 9, 12-14, 41, 42].
Here, we propose that optically driven chiral phonons in rare-earth trihalides produce giant effective magnetic fields that exceed those previously seen by several orders of magnitude. We predict, at the example of CeCl 3 , that effective magnetic fields of over 100 tesla should be achievable, which polarize the paramagnetically disordered spins, for laser energies well within the damage threshold of the crystal. The mechanism allows for bidirectional control of the induced magnetization through phonon chirality that in turn can be controlled by the polarization of the laser pulse.

II. PROPERTIES OF CERIUM TRICHLORIDE
Rare-earth trihalides are a class of 4f paramagnets with formula unit RH 3 . CeCl 3 (R = Ce, H = Cl) is a representative of this class of materials that crystallizes in the hexagonal P 6 3 /m structure [43] with an electronic band gap of 4.2 eV [44]. Its Ce(4f 1 ) valence-electron configuration remains paramagnetic for all temperature ranges relevant here, as spin ordering only occurs at very   [45]. We chose CeCl 3 as our model system, because the primitive unit cell consists of only 8 atoms ( Fig. 1(a)), resulting in a small number of 21 optical phonon modes characterized by the irreducible representations 2A g + 1A u + 2B g + 2B u + 1E 1g + 3E 2g + 2E 1u + 1E 2u in its 6/m point group. Early Raman studies have shown that the polarization of the 4f electrons in an external magnetic field leads to a splitting of the doubly degenerate E 1g and E 2g phonon modes into leftand right-handed circular polarization [46,47], therefore obtaining chirality, see Fig. 1(b). It has been suggested that also the infrared-active E 1u phonon modes split in the same way [48], yet no experimental infrared spectroscopy measurements had been performed at that time.
The infrared-active E 1u modes map into the same E representation at the local6 symmetry of the cerium ions as the Raman-active E 2g modes, for which phonon splittings have been measured, and should therefore have the same effect on the paramagnetic spins. Infrared-active phonon modes possess an electric dipole moment and can therefore be resonantly excited by the electric field component of a laser pulse to yield large vibrational amplitudes. We will explore in this work how optically driven chiral E 1u phonons act on the spins through the inverse of the spinphonon coupling.

III. SPIN-PHONON COUPLING AND COHERENT PHONON DYNAMICS
We begin by reviewing the theory of spin-phonon coupling in 4f paramagnets. Motions of the ions along the eigenvectors of doubly degenerate phonon modes modify the crystal electric field around the paramagnetic ions and induce virtual transitions between the ground-state energy levels and higher-lying states, see Fig. 2. The spin states in 4f paramagnets are close to those of the free ions and the total angular momentum (isospin), J, is a good quantum number. In CeCl 3 , the lowest energy level has J = 5/2, which splits into three Kramers doublets, of which m J = ±5/2 is the ground state. the interaction of chiral phonons with the isospin can be written as an effective "spin-orbit" type Hamiltonian [29,[49][50][51][52][53][54] where K is the coupling coefficient and L = Q ×Q is the phonon angular momentum, with Q = (Q a , Q b , 0) containing the normal mode coordinates of the two orthogonal components of a doubly degenerate phonon mode, Q a and Q b , in the ab plane of the crystal. m is the magnetic moment per unit cell, which, for components perpendicular to the ab phonon polarizations, is given by where g J is the Landé factor, e z is a unit vector along the c axis of the crystal and n ±J are Fermi-Dirac distributions describing the occupation of the groundstate doublet. The theoretical value of the prefactor in Eq. (2) for the m J = ±5/2 ground-state doublet, g ±5/2 = g J J(J + 1) = 2.54, is reasonably close to the experimental value of 2.02 [55], showing that most of the orbital angular momentum is unquenched. Please see the Appendix for detailed derivations.
We now look at the influence of the interaction on the phonons. The spin-phonon coupling modifies the off-diagonal terms of the dynamical matrix as D(m) = [9,13,46,47,[56][57][58][59][60][61]. As a result, the frequencies of right-and left-handed circular polarizations, Ω ± , of the doubly degenerate phonon mode split, where Ω 0 is the eigenfrequency of the doubly degenerate phonon mode. Without an external magnetic field, the energy levels of the ground-state doublet are degenerate, there is no net magnetic moment per unit cell, and the phonon frequencies in Eq. (3) remain degenerate. Applying a magnetic field, B c, to the paramagnet splits the ground-state doublet, and, using Eq. (2), induces a splitting of the form The prefactor in Eq. (4) directly corresponds to the saturation splitting, ∆Ω s = 4Kg ±5/2 µ B , and we can reciprocally extract the spin-phonon coupling from experimentally measured phonon splittings, K = ∆Ω s /(4g ±5/2 µ B ).
We now look at the inverse effect the interaction has on the magnetization, when phonons of one type of chirality are driven with an ultrashort laser pulse. The phonon angular momentum acts as an effective magnetic field, B,

0.
1. Phenomenologically, this type of interaction has recently been described as a phonon analog of the inverse Faraday effect in optics [13], which is known to induce magnetizations in paramagnets [62][63][64][65]. A first experiment demonstrating this effect with elliptically polarized phonons has induced spin and magnetization dynamics in a complex oxide in recent years [8]. If the spins reacted instantaneously to the effective magnetic field, the magnetization could be described statically as in Eq. (4). Experiments on the optical inverse Faraday effect have shown that this static limit of spin response holds for driving pulses on the order of nanoseconds [62,63]. For femtosecond pulse durations, additional diamagnetic effects from the unquenching of electronic orbital moments come into play [64][65][66][67], which cannot be described in the thermodynamic limit. Coherently driven phonons evolve over several picoseconds, which is also the timescale that spins and phonons have been shown to equilibrate through effective magnetic fields [11]. We therefore apply a rate-equation model to describe the dynamics of the spin population of the m J = ±5/2 ground-state doublet, n ±J [68,69], where the decay rates of spins in the respective states are described by , and γ 0 is the decay rate for zero level splitting, γ ±J (∆E → 0) = γ 0 . An ultrashort terahertz pulse can resonantly excite infrared-active phonons into a coherent quantum state, which allows us to treat the normal mode coordinate, Q, as semi-classical field amplitude [70][71][72][73][74][75]. We obtain Q by solving the equation of motion Here, κ is the linewidth of the phonon mode, Ω 0 is its eigenfrequency, and Z = m Z * m q m / √ M m is its mode effective charge, where Z * m is the Born effective charge tensor, q m is the eigenvector, and M m is the atomic mass of ion m. The sum runs over all ions in the unit cell. We model the circularly polarized terahertz pulse as ) cos(ω 0 t), E 0 is the peak electric field, ω 0 is the center frequency, and τ is the full width at half maximum duration of the pulse. Here, the two perpendicular components of the doubly degenerate phonon mode are excited with a quarter-period difference, resulting in circular polarization and therefore chirality. As light couples to phonon modes close to the center of the Brillouin zone, we may neglect any wavevector dependence in Eq. (7).

IV. COMPUTATIONAL DETAILS
We calculate the phonon eigenfrequencies and eigenvectors, and Born effective charges from first principles, using the density functional perturbation theory formalism [76,77] as implemented in the Vienna ab-initio simulation package (vasp) [78,79] and the frozen-phonon method as implemented in the phonopy package [80]. We use the VASP projector augmented wave (PAW) pseudopotentials with valence electron configurations Ce (6s 2 5s 2 5p 6 5d 1 4f 1 ) and Cl (3p 5 3s 2 ) and converge the Hellmann-Feynman forces to 25 µeV/Å. For the 8-atom unit cell, we use a plane-wave energy cut-off of 600 eV, and a 4×4×7 gamma-centered k-point mesh to sample the Brillouin zone. For the exchange-correlation functional, we choose the Perdew-Burke-Ernzerhof revised for solids (PBEsol) form of the generalized gradient approximation (GGA) [81]. We perform nonmagnetic calculations to obtain the structural and dynamical properties of CeCl 3 . A fully ab-initio treatment of paramagnetism in CeCl 3 would require supercell calculations and a description of 4f electron magnetism, which would make the computation of dynamical properties intractable for our purposes. Within the nonmagnetic treatment, the lattice constants of our fully relaxed hexagonal structure (space group P 6 3 /m, point group 6/m) of a = 4.21Å and c = 7.38Å with a unit-cell volume of V c = 199Å 3 agree reasonably well with experimental values [43]. Furthermore, our calculated phonon eigenfrequencies match the experimental values reasonably well [47,48], with a maximum deviation of ∼10%. Crystal structures are visualized using vesta [82].

V. PHONON-INDUCED EFFECTIVE MAGNETIC FIELDS AND MAGNETIZATIONS
We extract the magnitude of the spin-phonon coupling from experimental data of the phonon-frequency splitting according to Eq. (4), K = ∆Ω s /(4g ±5/2 µ B ). In rareearth trihalides, splittings of the Raman-active modes range between 0.3 THz and 0.75 THz [46,47]. Because the infrared-active modes change the local symmetry of the magnetic cerium ion in the same way, we expect a similar strength of the spin-phonon coupling as for the Raman-active modes and use an average of the experimentally found values of ∆Ω s /(2π) = 0.5 THz. This splitting is several orders of magnitude larger than the one induced by the magnetic moments of phonons in the phonon Zeeman effect [9,12,83,84], which we can therefore neglect here. Note that there are further microscopic origins of phonon angular momentum and magnetic moments, such as topological band features in semimetals [85][86][87][88][89][90] or thermal gradients [91,92], which however do not play a role here either.
In the following, we evaluate the effective magnetic fields produced by the two doubly degenerate infraredactive E 1u modes in CeCl 3 with eigenfrequencies of 5.9 and 4.8 THz. We find the mode effective charges of these modes to be 0.24e and 0.66e, respectively, where e is the elementary charge. For the phonon linewidth, κ, we assume a phenomenological value of 5% of the phonon frequency that matches those typically found in rare-earth trihalides [46,47]. Fig. 3 shows the coherent phonon dynamics following the excitation by a circularly polarized terahertz pulse with a duration of τ = 350 fs and a fluence of 10 mJ/cm 2 , as described by Eq. (7). The fluence F is connected to the peak electric field and the duration of the pulse through F = τ / √ 8 ln 2c 0 0 π/2E 2 0 , where c 0 and 0 are the speed of light and the vacuum permittivity. The center frequency, ω 0 , is chosen to be resonant with the eigenfrequencies of the respective phonon modes. In Fig. 3(a), we show the evolution of the phonon amplitudes Q a according to Eq. (7). The phases of the Q b components are shifted by a quarter period, respectively. The maximum amplitude of the E 1u (5.9) mode of Q a = 0.33Å √ amu, where amu denotes the atomic mass unit, is roughly three times smaller than that of the E 1u (4.8) mode of Q a = 1.1Å √ amu due to the smaller mode effective charge and higher phonon frequency. In Fig. 3(b), we show the evolutions of the effective magnetic fields produced by the two chiral phonon modes according to Eq. (5). We obtain a maximum effective magnetic field of B = 2.9 T for the E 1u (5.9) mode and 27 T for the E 1u (4.8) mode. This order-of-magnitude difference  comes from the quadratic scaling of the effective magnetic field with the phonon amplitudes. The direction of the effective magnetic field is determined by the handedness of the phonon chirality, which can straightforwardly be controlled by the handedness of circular polarization of the pulse. We now vary the strength of the excitation. We show the maximum amplitudes of the effective magnetic fields for a range of experimentally accessible fluences of the terahertz pulse [93] in Fig. 3(c), where we fix the pulse duration at τ = 350 fs. The effective magnetic fields depend linearly on the fluence and reach 11.4 T for the E 1u (5.9) mode and 107 T for the E 1u (4.8) mode at a fluence of 40 mJ/cm 2 . In order to ensure experimental feasibility, we evaluate the atomic displacements along the eigenvectors of the phonon modes. The Lindemann stability criterion predicts melting of the crystal lattice when the root mean square displacements reach between 10% and 20% of the interatomic distance [94]. We extract the maximum root mean square displacements as d = max n |d n / √ 2|, where d n = q n Q a (t)/ √ M n is the displacement of ion n. Even at fluences of 40 mJ/cm 2 , the largest root mean square displacements of the chloride ions reach only 1.3% of the interatomic distance of 2.97Å for the E 1u (5.9) mode and 3.8% for the E 1u (4.8) mode, well below the vibrational damage threshold. Note that other effects may occur, e.g. Zener tunneling, that are not accounted for here. At these high fields, nonlinear couplings between coherently excited infrared-active modes and other vibrational degrees of freedoms come into play [73,95]. These modes do not contribute directly to the spin-phonon coupling however, and we therefore neglect the effect of nonlinear phonon-phonon coupling in this context. Furthermore, the centrosymmetry of CeCl 3 prevents nonlinear optical effects, such as secondharmonic generation, to occur at high fluences.
Next, we look at the magnetization, M = m/V c , that can be induced in CeCl 3 according to Eqs. (2) and (6). In Fig. 4, we show the evolution of the magnetization in response to the effective magnetic field generated by the E 1u (4.8) mode when excited with a resonant terahertz pulse with a duration of 350 fs and a fluence of 10 mJ/cm 2 , as well as the dependence of the magnetization on the fluence of the laser pulse. In Figs. 4(a) and 4(b), we vary the decay rate, γ 0 , while keeping the temperature fixed at 4 K, and in Figs. 4(c) and 4(d) we vary the temperature, while keeping γ 0 = 0.1 THz. For fast decay rates and at low temperatures, even small fluences of <10 mJ/cm 2 are sufficient to fully polarize the spins of the material, yielding a transient saturation magnetization of M = 4 µ B /V c . The slower the decay rate and the higher the temperature, the higher the fluence of the laser pulse has to be in order to induce a significant magnetization. The influence of the decay rate on the achievable magnetization is hereby much larger than that of temperature. The slowest decay rate of γ 0 = 0.1 MHz that we look at here corresponds to the nanosecond timescale, on which the thermodynamic picture of spin polarization is known to hold [62,63]. Therefore, the corresponding values of induced magnetization for γ 0 = 0.1 MHz in Fig. 4(b) can be regarded as lower boundary.

VI. DISCUSSION
Our predictions can be experimentally realized in state-of-the-art tabletop setups that provide terahertz pulses in the required frequency range [93,96], where the phonon-induced magnetization of the material can be probed by Faraday rotation measurements. Tuning the frequency of the terahertz pulse in and out of resonance with the phonon modes can distinguish a possible contribution of the optical inverse Faraday effect to the magnetization from the phonon-induced mechanism. The effective magnetic fields calculated here reach 11.4 T for the E 1u (5.9) mode and 107 T for the E 1u (4.8) mode for terahertz pulses with a fluence of 40 mJ/cm 2 . We further predict that the subsequent spin polarization of the paramagnetic cerium ions can possibly reach full saturation at low temperatures. These magnetic fields and moments are several orders of magnitude larger than those demonstrated in previous experiments and calculations of the phonon inverse Faraday effect, where fields can be found around micro to milli tesla and magnetic moments around the order of nuclear magnetons [8,9,13,41]. Furthermore, the fields we predict are orders of magnitude larger than well-established experiments on the optical inverse Faraday effect, where fields of fractions of tesla have been reported for comparable fluences in the visible spectrum for different insulating and semiconducting systems [64,97]. We display a detailed comparison in Table I in the Appendix. A direct quantitative comparison of the optical and phonon inverse Faraday effects remains difficult, as no practical and general ab-initio theory of the mechanisms exists to date, and phonon-pumping experiments have only become feasible in very recent years. This comparison is further complicated by the breakdown of the effective magnetic field picture for pulse durations on the order of tens of femtoseconds [64][65][66][67]. In the future, explicit calculations of spin-phonon decay rates [98] and coupling strengths will therefore be necessary to further quantify the timescale and magnitude of the optical and phonon inverse Faraday effects and to make predictions for a broad range of materials. First steps have been made over the course of the last decade [66,67,[99][100][101], and first quantitative results for the optical inverse Faraday effect have been achieved in metals [102,103].
While we have chosen CeCl 3 as our model system, the mechanism described here should be general to the entire class of rare-earth trihalides [46,47] and possibly to 4f magnets in general, as similar magnitudes of the spin-phonon coupling have been found in ferromagnetic LiTbF 4 [104] and materials exhibiting the phonon Hall effect, such as paramagnetic Tb 3 Ga 5 O 12 [29,52]. A future question to answer is whether spin-phonon couplings in 3d magnets can reach similar magnitudes to those in 4f magnets. Potential giant phonon-induced effective magnetic fields in the paramagnetic phases of 3d magnets would directly impact a large variety of materials that are already being used in magnetoelectronic technologies [105,106]. Another future question is whether coherent chiral phonon excitation could stabilize an ordered spin phase above the equilibrium ordering temperature in paramagnets, similar to phonon-and light-induced superconductivity above the equilibrium critical temperature in superconductors [107][108][109][110].

APPENDIX: PHONON-FREQUENCY SPLITTING AND MAGNETIZATION IN 4f PARAMAGNETS
In this section, we provide detailed derivations of the equations used in the Spin-phonon coupling and coherent phonon dynamics section of the main text. Motions of the ions along the eigenvectors of doubly degenerate phonon modes modify the crystal electric field (CEF) around the paramagnetic ions and induce virtual transitions between the ground-state energy levels and higherlying CEF states, see Fig. (2) in the main text. The spin states of rare-earth ions in compounds are close to those of the free ions and the total angular momentum (isospin), J, is a good quantum number. In CeCl 3 , the lowest energy level has J = 5/2, which splits into three Kramers doublets, of which m J = ±5/2 is the ground state. the interaction of chiral phonons with the isospin can be written as an effective "spin-orbit" type Hamiltonian [29,[49][50][51][52][53][54], where J α is the total isospin of unit cell α, L αn is the phonon angular momentum generated by mode n, and k αn is the coupling coefficient. The index α runs over all unit cells of the crystal and n over all chiral phonon modes.
For optical phonons at the Brillouin-zone center, the phonon angular momentum is homogeneous across unit cells and we can drop the index α. It is given by L n = Q n ×Q n , where Q n = (Q na , Q nb , 0) contains the normal mode coordinates of the two orthogonal components of a doubly degenerate phonon mode, Q na and Q nb , in the ab plane of the crystal. We can further treat the isospin in a mean-field approximation and replace J α by the ensemble average, J α . Taking into account only isospin components perpendicular to the doubly degenerate phonon modes, the ensemble average of the isospin is J α = 2|J|e z ( n −J − n J ), where e z is a unit vector along the c axis of the crystal. J is the isospin of a single cerium ion, of which there are two per unit cell. n ±J is the Fermi-Dirac distribution describing the occupation of the ground-state Kramers doublet. The magnetic moment per unit cell, m, is then given by where g J is the Landé factor. Eq. (9) is the same as (2) in the main text. The theoretical value of the prefactor in Eq. (9) for the m J = ±5/2 ground-state doublet is g ±5/2 = g J J(J + 1) = 2.54, which is reasonably close to the experimental value of 2.02 [55], showing that the orbital angular momentum is mostly unquenched. We can now rewrite Eq. (8) in terms of the magnetic moment, yielding where we redefined the coupling as K n m = k αn J α . If we look at the interaction one mode at a time, we can drop the index n, which yields Eq. (1) from the main text. The phonon Lagrangian, L, including the interaction in Eq. (10), can be written as where Ω a = Ω b ≡ Ω 0 is the eigenfrequency of the doubly degenerate phonon mode and m = |m|. In frequency space, the Lagrangian can be transformed according to Q a → Q Ωa exp(iΩt),Q a → iΩQ Ωa exp(iΩt), and where Q Ω = (Q Ωa , Q Ωb , 0), and D is the dynamical matrix. The spin-phonon coupling modifies the dynamical matrix as D(m) = D (0) + iD (1) (m) [9,13,46,47,[56][57][58][59][60][61]. In order to evaluate the effect of the interaction on the phonon frequencies, we compute the determinant of the dynamical matrix, which contains the spin-phonon coupling in its off-diagonal components, Solving the determinant close to the Brillouin-zone center (Ω → Ω 0 ) yields Eq. (3) from the main text, describing a splitting of the frequencies of right-and left-handed circular polarizations, Ω ± , of the doubly degenerate phonon mode, Without an external magnetic field, the energy levels of the m J = ±5/2 ground-state Kramers doublet are degenerate, there is no net magnetic moment per unit cell, and the phonon frequencies in Eq. (14) remain degenerate. Applying a magnetic field, B c, to the paramagnet splits the ground-state doublet, ∆E = E −5/2 − E 5/2 = 2g ±5/2 µ B B, and induces a magnetization given by m = 2g ±5/2 µ B n −5/2 − n 5/2 In the above equation, we have used the relation 1 exp(−x) + 1 − 1 exp(x) + 1 = tanh Inserting Eq. (15) into (14) yields Eq. (4) of the main text, ∆Ω(B) = 2Km = 4Kg ±5/2 µ B tanh g ±5/2 µ B B 2k B T (17) where ∆Ω s is the saturation phonon-frequency splitting between right-and left-handed chiral phonon modes. This equation allows us to extract the spin-phonon coupling from experimentally measured phonon splittings, K = ∆Ω s /(4g ±5/2 µ B ).
TABLE I. Comparison of effective magnetic fields (B eff ) induced by the optical and phonon inverse Faraday effects (IFE) in different materials. We further display the pulse fluences, durations, and spectral ranges. Previous theoretical work has most of the time displayed induced magnetizations arising from phonon orbital magnetic moments in terms of nuclear magnetons per unit cell. We convert the effective magnetic field as B = µ0m/Vc, where µ0 is the vacuum permeability and m is the magnetic moment per unit cell. Note that direct quantitative comparisons between the mechanisms are complicated by the fact that the effective magnetic field picture breaks down for very short (tens of femtoseconds) pulses.