{'fraction_non_alphanumeric': 0.04085459251043152, 'fraction_numerical': 0.013825067342734908, 'mean_word_length': 4.415302109402932, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Jahn-Teller states mixed by spin-orbit coupling in an electromagnetic field

Spin-orbit coupling plays a pivotal role in condensed matter physics. For instance, spin-orbit interactions affect the magnetization and transport dynamics in solids, while spins and momenta are locked in topological matter. Alternatively, spin-orbit entanglement may play an important role in exotic phenomena, like quantum spin liquids in 4d and 5d systems. An interesting question is how electronic states mixed by spin orbit coupling interact with electromagnetic fields, which may hold potential to tune their properties and reveal interesting physics. Motivated by our recent discovery of large gyrotropic signals in some Jahn-Teller manganites, here we explore the interaction of light with spin-mixed states in a d4 transition metal. We show that spin-orbit mixing enables electronic transitions that are sensitive to circularly polarized light, giving rise to a gyrotropic response that increases with spin-orbit coupling. Interestingly, photoexcited transitions that involve spin reversal are behind such gyrotropic resonances. Additionally, we find that the interaction with the electromagnetic field depends strongly on the relative orientation of the propagation of light with respect to Jahn-Teller distortions and spin quantization. We suggest that such interactions offer the opportunity to use electromagnetic waves at optical wavelengths to entangle orbital and spin degrees of freedom. Our approach, which includes a group-theoretical treatment of spin-orbit coupling, has wide applicability and provides a versatile tool to explore the interaction of electromagnetic fields with electronic states in transition metals with arbitrary spin-orbit coupling strength and pointgroup symmetries.


I. INTRODUCTION
Orbital degrees of freedom are an essential ingredient of the physics and chemistry of transition metal compounds [2][3][4][5]. The coupling of orbitals to spin, charge or lattice determines many properties of solids and molecules. In the presence of orbital degeneracy, the symmetry of non-linear molecules and solid state systems is broken spontaneously through the Jahn-Teller effect [6][7][8]. This phenomenon has far-reaching consequences in spectroscopic and chemical properties [9][10][11], and is also responsible for the emergence of nontrivial quantum effects, associated with the appearance of rotational quantization of vibronic states and geometric phases [12][13][14]. On the other hand, spin orbit interactions are key to new developments in classical and quantum computation and lie beneath new discoveries in condensed matter physics related to topological matter, such as the quantum spin Hall effect or the realization of topological insulators and Weyl semimetals [15][16][17] and Kitaev physics in quantum spin-liquids [18][19][20][21]. At the same time, there is a highly nontrivial interplay between spin-orbit coupling and the Jahn-Teller effect when t 2g states are partially filled, where entangled quantum spin-orbital states may emerge [22][23][24].
An interesting question is how electromagnetic fields interact with spin-orbit mixed states, which could pave the ground to explore quantum physics in these systems. Motivated by our recent finding of large gyrotropic signals in La 2/3 Ca 1/3 MnO 3 (originated by the different optical response to light of opposite handedness in the presence of Jahn-Teller distortions) [25], we present here a group-theoretical analysis to study this problem. Our formalism has general applicability and provides a useful route to extend the analysis to heavy transition metals in arbitrary point-group symmetries. In the following, we describe in great detail the interaction with an electromagnetic field of spin-orbit mixed t 2g − e g states in a 3d metal, which provides the clues to its generalization to other transition metals.
We first note that when dealing with spin-orbit physics in light transition metals, the mixing between t 2g and e g orbitals is usually neglected, since crystal-field splitting and exchange energies are much larger than spin-orbit coupling [22]. However, this approximation breaks down under particular conditions. To illustrate this point, we consider the Tanabe-Sugano diagram for the case of an ion with d 4 configuration in O h symmetry [1,26]. We see that for values of the crystal field 10Dq and the Racah parameter B that fulfill the condition (Dq/B) c ≤ 2.7 the ground state term is 5 E g , whereas for large enough Dq/B the ground state is 3 T 1g (Fig. 1a). In both limits, a good approximation is that spin-orbit interactions act only on the t 2g manifold (the orbital moment is quenched for e g states), and the spin-orbit mixing of e g and t 2g states (and, therefore, between 5 E g and 3 T 1g ) can be ignored, at least to first order in spin-orbit coupling. This results in the so-called T −P equivalence, where the spinorbit physics of the t 2g manifold can be described with an effective angular momentum L = 1, like for p orbitals [1].
However, this approximation breaks down for specific situations. For instance, in 4d/5d transition metals, where several energy scales (including spin-orbit coupling and crystal-field) are comparable [27,28], the t 2g − e g mixing can not be generally ignored [29,30], which is relevant to predict magnetic excitations in heavy transition metal oxides. Alternatively, the T − P equivalence  [1]. Along the ordinate axis, the energy E is shown relative to the Racah parameter B. The abscissa displays the ratio between the "differential of quanta" Dq and B (for octahedral complexes the crystal field energy is 10Dq). b) Splitting of the spectroscopic lines in O h , D 4h and spin-orbit coupling (SOC) point symmetries. In this work, we consider that the reduction to D 4h symmetry is driven by Jahn-Teller interactions.
fails when the difference in energies between spectroscopic terms becomes comparable with spin-orbital coupling. This may happen for ratios (Dq/B) c ≈ 2.7 in the Tanabe-Sugano diagram [1], which may occur in 3dtransition metals ( Fig. 1(a)). Alternatively, a reduction from O h to D 4h symmetry (e.g., induced by a Jahn-Teller instability, as we discuss below), may lead to spin-orbital mixing between spectroscopic terms 5 A 1g and 3 E g away from (Dq/B) c ≈ 2.7, see Fig. 1(b). As we show below, spin-orbital mixing between t 2g and e g states -which may be also induced by electronic correlations or structural distortions-enables optical transitions that can be probed by circularly polarized light. The sensitivity to circular polarization stems from terms of different spin multiplicity mixed by spin-orbit coupling (e.g., S =1 for 3 E g and S =2 for 5 A 1g in D 4h point symmetry), which allows optical transitions between states with different spin projections, which, otherwise, are absent in the absence of spin-orbit mixing. For the sake of conciseness, we restrict our discussion to 3d 4 ions in which spin-orbital mixing is induced by symmetry reduction to D 4h induced by Jahn-Teller instabilities (Fig. 1b). At the end we discuss briefly how the group-theoretical analysis can be extended to other transition metals with arbitrary spinorbit coupling strength and their interaction with electromagnetic radiation.
To describe the physics of spin-orbital mixing we con-sider a Hamiltonian that has the following form: where H S stands for on-site interactions and W represents the interaction of electrons with an electromagnetic field. We consider that the interaction with light induces the hopping with amplitude α ij of the fourth electron in the d 4 ion to neighboring d 3 sites i, j in the lattice. The on-site Hamiltonian H S contains diagonal terms denoted by E ψ |ψ i ⟩⟨ψ i | and off-diagonal V ψϕ |ψ i ⟩⟨ϕ i | contributions, coming from vibronic couplings and spinorbit interactions. The Hamiltonian can be formally expressed in terms of irreducible representations |ψ i ⟩, |ϕ i ⟩ ∈ The terms of this basis are expressed as linear combinations of Slater determinants of monoelectronic orbitals t 2g ∈ (|xy⟩, |yz⟩, |xz⟩) and e g ∈ (|x 2 − y 2 ⟩, |z 2 ⟩) that respect the Pauli exclusion principle and the point-group symmetries in orbital (D 4h , due to Jahn-Teller instabilities) and spin spaces (in Appendix A we give a detailed derivation of this basis and the development of the corresponding Slater determinants). We note that (B 1g , B 2g ) and (A 1g , A 2g ) are degenerate (possibly this accidental degeneracy is broken if we consider higher-order relativistic corrections) and they are lumped together in the ba-sis |ψ i ⟩. Additionally, the presence of magnetic fields can lift the 2S + 1-fold degeneracy of terms 3 A 2g (S = 1), 5 A 1g (S = 2) and 5 B 1g (S = 2), giving a total of 19 states for the full dimensionality of the |ψ i ⟩ basis. The on-site Hamiltonian H S can be decomposed as: where H 0 includes the energy splitting between O h pointgroup terms 5 E g and 3 T 1g due to crystal field and exchange interactions (see Fig. 1b), while H JT takes into account interactions with Jahn-Teller modes and H SO is the spin-orbit coupling contribution. The paper is organized as follows. In Sec. II A, we describe the spontaneous breaking of orbital degeneracy driven by Jahn-Teller instabilities in a d 4 ion under O h symmetry. In this situation, 5 E 2g and 3 T 1g electronic states interact with doubly degenerate E g Jahn-Teller modes, resulting in E ⊗ e and T ⊗ e vibronic interactions. As a result, the point-group symmetry is reduced to D 4h and the terms split into 3 A 2g , 3 E g , 5 A 1g and 5 B 1g states (see Fig. 1b). In Sec. II B we study the point symmetries in orbital and spin spaces related to the spin-orbit operator. The combination of Jahn-Teller and spin-orbit interactions split further these terms into irreducible representations Fig. 1b and Appendix A. In Sec. III, we describe the interaction of electrons with electromagnetic fields. As a result of this interaction, we assume that electrons hop between neighboring sites. We focus on the problem of a lattice in which the fourth electron of an isolated Jahn-Teller d 4 ion hops to d 3 nearest-neighbours. We assume that the solid is a transition metal oxide with perovskite structure with ABO 3 chemical formulation, in which generally A is a rare-earth element, B a transition metal and O is oxygen. These systems form a large family of materials that includes La 2/3 Ca 1/3 MnO 3 , with a broad diversity of physical properties, including magnetism, ferroelectricity or superconductivity [31][32][33][34]. The undisturbed perovskite is formed by octahedral unit cells with the metal B sitting at the center of an octahedron formed by six ligand oxygen anions with O h point symmetry [35,36], see Fig.  2(a). Using the formalism of two-center Slater-Koster integrals, we derive analytic expressions for the lightinduced hopping amplitudes between lattice sites. The perturbative analysis discussed in Sec. IV demonstrates that spin-orbit coupling and intraatomic t 2g − e g mixing are essential to the appearance of gyrotropic responses, and that the latter involve photoexcitations in which one of the spins is inverted. Remarkably, this observation opens the possibility of using electromagnetic fields to manipulate spins via the mechanism described here. Subsequently, in Sec. V, we analyze the electronic response to circularly polarized electromagnetic waves. For that purpose, we analyze the density of 5 B 1g states from the imaginary part of quantum propagators of the different electronic orbitals and obtain expressions for their spectral functions. In Sec. VI we analyze these spectral func-tions in circularly polarized electromagnetic fields as a function of the relative strength of Jahn-Teller and spinorbit interactions. From this analysis, we extract information about the gyrotropic responses, by which the polarization of light is changed as a result of the interactions with t 2g −e g spin-orbit mixed states. Finally, in Sec. VII, we summarize the main results and discuss perspectives of further work, especially the possibility of entangling spin and orbital degrees of freedom using electromagnetic waves, which could be relevant in the framework of nontrivial quantum states in other systems, including heavy 4d-5d transition metals.

A. Jahn-Teller interactions
We first derive the Hamiltonian terms for the interaction of electron orbitals with Jahn-Teller modes. Under O h point-group symmetry they can interact with two degenerate representations for E g Jahn-Teller vibrational modes corresponding to tetragonal modes Q Egu = Q 3 = 2∆z − ∆x − ∆y and orthorhombic modes Q Egv = Q 2 = √ 3(∆x − ∆y), respectively (see Fig. 2). We thus need to solve the E ⊗ e and T ⊗ e problems to derive analytic expressions for the corresponding vibronic interactions between Jahn-Teller E g modes and doubly degenerate E g and triply degenerate T 1g electronic states in O h symmetry [37]. A convenient way to derive these expressions is to write the Jahn-Teller modes in terms of an angle ϑ as Q 3 = cos ϑ and Q 2 = sin ϑ [38,39]. Using Pauli matrices υ i in the pseudospin space of [ 5 A 1g , 5 B 1g ] states in D 4h symmetry, the E ⊗ e Jahn-Teller interaction can be expressed as: where F E and G E are linear and quadratic vibronic constants, and δϑ represents perturbative Q 2 orthorhombic distortions that will be described below. The dependence of Jahn-Teller modes on ϑ defines a potential energy surface, which, in the case of harmonic approximation (i.e., G E = 0), defines a "Mexican hat" [40]. In solids, however, anharmonic contributions are usually relevant and quadratic constants (so that G E ̸ = 0) must be included in Eq. (3). As a result, the surface potential warps producing three minima at ϑ n = 2nπ/3, which correspond to tetragonal elongations of the octahedral cell along z, y and x axes for n = 0, 2, 1, which stabilize the occupation of d z 2 , d y 2 and d x 2 orbitals, respectively [41][42][43] (see also Appendix B for a detailed description of the Jahn-Teller Hamiltonian and the vibronic interactions). Actually, the stabilization of tetragonal Jahn-Teller distortions in solids has been confirmed experimentally in a large number of compounds [44], including manganites [45][46][47][48].
FIG. 2. a) Depiction of a transition metal in an octahedral coordination with six oxygen ligands. b) Tetragonal Jahn-Teller distortion corresponding to the mode Q = Q3. c) Orthorhombic Jahn-Teller distortion corresponding to the mode Q = Q2. d) Superposition of the two distortion modes corresponding to Q = Q3 cos ϑ + Q2 sin ϑ.
Now we derive the analytic expressions for vibronic interactions involving T 1g electronic states, which can interact with E g and T 2g representations of Jahn-Teller modes [37]. Consequently, such derivation requires solving the T ⊗ e and T ⊗ t problems [41]. We assume, however, that the T ⊗t contribution is negligible, since in our case lattice deformations are predominantly driven by Jahn-Teller instabilities of electrons in e g states. There-fore, we only consider the contribution of T ⊗ e to the Hamiltonian as follows: where F T is the vibronic coupling constant for T ⊗ e and λ 2 and λ 8 are Gell-Mann matrices, defined as As mentioned above, anharmonic lattice contributions stabilize tetragonal elongations along the three main axes of the ABO 3 octahedral cell units, denoted as ϑ n = 2nπ/3. This entails breaking the degeneracy of 5 E 2g electronic states in O h into 5 A 1g and 5 B 1g terms in D 4h symmetry, being the latter lower in energy for the elongated tetragonal distortions (Fig. 1b). In this work we also consider small orthorhombic Jahn-Teller distortions corresponding to Q 2 modes, which modify perturbatively the elongated tetragonal distortions, so that the angle in Q 2 − Q 3 space is mapped to ϑ n → ϑ n + δϑ n . As can be inferred from the E ⊗ e Hamiltonian ((2)), these orthorhombic perturbations induce non-diagonal transitions between 5 A 1g and 5 B 1g terms. On the other hand, for the t 2g sector, the reduction to D 4h symmetry splits the 3 T 1g term into single-degenerate 3 A 2g and doubly degenerate 3 E g terms, being the latter lower in energy. In this case, it can be shown that, to first-order, ϑ n +δϑ n perturbations do not mix 3 A 2g and 3 E g terms.

B. Spin-orbit coupling
In order to compute the matrix elements of spin-orbit coupling in D 4h symmetry we use the operator equivalent method [1,49]. In this approach, the spin-orbit operator is defined by linear combinations V Λ λq = T Λ λ S 1 q , where T Λ λ corresponds to irreducible representations Λ in the orbital space with basis λ and S 1 q corresponds to irreducible representations D (1) q in the spin-rotation group. Since the T ⊗e Hamiltonian term in Equation (3) is relatively small -i.e., F T δϑ ⪅ ξ SO , where ξ SO is the spin-orbit coupling constant-we can study the orbital space T Λ λ in the D 4h point-group. Therefore, while T Λ λ transforms as T 1g in O h symmetry, the reduction to D 4h implies that T Λ λ transforms according to irreducible representations A 2g with spatial symmetry ν and E g with spatial symmetries κ, µ (see Appendix A 2 b for a definition of these symmetries). On the other hand, the spin part S 1 q is expressed using spherical coordinates, with quantum numbers q = 0, ±1. Taking this into account, it can be demonstrated that the spin-orbit coupling Hamiltonian can be expressed in terms of operator equivalent matrices V Eg κ±1 , V Eg µ±1 and V A2g ν0 as follows: where the sign of q is denoted by a breve symbol, i.e., q = −q. Then we apply the Wigner-Eckart theorem to V Λ λq , which implies working with reduced matrices ⟨ΓS||V Λ ||Γ ′ S ′ ⟩ [50]. In our case, these are 4×4 matrices defined in terms of the irreducible representations of the D 4h point-group { 3 E g , 3 A 2g , 5 A 1g , 5 B 1g }, which are expressed as follows (see Appendix C for the details of this derivation): Once the reduced matrices are computed, the spinorbit elements can be found using Clebsch-Gordan coef-ficients as follows: With these relations, and taking the basis where again the breve symbol denotes the sign of spin quantum numbers, one can write the full 19x19-dimensional spin-orbit matrix as: In this expression, solid lines separate the matrix elements corresponding to spectroscopic terms We note that matrix Eq. (9) is represented for the quantization of ⃗ L and ⃗ S along the same axis. However, in general, the quantum spin axis can be oriented along arbitrary directions with respect to ⃗ L. Therefore, it is convenient to apply appropriate rotations R in the spin space to orient the spin quantization axis along arbitrary directions defined byn as follows: This rotation is characterized by an axist = (ẑ×n)/|ẑ×n| and a rotation angle θ = arccosẑ ·n.
If we definen = (sin θ cos ϕ, sin θ sin ϕ, cosθ), we havê t = (− sin ϕ, cos ϕ). Then, we obtain the following spinrotation matrices for the cases S = 1 and S = 2: We use these matrices to compute the spin-orbit elements of Eq. (8) for arbitrary directions of the quantized spin axis.

III. THEORY: INTERACTION WITH ELECTROMAGNETIC FIELDS
A. Light-induced electron transfer between lattice sites FIG. 3. Diagram of the electron transfer between neighbouring sites mediated by oxygen ions, induced by the interaction with an electromagnetic field. In the process, there is a first transition from oxygen p-orbital to neighbouring manganese d'-orbital. A second transition involves a transfer from a manganese d-orbital to an oxygen p-orbital.
So far we have considered on-site interactions of electronic orbitals with the crystal field, Jahn-Teller vibra-tional modes and atomic spin-orbit coupling. In the following, we describe their interaction with electromagnetic fields, which we assume induce electron transfer between neighboring sites in the perovskite lattice. The idea of light-induced electron transfer has been proposed, e.g., in some manganites, where optical energy excitations have been associated to polaronic transport due to intersite e g − e g photoinduced transitions [25,51,52]. Since the cation separation in perovskites (≈ 4Å) is large for significant direct overlap [53][54][55], we consider the electron transfer dominated by hopping through p orbitals of oxygen. As depicted in Fig. 3, we consider the transfer between neighbouring d 4 and d 3 ions, which can be described as: In order to describe this transfer, we define the following many-electronic wavefunctions for the two neighbouring d ions and the oxygen ligands: Eq. (14a) corresponds to the initial configuration d 4 p 6 d 3 , where the d 4 ion is described by some of the 2S+1 Γ representations discussed in Sec. I, whereas 4 A 2g is the ground state for the d 3 ion, according to the Tanabe [44]. Therefore, in the presence of an electromagnetic field, the orbitals p and d are coupled, so that the states |Ψ⟩ and |Ψ ′ ⟩ are perturbed by the intermediate states |Φ pw ⟩ as follows: where t pd /2 is the p − d hopping amplitude induced by the electromagnetic field, which allows nonzero matrix elements between states |Ψ⟩ and |Ψ ′ ⟩.
We treat the interaction with light to first order, so that the amplitude of the light-induced transfer requires the computation of electromagnetic matrix elements that involve two-center integrals including the vector potential −ı ⃗ ∇ (defined in the Coulomb gauge): whereê q indicates the hopping direction in the lattice,ε is the unit vector along the orientation of the vector potential, a is the lattice parameter and ψ, ϕ pw describe monoelectronic orbitals in the transition metal and oxygen, respectively, that are involved in the photoinduced transfer. We note that although the spectroscopic terms are given as combinations of Slater determinants, the vector potential in Eq. (16) is a one-body operator that acts only on the monoelectronic orbital where the transferred electron resides (see Appendix D for a detailed discussion of how one-body operators act on the many-electron wavefunctions). The matrix elements shown in Eq. (16) are therefore expressed in terms of monoelectronic functions ψ and ϕ pw . This derives from the properties of the one-body potential, whereby matrix elements such as −ı⟨Φ pw | ⃗ ∇|Ψ⟩, where |Φ pw ⟩, |Ψ⟩ are many-electron functions described by Eq. (14), can be rewritten as −ı⟨ψ| ⃗ ∇|ϕ pw ⟩, where |ϕ pw ⟩, |ψ⟩ describe monoelectronic orbitals.
On the other hand, while the expression in Eq. (16) corresponds to a transfer from a p to a d orbital, the d to p transition is described by its complex conjugate (P ψ qεw ) * . Interestingly, it can be shown that expressions like ∂εφ pw (withε alongx,ŷ, orẑ) appearing in Eq. (16) can be expressed as linear combinations of Slater-Koster coefficients (see Chapter 7 in Ref. [31] for a derivation). For instance, for the vector potential alongε||x, we make use of the following coefficients: where ψ z 2 , ψ xy are wavefunctions for the monoelectronic states |z 2 ⟩ and |xy⟩, andψ s ,ψ z 2 andψ xy are effective wavefunctions which have the same symmetries as s, d z 2 and d xy orbitals (see Ref. [31]). Table I displays all nonzero matrix elements for the vector potential along the three directions in space in terms of the coefficients (sdσ), (ddσ), (ddπ) and (ddδ). The hopping amplitudes α ψiϕj q between ψ i and ϕ j orbitals located at neighboring sites (i, j such that ⃗ r i −⃗ r j ∥ e q ) are calculated perturbatively, taking into account the p − d hopping t pd and the charge transfer energy ∆ CT between p and d orbitals [44]: Since the electromagnetic field cannot interact directly with spins, the matrix elements of the electromagnetic operator Wε have the following form: where i, j refer to neighboring locations in the lattice. In the next section we explain how the hopping amplitudes depend on the light polarization, which is described by the unit polarization vectorε along an arbitrary direction.

B. Cooperative Jahn-Teller effects
Although we address the dynamics of electron transfer from isolated Jahn-Teller ions, we incorporate cooperative effects, known to be relevant in solids [39,44,[56][57][58][59][60][61]. The reason is that the dynamics of ions is much slower than the electronic transfer rates, so that we assume that cooperative effects restrict the possible Jahn-Teller deformations of the neighbouring sites where the transferred electron can jump into (see Fig. 4). As discussed in Sec. II A, we consider Jahn-Teller modes of the d 4 ion described by angles ϑ n = 2nπ/3 + δϑ, n = 0, 1, 2 and δϑ ≪ 2π/3. In consequence, there are three possible orientations for the transfer across the six oxygen anions surrounding the initial d-site, namely along ±x, ±ŷ or ±ẑ. Then, cooperative effects are incorporated by imposing restrictions on the hopping from an initial d 4 ion with tetragonal distortion alongẑ (Jahn-Teller mode with angle ϑ 0 ) to neighboring sites along the three directions (see Fig. 4). For instance, when δϑ > 0 there is a contraction along y axis, forcing neighbours on the xy plane to be distorted along the y direction with δϑ < 0. On the other hand, when a site distorted along z has δϑ < 0, there is a slight contraction along x, so their neighbours are distorted along x with δϑ > 0. Our model considers the orbital ordering with maximum entropy, which consists in having half of the octahedra distorted along a particular direction (chosen to be z) and the rest is equally distributed among elongations along x and y [62]. The resulting orbital ordering is depicted in Fig. 4. One of such transfers is illustrated by labels "1" and "2", where the electron jumps through an intervening oxygen. Each vertex represents a transition metal, around which the octahedron elongates along the solid lines. The electron is initially located at d 4 sites whose distortion is along z. Due to cooperative effects, the d 3 sites around the initial d 4 site are elongated along directions perpendicular to z.
C. Hopping amplitudes for circularly polarized light

Left-and right-handed basis for circularly polarized light
We describe the polarization of light by a complex vectorε ∈ C 3 normalized toε ·ε * = 1: Here E 0i is the amplitude of the i-th component of the electric field and ϕ i the correspondent phase, defined in this expression to keep ϵ x real. Since light is a transverse wave, we need only two complex vectors to define a basis for the polarization, both orthogonal to propagation directionk. For circularly polarized light we use left-handed and right-handed polarizations, which, in the case of wave propagation alonĝ k =ẑ are defined as: For arbitrary orientations of the propagation of light, we use the rotation matrix R to find the new basiŝ ϵ L,R = Rς L,R for the polarization. This rotation is characterized by an axisû = (ẑ ×k)/ sin α and an angle of rotation cos φ =ẑ ·k. Since every unit vectork can be described using polar and azimuth angles, α, β, then k = (sin α cos β, sin α sin β, cos α),û = (− sin β, cos β, 0) and φ = α. Then, the rotation matrix for arbitrary wave propagation can be defined as: where a contracted notation for trigonometric functions is used, namely, s x = sin x and c x = cos x. Finally the polarization vector for arbitrary wavevector orientation has the following expression:

Electromagnetic response and time-reversal symmetry
Left-and right-handed polarizations are related to each other by complex conjugation,ε * L =ε R . With this relation we can deduce that (P ψ qLw ) * = −P ψ qRw (see Eq. (16)) if ψ(⃗ r) ∈ R, ∀⃗ r ∈ R 3 , in other words: ψ has real spatial symmetry, which is the case of the basis used here. Time reversal involves complex conjugation (since light acts only on the orbital angular momentum) and the interchange of the initial and final orbitals in the hopping, so that we have (α ψϕ qL ) * = α ϕψ qL = α ψϕ qR (27) which means that KW L K † = W R , being K the complex conjugation operator and W L,R the electromagnetic operator for left-and right-handed light. We note that the Hamiltonian terms (see Eq. (2)) H 0 and H JT are expressed as real matrices, while H SO is complex. Thus, in the absence of spin-orbit coupling (ξ SO = 0), we have KH L K † = H R and, as a consequence, the gyrotropic signal is zero. A nonzero gyrotropic response (a different response to electromagnetic waves of opposite handedness) arises only when ξ SO ̸ = 0, which implies KH R K † ̸ = H L . This conclusion does not depend on the basis, since it holds even when the wavefunctions are not real, for instance, when they are expressed in spherical basis. Indeed, transforming from spherical to a real basis involves a unitary transformation U, so that the relation between W L and W R is whereK is also an antiunitary transformation that keeps invariant H 0 and H JT . This confirms that a change of basis does not break the time-reversal invariant relation between H L and H R when ξ SO = 0. A2g. The shadowed area corresponds to the parameter space where the spin-orbit mixing between t2g and eg states and, consequently, the gyrotropic response, are both strong.

Orbital-selective gyrotropic responses in broken time-reversal symmetry
The origin of the gyrotropic responses can be traced back to the transfer induced by light between specific orbitals in the d manifold. To shed light on this issue, it is convenient to express the polarization asε L = (a, b, c) andε R = (a * , b * , c * ). With this, the transfer amplitudes t ψϕ qε in Eq. (27) can be expanded as products of P ψ qεŵ integrals (see Eq. (16)) as By inspection of Eq. (22) one realizes that any transfer involving hopping between e g and t 2g orbitals at neighboring sites is forbidden, since S = 2 and S ′ = 1 (or viceversa). Then, by taking into account Eq. (29) and the nonzero transfer integrals P ψ qεw displayed in Table  I, one can verify that all hopping amplitudes involving hopping between neighboring e g − e g orbitals are real and are consequently time-reversal invariant. Therefore, the transfer between e g − e g orbitals cannot give a gyrotropic response, at least to first order perturbation in the electromagnetic field. On the other hand, the lightinduced transfer between neighboring t 2g − t 2g orbitals has complex amplitude and breaks time-reversal symmetry, causing distinct electromagnetic responses for light of opposite handedness. As a consequence, both the presence of spin-orbit coupling and intersite t 2g − t 2g transfer are key ingredients to have a gyrotropic signal.
In the light of the previous discussion, one expects a strong influence of the t pd hopping integral on the gyrotropic signal. Since the overlapping integrals (sdσ) and (ddσ) are significantly larger than (ddπ) and (ddδ), the energy of t 3 2g e 1 g states (corresponding to 5 A 1g and 5 B 1g representations) is influenced much more strongly by t pd than t 4 2g states with 3 E g and 3 A 2g representations (see Fig. 5). It is then expected that as the value of t pd grows, the eigenvalues of t 4 2g states will cross eventually the eigenvalues of t 3 2g e 1 g states, producing a strong spinorbit mixing and enhancing the gyrotropic signal (shadowed area in Fig. 5). In the present problem, we have verified numerically that this condition is fulfilled for values in the range t 2 pd /∆ CT ∼ 0.3eV − 1eV.

IV. PERTURBATION ANALYSIS OF THE GYROTROPIC RESPONSE
As discussed previously, a gyrotropic signal requires photoinduced transfer between adjacent t 2g −t 2g orbitals. In consequence, the unperturbed ground state of the isolated Jahn Teller ion 5 B 1g (t 3 2g e 1 g ) has to be excited to a t 4 2g configuration to activate this transfer channel. Here we derive a perturbation analysis in spin-orbit coupling and orthorhombic modes to understand the electronic transitions that contribute to the gyrotropic signal. First of all, we introduce the notation |i 2S+1 ΓγM ⟩, which indicates the irreducible representation of the wavefunction at the i-th site in the lattice. The introduction of intersite hopping by interaction with the electromagnetic field breaks the degeneracy between the same state at different sites, originating 0-th order eigenstates denoted by: To continue with the perturbed states we have to understand how the Hamiltonian acts on those 0-th order eigenstates. Recalling that the Hamiltonian has on-site H S and inter-site electromagnetic W terms, we have: A difficulty arises to compute such matrix elements due to cooperative Jahn-Teller effects. In particular, a given irreducible representation may contain different wavefunctions at adjacent sites in the lattice. For instance, the 3 A 2g term corresponds to a |ζηττ | Slater determinant if the distortion is along the z direction, but if this distortion is along y it corresponds to |ζητη|. Thus, the computation of matrix elements described by Eq. (31), necessary to determine the perturbed eigenstates, is challenging. We sort out this difficulty by approximating the matrix elements, for instance, In this example, which can be generalized to arbitrary elements, every term is not determined exactly, but it gives a reasonable estimate of the contributions coming from intersite hopping, spin-orbit coupling and Jahn-Teller interactions.
Now we develop the perturbative analysis. Since we consider the dynamics of an electron initially located in a tetragonally elongated site, we are therefore interested in calculating the transition rates between 5 B 1g (corresponding to a t 3 2g e 1 g configuration) and 3 E g or 3 A 2g terms (both corresponding to a t 4 2g configuration). By the effect of spin-orbit coupling and orthorhombic distortions, the non-perturbed wavefunctions |α 3 A 2g νM ⟩ and |α 3 E g t q M ⟩ become, respectively, |α 3 A 2g νM ⟩ • and |α 3 E g t q M ⟩ • , where spherical harmonics are used to describe the orbital components t q of the wavefunctions (see Eq. (A10) for the definition of t q ). The matrix elements are then approximated to first order in spin-orbit coupling and orthorhombic modes as follows: where E EB and E EA are, respectively, the energy gaps between 3 E g and 5 B 1g and between 3 E g and 5 A 1g , while E AB and E AA are the analogous gaps corresponding to 3 A 2g instead of 3 E g . These energy gaps determine the degree of orbital mixing between e g and t 2g states. The vibronic constants (F E , G E , F T ), spin-orbit coupling (ξ SO ) and intersite hopping (t pd ) are also included in the expressions above. According to this perturbational analysis, the different transitions contributing to the gyrotropic signal are sketched in Fig. 6. We first note that spin-orbit corrections connect 5 B 1g with 3 E g , giving rise to the matrix element of Eq. (32a). On the other hand, Eq. (32b) stems from inter-site hoppings and orthorhombic corrections connecting 5 B 1g and 5 A 1g followed by spin-orbit mixing of 5 A 1g with 3 E g , while Eq. (32e) takes account of inter-site hoppings and orthorhombic corrections connecting 5 B 1g and 5 A 1g plus spin-orbit coupling between 5 A 1g and 3 A 2g . Finally, Eq. (32c) and Eq. (32d) come from spin-orbit interactions within the 3 E g subspace with M = 0, where the degeneracy of the wave-functions, which is preserved by spin-orbit coupling, lifts under the action of orthorhombic modes. An inspection of these expressions allows to extract the following conclusions: • A relevant gyrotropic signal appears when the e g − t 2g spin-orbit mixing is large. According to Fig. 5, this happens when the gap E EA between 5 A 1g and 3 E g is reduced by the effect of light induced transfer through the t pd hopping integral. In this case, the strong reduction of E EA entails an enhancement of contributions described by Eqs. Eq. (32b)-Eq. (32e).
• All amplitudes involving 5 B 1g , 3 E g and 3 A 2g described by Eq. (32) imply transitions between t 3 2g e 1 g and t 4 2g configurations, where one of the spins is inverted during the transition. The only exception is the transition described by Eq. (32e), which is a second-order correction in ξ SO δϑ, i.e., it requires the simultaneous action of spin-orbit and orthorhombic interactions. Since δϑ is small, and considering typical values for the vibronic constants (F T , F E , G E . see section Sec. VI), the contribution of this term is negligible. Therefore, we conclude that the observation of a large gyrotropic signal is fundamentally contributed by transitions that in-volve a spin reversal.
• The perturbative influence of orthorhombic Jahn-Teller modes is described by the parameter δϑ. For small values of the 5 A 1g − 3 E g gap, i.e., E EA ⪅ ξ SO , the predominant transition contributing to the gyrotropic signal is given by Eq. (32b). In this case, in addition to spin-orbit coupling, the hopping between neighboring t 2g − t 2g states and orthorhombic modes enhance the gyrotropic signal. However, since the inter-site hopping is far larger than the energy of the orthorhombic distortions, the dependence of the gyrotropic response on δϑ is very weak. On the other hand, for large enough values of the gap E EA > ξ SO , the transition described by Eq. (32a) becomes predominant, but its amplitude is significantly smaller than for the case E EA ⪅ ξ SO . We can then conclude that the role of orthorhombic perturbations is minor, at least in the regime where F T δϑ ⪅ ξ SO and, therefore, the gyrotropic response is dominated by transitions between 5 A 1g and 3 E g , where the wavefunctions are perturbed by spin-orbit coupling.
We end this section by discussing the effects of the geometry on the gyrotropic signal, stemming from the relative orientations of light propagation and spin quantization, takingẑ as the orientation along the tetragonal distortion. An inspection of Eq. (29) reveals that for light propagating alongk =ẑ, namely, when light propagates along the distortion, the allowed gyrotropic hopping channel is mediated by ζ − η orbitals. In contrast, when light propagates perpendicular to the Jahn-Teller distortion the allowed gyrotropic hopping is η − τ for k =x, while for propagation alongk =ŷ the gyrotropic hopping is mediated by τ − ζ orbitals (see Eq. (A3c), Eq. (A3c) and Eq. (A3e) for a definition of the t 2g orbitals ζ, η and τ ). As a result, the magnitude of the gyrotropic signal strongly depends on both spin axis and light propagation. The reason is as follows. In general, for a given couple of t 2g orbitals in the hopping channel, the matrix elements of the angular momentum are nonzero only if the direction of the momentum component is contained in both spatial symmetries of the orbitals. For instance, for τ ∼ xy and η ∼ xz orbitals, the only non-vanishing element is ⟨τ |l x |η⟩ ̸ = 0. In addition, for a given pair of coupled orbitals in the hopping channel, it can be shown that the spin axis has to be oriented along the component of the nonzero matrix element to have a gyrotropic signal. For instance, for light propagating along the Jahn-Teller distortion, i.e., fork =ẑ, the only gyrotropic channel is ζ − η, Therefore, if the spin is quantized along x, then ⟨ζ|l x |η⟩ = 0 and the gyrotropic signal is completely extinguished. Numerical calculations, discussed below, have been performed to study systematically the effect of geometry on the response to circularly polarized light.
In the following, we introduce the formalism to evaluate the gyrotropic response (section Sec. V), which we use to perform numerical calculations based on the exact diagonalization of the full Hamiltonian. In Sec. VI we analyze the influence of spin-orbit coupling and orthorhombic Jahn-Teller modes on the gyrotropic response, which confirms the general tendencies discussed in this section. 6. Sketch of the transitions allowed by intersite hopping induced by light (green solid lines), spin-orbit coupling (brown dashed lines) and Jahn-Teller orthorhombic distortions (blue dashed-dotted lines). Thicker arrows indicate stronger interactions, corresponding to transitions between 5 B1g and 5 A1g mediated by intersite hopping and between 5 A1g and 3 Eg mediated by spin-orbit coupling.

V. RESPONSE TO CIRCULARLY POLARIZED ELECTROMAGNETIC WAVES
We consider the excitation of an electron located initially in a d 4 site distorted along z, see Fig. 4. In the presence of an electromagnetic field, this electron is transferred to any of the six nearest neighboring d 3 sites in the lattice. As explained in section Sec. III A, we assume that cooperative effects induce orbital ordering around the initial d 4 site, so that the site that receives the transferred electron can only deform along particular orientations, as shown in Fig. 4. In the calculations, the orbital ordering is extended periodically throughout the solid. To compute the dynamics, we suppose that the system is prepared in an ensemble where ψ refers to state 5 B 1g in the central site, which has the lowest energy (see Fig. 1b). Here, p ψ is the relative weight assigned to each spin projection allowed by the irreducible representations of the corresponding many-electron wavefunctions. The values of p ψ are indeed obtained for each specific case after diagonalization of the full Hamiltonian. The ensemble in Eq. (33) evolves in time as which allows us to compute the quantum propagator [63,64] where the Heaviside function Θ(t) accounts for causality and ⟨O⟩ is the thermal average of the operator O.
Here ρ is the density matrix for a thermal bath where Z is the partition function and the second equality holds if {|k⟩} is an eigenbasis of the Hamiltonian. We can then express the quantum propagator in the following way Here we have introduced the identity using the Hamiltonian eigenbasis labeled with m. Defining Ω km = |ω k − ω m | and approximating e −βω ≈ 1 for ω < β −1 and e −βω ≈ 0 for ω > β −1 we can rewrite the propagator as follows preserving only the terms such that ω k ≫ ω m : The first term in Eq. (39) corresponds to the retarded propagator, while the complex conjugate term is the advanced propagator. In order to compute the spectral response, we use the Heaviside function in the frequency domain Using the expression Eq. (40) we can write the spectral representation of the propagator in frequency domain where Ω km denotes the frequency eigenvalues of the full Hamiltonian. In the limit η → 0 + , we have which has the form of a density of states, which we denote as ς(ω), while the parameter η is related to the lifetime of the excited states. After some algebra, the spectral function can be rewritten as In the next section, we use this function to evaluate the gyrotropic response when time-reversal symmetry is broken.

A. Calculation of the spectral functions for circularly polarized light
For the calculation of the spectral functions defined by Eq. (43) we have solved the full Hamiltonian (Eq. (1)) to compute the eigenvalues. The spectral functions have been obtained for left-(ς L ) and right-(ς R ) circularly polarized light, by calculating the hopping amplitudes as described in Sec. III C. From these functions, we have built non-gyrotropic (ς ng ) and gyrotropic (ς gy ) spectral functions in frequency space, which give account, respectively, of the dynamic responses that are insensitive and sensitive to the handedness of the polarization of light. These functions are defined as follows: Finally, we define a function that integrates the gyrotropic signal over the analyzed spectral range (0eV − 3.5eV): The numerical calculations were carried out by setting the vibronic constants to F E = 450meV, F T = 130meV and G E = 20meV. These values are in agreement with the Jahn-Teller splitting observed for e g and t 2g electrons in 3d elements [44]. On the other hand, the charge transfer gap has been set to ∆ CT = 4eV [44], the damping factor to η = 180meV (see Eq. (43)) and the p − d hopping to t pd = 1.2eV (Eq. (21)). Finally, the Slater-Koster coefficients were set to (sdσ) = 1, (ddσ) = 0.82, (ddπ) = 0.29 and (ddδ) = 0.07.
We studied different geometries by varying the relative orientation of the light propagation, Jahn-Teller distortions and spin quantization. By way of illustration, the spectral functions ς ng and ς gy displayed in Fig. 7 were calculated for three different geometries, which are sketched in the top panels of each column. The spectral functions were computed for four different values of  Fig. 1b). The spectra computed for ∆ = 1.36eV, 0.89eV, 0.44eV, 12.5meV are displayed in descending order for each column (labelled a), b) and c)). The different functions have been computed in different geometric configurations, as sketched on the top of each column. The visible part of the spectrum has been shadowed and divided in two parts, below and above 550nm. This division enables an easier comparison with the experimental magneto-optical spectra reported in Ref. [25]. Those experiments show two absorption peaks centered, respectively, at wavelengths < 550nm and > 550nm, where only the latter gives rise to a gyrotropic response. The spectra computed in column b) for a gap ∆ = 1.19eV are in agreement with the experimental spectra reported in [25].
the energy gap ∆ (as indicated in the panels of Fig. 7), while the spin-orbit coupling was set to ξ SO = 20meV. The gap ∆ is defined as the energy difference between the unperturbed 5 A 1g and 3 E g terms (see Fig. 1b), which gives an estimation of the degree of t 2g −e g mixing before the introduction of the electromagnetic field. All parameters, including ∆, were chosen to work in a region of the Tanabe-Sugano diagram appropriate for manganese ions, for which the crystal field is 10D q ≈ 2eV and the Racah parameter is B ≈ 0.11 − 0.13eV [26,65,66]. On the other hand, in Fig. 8, the integrated gyrotropic signal described by N gy (Eq. (45)) is mapped as a function of spin-orbit coupling ξ SO and orthorhombic perturbations δθ for each value of ∆. Panels in Fig. 8 are organized in the same way as in Fig. 7, i.e., each column corresponds to each of the geometries sketched on the top.

B. Discussion of gyrotropic and nongyrotropic responses
We first discuss the nongyrotropic spectra described by functions ς ng shown in Fig. 7. First of all, we observe that the structure of resonances remains virtually unchanged, as long as the geometry is fixed, regardless of the values of the other parameters. In addition, a comparison between the ς ng spectra displayed in Fig. 7 (a) and (c) shows that the nongyrotropic response does not depend on the direction of the spin quantizationŝ, provided that the relative orientations of light propaga- tion and distortions remain the same. This is an indication that the nongyrotropic spectra are contributed essentially by transitions between 5 B 1g and 5 A 1g states. This observation is supported by the fact that the corresponding spin-orbit elements are zero for these states (see Eq. (9)), which explains why the ς ng spectra remain unchanged as the axis of spin quantization changes. Therefore, the structure of resonances observed in ς ng arises basically from the interactions between 5 B 1g and 5 A 1g states, mostly through intersite hopping induced by light (Fig. 5). More specifically, such transitions connect 5 B 1g and 5 A 1g bonding/antibonding states emerging from light-induced hybrization, which explains the structure of the peaks in the spectra (Fig. 7). The resonances located at lower energy, below the visible range (shadowed areas in Fig. 7) are mostly contributed by transitions between hybridized 5 B 1g orbitals, while the resonances located in the visible region correspond to transitions involving 5 A 1g states.
We turn now our attention to the gyrotropic ς gy spectra. As observed in Fig. 7, their structure is much simpler, with a main resonance located in the visible or nearinfrared, depending on the geometry and value of ∆. The contributions to this resonance come mainly from transitions between 3 E g and 3 A 2g states that are perturbed by spin-orbit coupling (see section Sec. IV). The spectral weight of these transitions is too small to be observed in the nongyrotropic spectra. Nonetheless, their effect on the gyrotropic response is crucial, via the orbital mixing induced by spin orbit coupling. One way to evaluate this mixing is by varying the energy gap ∆ defined above. In particular, we observe in Fig. 7 that the gyrotropic signal is the smallest for the largest value of ∆, as expected from the smaller orbital mixing in this case.
Let us now discuss the effect of orbital hybridization induced by the coupling to the electromagnetic field. As discussed in Sec. IV, the coupling to light induces p − d hybridization, which causes an evolution of the eigenstates as a function of the overlapping between oxygen and transition metal states, as sketched in Fig. 5. In our numerical calculations we observe that for values t pd ≈ 1.2eV and ∆ ⪅ 0.89eV there is a strong orbital mixing. Indeed, for this choice of values, the calculated spectra are in agreement with the experimental spectra reported in Ref. [25]. In particular, in the calculated spectra we observe two nongyrotropic resonances in the red and blue parts of the visible range, respectively, while a main gyrotropic resonance is seen in the blue region, in agreement with the experiments [25].
Next we discuss the effects of the geometry on the gyrotropic response. The data shown in Fig. 7 and Fig. 8 reveals the strong dependence of the spectral functions on geometric factors. The effect is particularly evident for spectra shown in panels (c) of both Figures, which show that the gyrotropic signal is completely extinguished for this particular geometry. As discussed in Sec. IV, the reason for this extinction is that the matrix elements of the angular momentum that connect the t 2g orbitals in adjacent sites are null in this case, because the direction of the momentum component is not contained in the space symmetries of the t 2g orbitals involved in the transfer. On the other hand, the gyrotropic signal for the geometry sketched in panels (b) of both Figures is significantly larger than for the spectra and maps shown in panels (a). The reason is that for the geometry of spectra and mappings shown in Fig. 7(a) and Fig. 8(a), the photoinduced transfer that contributes to the gyrotropic response happens only between adjacent η − τ orbitals. In contrast, for the geometry studied in Fig. 7(b) and Fig. 8(b) both ζ − η and ζ − τ hopping channels are allowed, increasing the gyrotropic signal.
Finally, the data shown in the maps of Fig. 8 reveals the dependence of the gyrotropic response on spinorbit coupling. Indeed, when the latter tends to zero, the gyrotropic signal becomes vanishingly small, while it becomes progressively more intense as the spin-orbit coupling increases. On the other hand, we see that orthorhombic perturbations barely modify the gyrotropic spectra, as expected from the fact that these perturbations are much smaller than the intersite hopping induced by coupling to the electromagnetic field. Both observations are in agreement with the conclusions drawn from the perturbative analysis discussed in Sec. IV.

VII. CONCLUSIONS AND PERSPECTIVES
We used a group theoretical approach to study the interaction of transition metals with electromagnetic fields. For that purpose, we described the relevant electronic states by irreducible representations of pertinent pointgroup symmetries, which were constructed from manyelectron wavefunctions based on Slater determinants. The energetics of the problem was established to comply with Tanabe-Sugano diagrams corresponding to the particular ion under study [1]. Starting from an initial O h symmetry, we analyzed the effect of symmetry reduction due to Jahn-Teller interactions and spin-orbit coupling. The interaction with electromagnetic fields was assumed to produce photoinduced transfer of electrons across the lattice. The model predicts an electronic response that depends on the handedness of circular polarization of light. Key ingredients to this gyrotropic response are spin-orbit coupling and intraatomic e g − t 2g orbital mixing. Remarkably, the gyrotropic resonances are related to photoexcitations where one of the spins is inverted, enabling the use of electromagnetic fields to manipulate spins. We also analyzed the dependence of the gyrotropic response on the relative orientation of Jahn-Teller distortions, light propagation and spin quantization. In particular, we found specific conditions for which the gyrotropic response is largely reduced or even extinguished. We focused our analysis on 3d ions, with the aim of understanding our recent observation of a gyrotropic response associated with Jahn-Teller polarons in La 2/3 Ca 1/3 M nO 3 [25]. Using standard values for Jahn-Teller constants, spin-orbit coupling and charge transfer our model replicates a gyrotropic response in the blue region of the visible range, in agreement with the experiments [25].
Beyond this particular case, a similar approach may be generalized to study the interaction with electromagnetic fields of transition metals with arbitrary point symmetries and spin-orbit couplings. One perspective is the entanglement of spin and orbital degrees of freedom using light at optical wavelengths. One could think, for instance, of studying quantum tunneling of E ⊗ e or T ⊗ e Jahn-Teller vibronic states [37]. In particular, the interaction with light may drive photoexcited states, whose ground state is formed by coherent superpositions of those vibronic states through quantum tunneling, which may be detected with polarized light. Such excitations could form the basis for quantum states of interest for quantum technologies [22,67]. Another prospect may be the study of 4d-5d transition metals hosting (quantum) spin liquids [18][19][20][21]. In this case, magnetic interactions would compete with the coupling to the electromagnetic field, which could lead to a rich diagram of quantum phases as a function of the wavelength of the electromagnetic radiation. To solve this problem, a group-theoretical approach would require working in appropriate regions of Tanabe-Sugano diagrams [1,49]. For instance, for heavy metal d 4 ions the ground state is 3 T 1g instead of 5 E g , due to the larger crystal field that leads to the condition (Dq/B) c > 2.7, see Figure 1a. Additionally, since in this case Jahn-Teller interactions in the t 2g manifold are typically smaller than spin-orbit coupling [20,37,44], the group-theoretical analysis should consider lowering the point symmetry by spin-orbit interactions in the first place.
We note that our model Hamiltonian considers electrons that are subject to Jahn-Teller instabilities. In general, these may coexist with electrons in delocalized bands. This is indeed the case of many oxides, including La 2/3 Ca 1/3 M nO 3 , where both Jahn-Teller polarons and delocalized electrons participate in transport [68][69][70]. We stress that our model captures the essential physics of electrons that are affected by Jahn-Teller interactions, neglecting contributions from delocalized electrons. As we demonstrate here, this is enough to describe the specific contribution of Jahn-Teller vibronic states to the gyrotropic response, which is experimentally distinguish-able from the conventional response arising from delocalized bands [25]. On the other hand, although the analysis based on spectral functions gives a fundamental understanding of light-matter interactions in these solids, further developments can address linear response theory to obtain responses like optical conductivity and permittivity that can be matched with experiments.
We also note that our assumption of photoinduced electron transfer implies an enhanced conductivity at resonant frequencies, which, eventually could be tested experimentally by measuring electronic transport under illumination at relevant wavelengths. These experiments could be done in La 2/3 Ca 1/3 M nO 3 , but other candidates would also comprise materials like (P r x La 1−x ) 2/3 Ca 1/3 M nO 3 [71] or magnetite [72], where optical signatures of Jahn-Teller polarons have been observed [73]. Generally, materials prone to Jahn-Teller instabilities, including colossal magnetoresistance manganites [68][69][70], could be worth exploring in search for gyrotropic responses arising from spin-orbital mixing. On the theoretical side, other models can extend the analysis to the optical responses of clusters of Jahn-Teller ions rather than isolated ions. Alternative models may also explore these responses in the absence of phototransfer, e.g., the photoexcitation of Jahn-Teller states in molecules [37,74], in which the group-theoretical approach should be applied at the level of molecular orbitals. Finally, in the present model, the electromagnetic radiation and the lattice modes are treated as classical fields. Further extensions would require a full quantum approach to describe these fields, especially relevant for the application of the aforementioned ideas to concepts like cavity quantum electrodynamics [67,75].
Appendix A: Wavefunctions of many-electron states

Notation
We use group theory to construct many-electron wavefunctions. We work in a product basis between orbital and spin momenta, so that wavefunctions are defined by kets like: Here S is the spin magnitude, so that 2S + 1 is the spin degeneracy, and Γ is the irreducible representation in orbital space of group G, expressed in the basis γ (which is omitted for unidimensional representations). Finally M is referred to the spin quantum number, M = −S, ..., S.
In the description of spin-orbit coupling, wavefunctions are denoted by |Γγ⟩ since orbital and spin angular momenta are coupled and double group representations are best suited to take into account the electron spin. The wavefunctions describe multi-electronic states, so they are linear combinations of Slater determinants: Here ψ i = ϕ i χ i is the i-th mono-electronic wavefunction, where the orbital angular momentum part is described by ϕ i = ζ, η, τ, u, v (labels are defined in EqEq. (A3)) and χ i = α, β is the spinor part. Beingn the quantization axis for the spin, thenn ·ŝα = +1/2 andn ·ŝβ = −1/2. For brevity, the notation inside the Slater determinants is written as ϕα → ϕ and ϕβ →φ, which is taken from Ref. [1]. Finally, the sign in bras and kets for wavefunctions is denoted by a breve symbol, i.e.,M = −M .

Construction of the wavefunctions a. O h point symmetry
In the case of one electron in a d shell in a cubic crystal field the ten-fold degeneracy of the free ion is broken into a e g shell (with four-fold degeneracy) and a t 2g shell (six-fold degeneracy). The basis angular functions of these two shells can be expressed as linear combinations of spherical harmonics Y m l to get real spatial symmetries: u and v are the basis for the e g shell and ζ, η and τ are the basis for the t 2g shell.
When there is more than one electron in the d shell we construct many-electron wavefunctions using Slater determinants. Since the spin-orbit interaction is small compared to exchange interactions, the many-electron wavefunctions are built by coupling separately the orbital and spin momenta, following the Russell-Saunders coupling scheme and using Clebsh-Gordan coefficients. For the orbital part, the coefficients are adapted for the O h point-group symmetry. The many-electron wavefunctions are expressed as linear combinations of Slater determinants, which take all possible permutations of electrons sitting on the different monoelectronic orbitals and having all possible spin orientations [1].
In the d 3 configuration the ground state for all values of the Racah parameter ∆/B in the Sugano-Tanabe diagrams is 4 A 2g . It can be shown that the wavefunction that describes this four-fold degenerated term is built from Slater determinants with the three monoelectronic orbitals of t 2g as follows [1]: When a fourth electron is added to the t 2g shell, the many-electron wavefunction transforms as a T 1g representation with total spin S = 1 and is built from nine degenerate Slater determinants as follows: The orbital basis functions for T 1g are κ, µ, ν which have the same relation under rotations as the Cartesian coordinates x, y, z but being even under parity. Eq. (A5) give us three determinants. The remaining six are obtained for the µ and ν basis of the 3 T 1g representation. The latter are obtained from κ by rotating in the orbital space by ±2π/3 around the [111] axis.
On the other hand, if the fourth electron is in the e g shell the many-electron wavefunction corresponds to an E g representation. In this case, the spin number is S = 2 and the many-electron wavefunctions can be expressed as ten linear combinations of determinants as follows: Here the bases for the irreducible representation 5 E g are γ, γ ′ = u, v being γ ̸ = γ ′ and the positive (+) signs corresponding to γ = u and the negative (−) signs corresponding to γ = v.
b. D 4h point symmetry When the cell is tetragonally distorted by a Jahn-Teller instability (reducing the symmetry to D 4h ), the T 1g representation is broken into a representation A 2g with orbital symmetry γ = ν, and E g with symmetries γ = κ, µ. On the other hand, the reduction to D 4h symmetry splits the E g representation into A 1g with basis γ = u and B 1g with basis γ = v. Therefore, in D 4h point symmetry the terms are split in the following representations (we show only the terms with maximum spin quantum number): We note that under a tetragonal elongated distortion the term 5 B 1g is lowest in energy. Although in the mono-electronic picture the fourth electron that drives the Jahn-Teller instability sits in an orbital with u symmetry, the many-electron wave-function of the 3 B 1g term has v symmetry. We also note that we introduce a global phase −1 in the term 5 B 1g to eliminate a minus sign.
c. Spin-orbit coupling Finally, we analyze how spin-orbit coupling splits the representations expressed in Eq. (A7). We remind that under spin-orbit coupling, wavefunctions are expressed as double group representations |Γγ⟩. In Eq. (A8) and Eq. (A9), double group representations (on the left side) are expressed in terms of the irreducible representations in D 4h point symmetry (right side). First of all, note that the only term that splits under spin-orbit coupling is 3 E g . This can be understood by observing the reduced matrices in Eq. (7a) of the main text: only the reduced matrix V A2g has one non-zero diagonal element corresponding to this term. The spin part of 3 E g cannot be described by a j = 1 representation in continuous rotation group because spins interact with the orbital space, which, in this case, is described by the E g representation in D 4h . In this symmetry, the continuous spin rotation group D (S=1) decomposes into A 2g + E g representations -which are gerade, since spinors are even under parity inversion [76]-, with representations in spherical basis A 2g with q = 0, and E g with q = ±1. Therefore, we need to obtain the representation of the composite product (A 2g + E g ) ⊗ E g = A 2g ⊗ E g + E g ⊗ E g . Thus, for A 2g ⊗ E g , i.e., when the orbital component E g couples to the A 2g representation in the spin space, it generates two functions that transform as E g representations: On the other hand, the E g representation of the orbital part combines with the E g representation of the spin part, giving the following irreducible representations E g ⊗ E g = A 1g + A 2g + B 1g + B 2g , which are expressed as follows.
Note that we used the Clebsch-Gordan coefficients displayed in Table II in Appendix C to obtain the expressions in Eq. (A8) and Eq. (A9). On the other hand, for some calculations it may be convenient to express Eq. (A8) and Eq. (A9) in spherical basis for the orbital angular momentum of the 3 E g term, which contains components with quantum numbers M L = ±1 labeled as t ± , and t 0 . We can then rewrite the corresponding terms as: Finally, we discuss how spin-orbit coupling splits the representation 3 E g in D 4h symmetry, see (Fig. 1b). For that purpose, we use the Wigner-Eckart theorem to compute the eigenenergies of the spin-orbit matrix elements. First, since ⟨E g κ|E g γE g κ⟩ = ⟨E g µ|E g γE g µ⟩ = 0 for γ = κ, µ, ν, the spin-orbit eigenenergies of the E g terms (see Eq. (A8)) are 0. On the other hand, by virtue of the expressions in Eq. (A11), the terms in Eq. (A9) can be expressed in the following way: where c γq can be obtained from Eq. (A11). The Clebsh-Gordan coefficients necessary to apply Wigner- 2, so that the matrix elements can be computed as: As a result, spin-orbit coupling does not change the energy of the doubly degenerated E g spin-orbit term, while symmetric (A 1g , A 2g ) and antisymmetric (B 1g , B 2g ) representations split by ∓ξ SO /2 with respect to the E g term (Fig. 1b). We also note that the accidental degeneracy of (A 1g , A 2g ) and (B 1g , B 2g ) terms may be eventually lifted if one considers developments beyond first-order relativistic contributions.

Appendix B: Jahn-Teller Hamiltonian
Atoms or ions in a molecule or a unit cell have a position where the energy of the system is minimized. Sufficiently small deviations from these equilibrium positions can be described through a force constant: Here QΓ are the vibronic coordinates that transform under irreducible representationΓ, which can be described in the frame of group theory being linear combinations of the displacements of the atoms in Cartesian coordinates ∆X n , ∆Y n , ∆Z n , and E Γ is the energy of the system which depends on the irreducible representation Γ of the electronic wavefunction. In the presence of orbital degeneracy, the equilibrium positions change spontaneously, reducing the symmetry through the Jahn-Teller theorem [6]. This situation can be described by the addition of another term in the energy of the system that includes the potential energy of the nuclei in the field of the electrons in the state defined by the representation Γ and basis γ, i.e., the adiabatic potential energy surface (APES) ε Γ γ ( ⃗ Q) [9]: where ε Γ γ ( ⃗ Q) is obtained by solving the secular equation for the vibronic coupling matrix operator W which, to second order, is defined as: where V refers to the electron-ion interaction potential. We can thus define first-order Eq. (B4a) and secondorder Eq. (B4b) vibronic coupling terms as follows [9]: These operators transform as the representation of the group corresponding to the lattice distortions [9], so for the computation of the matrix elements we can use the Wigner-Eckart theorem.
For the E ⊗ e problem, the following matrix elements can be derived using the functions defined in Eq. (A3): By using the Wigner-Eckart theorem we can develop the corresponding Hamiltonian as: where υ i are the Pauli matrices in the pseudospin space of {v, u} and the vibronic coordinates have been normalized as Q 2 = ρ sin ϑ and Q 3 = ρ cos ϑ. The eigenstates of this Hamiltonian are: with the following eigenstates: The energy minima are found when ϑ = 2nπ/3, n = 0, 1, 2 -which corresponds to tetragonal elongations along the three axes-and ρ = F E /(K E − 2G E ). In order to simplify the computation, the radial variable is normalized to ρ = 1. Around the tetragonal elongations, the parameter Ω in Eq. (B9), which is defined as: can be approximated as Ω ≈ ϑ. In this situation, w − has a u-like symmetry and w + has a v-like symmetry.
Since we consider tetragonal elongations with orthorhombic perturbations (ϑ = 2nπ/3 ± δϑ), we can work with the following basis: With this basis, and taking into account δϑ, the expression Eq. (B7) is transformed as: As argued in the main text, in the t 2g shell we only consider the T ⊗ e problem. For this problem we also neglect second order vibronic constants. The T ⊗e vibronic constant is defined as: Using again the Wigner-Eckart theorem, the Hamiltonian of the T ⊗ e Jahn-Teller interaction is derived as a function of Gell-Mann matrices, λ k , in the basis {ζ, η, τ }: The energy minima correspond again to the tetragonal elongations with ρ = F T /K T . Since the nuclei motion is much slower than the electronic transitions, we can make the assumption that these minima are the same as the ones for the E ⊗ e problem, so we normalize again ρ = 1, so that F T = K T . These minima correspond to the basis defined before. We can generalize the expression Eq. (B14) to the local basis at each value ϑ n , denoted as {ζ,η,τ }, which is defined by rotations around the [111] axis in the orbital space, i.e.,τ =R n 3 (xyz)τ , whereR n k (x...) defines a rotation of k-th order (angle 2π/k) executed n times along the axis defined by the coordinates in the parentheses. Then, the orthorhombic distortions in T ⊗ e are described by: The description so far is done with monoelectronic orbitals. We can generalize these results to many-electron wavefunctions. The vibronic constants are calculated using a one-body potential, so that for the non-diagonal matrix elements we only need to check the orbitals that are different (see also Appendix D). All the wavefunctions described by Eq. (A5) and Eq. (A6) are defined in the orbital part in terms of determinants of the type |ζητ γ|. Since there is just one different orbital in each Slater determinant, the off-diagonal elements are not modified. Then, for the first order vibronic constants, since the sums for ⟨t|X Eg u |t⟩ and ⟨t|X Eg v |t⟩ for t = ζ, η, τ are null, the results for many-electron wavefunctions are the same as for the monoelectronic orbitals. We note that the same arguments apply for the T ⊗ e problem, since second order vibronic constants are neglected.
Appendix C: Reduced matrix elements of the spin-orbit coupling operator For the computation of the matrix elements of the spinorbit coupling, we use the Wigner-Eckart theorem applied to the spin-orbit operator V Λ λq defined in Sec. II B, see also Ref. [1]. This operator transforms according to irreducible representations Λ in the orbital space with basis λ and S 1 q corresponds to irreducible representations in the spin-rotation group. To calculate a given matrix element, we apply the Wigner-Eckart theorem as follows: where |ΓγSM ⟩ and |Γ ′ γ ′ S ′ M ′ ⟩ correspond to wavefunctions that transform as irreducible representations Γ, Γ ′ in bases γ, γ ′ with spin S, S ′ and spin quantum numbers M , M ′ , while g Λ and g Γ are the dimensionality of representations Λ and Γ. The application of the Wigner-Eckart theorem requires the computation of the reduced matrices ⟨ΓS||V Λ ||Γ ′ S ′ ⟩. The latter have to be hermitic, which, as will be shown below, is guaranteed by the following expression [1]: where, since ⃗ S is expressed in spherical coordinates, we have q = +1, 0, −1. One has to consider also the fol- lowing relations between the Clebsh-Gordan coefficients: The factors ϵ(ΓΛΓ ′ ) and χ(ΓΛΓ ′ ) depend on the phase convention [50]. We have fixed this convention by imposing ⟨Γγ|A 1g uΓγ⟩ = 1 for any representation Γ and basis γ. Assuming this convention, we have: One can verify that Eq. (C3) and Eq. (C4) imply that ⟨ΓS||V Λ ||Γ ′ S ′ ⟩ = ⟨Γ ′ S ′ ||V Λ ||ΓS⟩, which, as mentioned above, guarantees the expected hermiticity of the spinorbit operator.
As shown in Sec. II B, the reduced matrices for the spin-orbit operator are expressed through irreducible representations V A2g and V Eg . To derive the matrix elements of these matrices, we use the Clebsh-Gordan coefficients expressed in Table II and the ladder operators defined as: which are related to the spherical components of the angular momentum operators through: We first derive the matrix elements corresponding to the representation V A2g .
We remind that the reduced matrices are expressed in the basis { 3 E g , 3 A 2g , 5 A 1g , 5 B 1g } (see discussion in Sec. II B). We first note that the direct product A 2g ⊗ A 2g = A 1g implies that the spin-orbit operator in representation A 2g has nonzero matrix elements connecting 5 A 1g and 3 A 2g . Let us find such elements by applying the operators of angular and spin momenta to the wavefunctions of representation 3 A 2g . We choose a spin-orbit operator in representation q = 1 (corresponding to spin operator s + ) and | 3 A 2g ν1⟩ expressed in terms of the corresponding Slater determinants (Eq. (A7a)). By applying the operators directly to the wavefunctions it follows that: On the other hand, by applying the Wigner-Eckart theorem (Eq. (C1)), we obtain: Combining Eq. (C7) and Eq. (C8) we obtain the reduced matrix element ⟨ 5 A 1g ||V A2g || 3 A 2g ⟩ = ı √ 10. Next, we note that the direct product E g ⊗ A 2g = E g implies that E g wavefunctions can be connected through the spin-orbit operator only to wavefunctions of the same representation. This gives a diagonal element in the V A2g matrix. To find such element we apply, as before, the operators of angular and spin momenta to wavefunctions 3 E g (Eq. (A7b) and Eq. (A7c)) and choose a representation q = 0, involving the spin operator s 0 . The application of the operators to the wavefunctions gives: On the other hand, the application of the Wigner-Eckart theorem gives: We therefore obtain ⟨ 3 E g ||V A2g || 3 E g ⟩ = ı √ 3. Using similar arguments, it can be shown that the rest of matrix elements of V A2g are zero, resulting in the reduced matrix described by Eq. (7a).
We derive now the matrix elements corresponding to the representation V Eg . We first note that the direct product E g ⊗ E g = A 1g ⊕ A 2g ⊕ B 1g ⊕ B 2g means that the term 3 E g can be connected by the spin-orbit operator to all other representations. As done for V A2g , we combine the application of the operators of angular and spin momenta to the pertinent wavefunctions with the application of the Wigner-Eckart theorem. We obtain the following expressions: which allow us to obtain all the matrix elements for the V Eg matrix as follows: 15. These elements give the reduced matrix V Eg expressed in Eq. (7b).

Appendix D: One-body operators
In this work, Jahn-Teller and spin-orbit Hamiltonians contain one-body operators. In the same way, lightinduced transfer requires also one-body operators in the electromagnetic Hamiltonian. In the following, we explain how one-body operators act in the formalism of the many-electron wavefunctions defined in Appendix A.
For that purpose, we recall that to comply with Pauli exclusion principle we need to define a multielectronic wavefunction Ψ through the antisymmetrization operator A acting on the product of monoelectronic states occupied by electrons: where σ is an element in the permutation group P, N is the number of fermions of the system and in (−1) σ represents the parity of the permutation. This results in the formation of Slater determinants.
We can define a one-body operator O as the sum of operators o k acting over the k-th fermion as follows: To find the matrix elements, we take into account the following properties of the antisymetrization operator: • Applying A to a Slater determinant returns the same Slater determinant, so that A 2 = A.
• Since A is a real operator A † = A.
• Since any one-body operator is even under permutations it always commute with A, [O, A] = 0.
Consequently, the matrix elements involving one-body operators between many-electron wavefunctions can be found as follows: By orthogonality, the only permutation that does not vanish is the identity. Note that diagonal elements (Φ = Ψ) do not vanish: On the other hand, off-diagonal elements (Φ ̸ = Ψ) can be nonzero only if the many-electron wavefunctions differ only by one one-body wavefunction. Otherwise, if they differ by more than one one-body wavefunction, the inner product vanishes.