Optical control of competing exchange interactions and coherent spin-charge coupling in two-orbital Mott insulators

In order to have a better understanding of ultrafast electrical control of exchange interactions in multi-orbital systems, we study a two-orbital Hubbard model at half filling under the action of a time-periodic electric field. Using suitable projection operators and a generalized time-dependent canonical transformation, we derive an effective Hamiltonian which describes two different regimes. First, for a wide range of non-resonant frequencies, we find a change of the bilinear Heisenberg exchange $J_{\textrm{ex}}$ that is analogous to the single-orbital case. Moreover we demonstrate that also the additional biquadratic exchange interaction $B_{\textrm{ex}}$ can be enhanced, reduced and even change sign depending on the electric field. Second, for special driving frequencies, we demonstrate a novel spin-charge coupling phenomenon enabling coherent transfer between spin and charge degrees of freedom of doubly ionized states. These results are confirmed by an exact time-evolution of the full two-orbital Mott-Hubbard Hamiltonian.


Introduction
The exchange interaction J ex between microscopic spins is the strongest interaction in magnetic systems. Therefore, the control of exchange is a very promising way for ultrafast control of magnetic order, with potentially high energy efficiency. Recently, the ultrafast control of J ex has received significant interest both in experiments with cold atoms as well as in condensed matter systems [1][2][3][4][5][6][7][8][9][10][11]. An appealing way to achieve a control of J ex is to use periodic driving with off-resonant pulses as was extensively investigated theoretically [12][13][14][15][16][17]. In particular, it was predicted theoretically [12,13] and recently confirmed experimentally [11] that by tuning the strength and frequency of the driving, J ex can be reduced, enhanced and even reverse sign in a reversible way. However, so far, most theoretical studies rely on single-orbital models, while multi-orbital physics is important in many materials. Moreover, existing studies [18][19][20][21] on multi-orbital systems did not reveal the role of orbital dynamics on the control of exchange interactions.
In order to have a better understanding of the influence of orbital dynamics on the ultrafast and reversible control of exchange, we report the study of a two-orbital system at half filling under the effect of a periodic electric field. There are two main differences between single and multi-orbital systems which are already captured in the two-orbital case. First, there is the Hund interaction J H that directly arises from inter-orbital exchange on the same site. At half filling and for J H >0, each orbital is singly occupied and the low-energy degrees of freedom are spin-one states which interact both via a normal Heisenberg exchange J ex S i · S j and with a biquadratic exchange interaction B ex ( S i · S j ) 2 . While J ex favors collinear spin order at neighboring sites, for B ex >0 , non-collinear spin order can become preferential. For calssical spins, the presence of biquadratic exchange interaction can lead to spin spiral states [22]. For quantum spins in low dimentions systems, the presence of B ex can give rise to disordered phases such as dimerized or quadruolar phase [23][24][25]. Second, as illustrated in Figure 1, the two-orbital model has excited states which are doubly ionized and strongly gapped with respect to states with only one electron in each orbital (we will refer to configurations with one electron in each orbital as singly occupied states). The doubly ionized states are charge states, which are coupled to singly occupied states by two subsequent hopping processes.
Below we demonstrate that there exist two distinct regimes for the non-resonantly driven two-orbital model. First, a regime for which the control of intersite exchange interactions dominates. We recover a Heisenberg exchange interaction J ex (E, ω) that is similar as in singleorbital systems, where E is the driving strength and ω the driving frequency. We find that analogous to J ex (E, ω), also the biquadratic exchange interaction B ex (E, ω) can be reduced, enhanced and reverse sign by tuning the strength and frequency of the driving field. In addition, we find a regime for which the exchange interactions compete i.e. J ex (E, ω)∼B ex (E, ω)>0. Second, we elucidate a regime for which a new type of spin-charge coupling phenomenon dominates over the exchange interaction. In this regime, a reversible transfer between spin and charge degrees of freedom is feasible. The paper is organized as follows: in Section 2 we introduce the two-orbital Hubbard model, define projection operators, and introduce a generalization of the time-dependent canonical transformation [16,18,26,27]. In Section 3 we derive the effective Hamiltonian, study its low energy part, and show how to map it onto a spin-one model. From this spin-one model, the Heisenberg exchange interaction as well as the additional biquadratic exchange interaction are extracted. Beyond the spin model, we study the spin-charge coupling phenomenon. Moreover, we confirm the analytical results on the spin-charge coupling by computing the time evolution of the full two-orbital Mott-Hubbard model for a two-site cluster. Finally, in Section 5, we draw conclusions.

Electronic model
To study the role of orbital dynamics on the electrical control of exchange, we investigate a two-orbital model at half-filling. This can be associated with the e g band of an oxide compound. The Hamiltonian is given byĤ(t) =Ĥ U +Ĥ kin (t), whereĤ U =Ĥ nn +Ĥ sf .Ĥ nn , H sf andĤ kin contain the density-density interaction, the spin-flip and pair hopping, and the intersite hopping, respectively: Hereĉ † iασ (ĉ iασ ) are fermionic creation (anihilation) operators for site i, orbital α=a, b, spin σ= ↑, ↓, andn iασ =ĉ † iασĉ iασ . U is the on-site Coulomb interaction and J H is the Hund exchange interaction.
The time-dependence of the hopping term originates from the external electric field which is described using the Peierls substitution t ij (t)=t 0 e ieA ij (t) [12,28,29], where e is the electronic is the projection of the vector potential along the direction from site i to j, where E 0 is the amplitude of the field. Since both e g orbitals originate from d orbitals, no on-site electric dipole transition are allowed. We define the parameter E = eaE 0 /ω which represents the driving strength, whith a=|R i −R j | and we take t 0 =1 for the numerical calculations.

Projection operators
The conventional way to derive the exchange interaction is to use a canonical transformation also known as Schrieffer-Wolff transformation [16,18,26,27]. For the two-orbital case, this is more involved due to the Hund interaction J H . To deal with this additional complexity, we restrict the Hilbert space to blocks involving only two sites (ij). For all states |φ k on the bond (ij), we then define projection operatorsP ν d (N, M ) onto the following quantum numbers: • Particle number: whereN = iασn iασ and N =0, ..., 8 the number of electrons which occupy the system.
It is shown in Appendix A that explicit expressions forP ν d (N, M ) in terms of single-electron operators can be derived usingp whereĥ iασ =(1−n iασ ). With these definitions, the identity reads The hopping term Eq. (3) connectsP ν d with different d and can be re-written in terms of operatorsT +1 (t),T −1 (t) andT 0 (t) that change d by +1, −1 and 0 respectivelŷ whereT andT Expressions for forT +1 (t),T −1 (t) andT 0 (t) in terms of single electron operators are given in Appendix A. The projection operatorsP ν d and hopping operatorsT +1 (t),T −1 (t),T 0 (t) play an important role in the canonical transformation described below.

Generalized time-dependent canonical transformation
The canonical transformation is a technique which enables the derivation of an effective Hamiltonian for the subspace of statesP ν d [16,18,26,27,30]. Formally, this is achieved by unitary transformationV (t)=e −iŜ(t) that transforms the HamiltonianĤ(t) to a rotated frame. The effective Hamiltonian in the rotated frame readŝ The aim is to identify a suitable subspace (defined by values of d and ν) and determineV such thatĤ eff leaves this subspace invariant. To do this, we perform the unitary transformation perturbatively, treating the hopping parameter t 0 U as a perturbation. We expand iŜ(t) andĤ eff (t) in terms of a Taylor series For deriving a pure spin model, one could construct the unitary transformation such thatĤ (n) eff does not contain terms that change d [16,[30][31][32], and obtain an effective Hamiltonian in the subspace d = 0. Here we use a more general requirement which will allow us to derive an effective Hamiltonian in a subspace different from that without doublons. This turns out to be crucial for a description of multi-orbital systems. We enlarge our effective model and keep terms that change d by ±2, while we designP ν d iŜ (n)P ν d such thatP At half filling and without inter-orbital hopping (t α =β =0), only odd orders of iŜ (n) (t) ∝ t n 0 remain, Eqs. (17) and (18) not only allow us to obtain an effective description of the low energy stateŝ P ν 0 , but also enable us to keep track of the coupling between the low energy space (spin:P ν 0 ) and the space with the highest excited states (charge:P ν 2 ). Eqs. (17) and (18) yields the zeroth order contribution to the effective Hamiltonian Using the projection operators, we obtain the following equation for iŜ (1) In contrast to the zeroth order contributionĤ  [18]. Here we use a simpler algebraic solution that is feasible for time periodic driving and which is closely related to Floquet theory [12,16,31] and the high frequency expansion [13,26,32]. Given a time periodic electric field E(t)=E(t + T ) with a period T = 2π ω , we can expandT ±1 (t) and iŜ (n) (t) in a Fourier series as followŝ where m is the Fourier index, which can be seen as the number of virtual photons absorbed by the system [31]. Using Eqs. (20) and (21), we obtain: For ν = 1, E ν d is a matrix and we would have to further decompose P ν d for the procedure to be exact (see also Section 2.2). Here instead we use an approximation This is a generalization of the energy approximation employed in [33], where E d is approximated by the mean energy of all states for given d. The present approximation is accurate for where the number of doublons d =d and the Hund rule violation index ν =ν . We find that this condition is satisfied for the calculations presented in Section 3. The first order effective HamiltonianĤ (1) eff (t) vanishes becauseT 0 (t)=0 for orbital-diagonal hopping t α =β =0. Eq. (22) allows us to compute higher order contributions toĤ eff (t) in a straightforward way. The second order contribution readŝ where, d, d =0, 2 and ν, ν =0, 1. The third order contribution toĤ eff (t) gives us an expression forP ν d iŜ This yields the following fourth order contribution to the effective Hamiltonian: With Eqs. (25) and (27) we have derived the central result of this section, namely an effective Hamiltonian up to fourth order in the hopping.

Results
In this section we present the results obtained with the projection operators and the effective Hamiltonian derived above. First we show that the d = 0 part can be mapped onto an effective spin-one (S=1) model. This requires two additional unitary transformations: one to reduceĤ eff to the d=0 sector and a second to specialize to the S=1 states only. In the effective spin-one Hamiltonian, we extract the Heisenberg J ex (E, ω) and biquadratic B ex (E, ω) exchange interactions. Second, we focus on the coupling terms between sectors d=0 and d=2 by taking both of them into account in the low energy description. This goes beyond the spin model and captures the spin-charge coupling dynamics.

Spin-one model
According to condition Eq. (17), the full effective Hamiltonian yields n=2,4Ĥ (n) In this subsection we study the low energy effective Hamiltonian up to fourth order in the hopping. In the derivation of the spin-one model, we have to consider the sectorP ν 2 as a high energy sector and perform a second time-dependent canonical transformation in order to project out states for which d=2: where ν,ν P ν 0H (4) We used that C ν0,m d0 =C ν1,m d0 . Note that this canonical transformation includes all modes m from the first canonical transformation. Details of the second canonical transformation are given in Appendix B and illustrated in Figure 2b. We would like point out that in the full lattice, additional 4 th order interactions occur, such as ring-exchange terms, spin chirality terms [16] as well as additional 4 th order contribution to the Heisenberg and biquadratic exchange interactions. Since we restrict ourselves to a two site model, such processes are not taken into account in our calculations.
Hamiltonian Eq. (29), can be written in terms of spin-one (S=1) operators as described before in [34,35]. In general, S=1 operators can be defined using many-electron operators [36]. Here, we define their projection onto local spin states |S, M S Then, we can write the spin-one states in terms of single electron states using suitable Clebsh-Gordan coefficients where one can writeŜ q in terms of single electron operators (index i is omitted for brevity), which yieldsŜ Using the definitionŜ ±1 =∓ 1 √ 2 (Ŝ x ±iŜ y ) for the spin-one spin flip terms [36], one can compute the product S i · S j as well as ( S i · S j ) 2 in terms of single electron operators and identify them with the terms of Eq. (29). This procedure leads to an effective Hamiltonain written in terms of S=1,R 1 ij andR 2 ij operators, see Appendix D. Subsequently, by time averaginḡ d=0 eff (t)dt, we obtain an effective time-independent Hamiltonian: K 1 (E, ω) corresponds to the exchange J ex (E, ω) up to second order in the hopping. K 2 (E, ω) gives a fourth order contribution to J ex (E, ω) as well as the biquadratic exchange B ex (E, ω). K 3 (E, ω) gives a contribution directly to B ex (E, ω).
The remaining termsR 1 ij andR 2 ij describe orbital resolved spin dynamics that strictly go beyond a spin-one model. Their expression in terms of fermionic operators can be found in Appendix D. To arrive at an effective spin-one model only, we perform a third time-dependent canonical transformation between the spin one sectorP S and the S =1 sectorP R from theP ν 0 sector. In this process, ilustrated in Figure 2c,P R is taken as a high energy sector. Details of the calculations can be found in Appendix D. Eventually, we obtain an effective spin-one where J m is a Bessel function of order m. The first term of Eq.(38) corresponds to K 1 (E, ω) and the second term is a contribution from K 2 (E, ω).
We now would like to compare the behavior of the second order J ex (E, ω) in single and two-orbital systems. J ex (E, ω) for single-orbital systems reads We can see that J ex (E, ω) in the two-orbital model, Eq. (38), has an additional factor 1/2 as compared to J single ex . This is due to the inter-orbital hopping t αβ =0 which changes the prefactor of J ex (E, ω) as compared to the single-orbital case. For t αβ =t αα =t 0 , the second order contribution of the single and two-orbital model would have the same prefactor. However, the relative modification of the exchange ∆J ex (E, ω)/J ex (E, ω) is the same and therefore, orbital dynamics does not change the control of J ex (E, ω).
The biquadratic exchange interaction can be writen as a sum of six contributions: where with We used the notation B [i] ex (P ν d ) to denote the first, second and third canonical transformation via theP ν d sector, for i=1, 2, 3 repsectively. Since the energy approximation Eq.(24) leads to a same energy for bothP 0 0 andP 1 0 , we group the biquadratic paths via these two sectors into one path via theP 0 sector, B In deriving Eqs. (41)(42)(43)(44)(45)(46), one obtains factors which contain k, l, m and n indices, see Eqs. (47) and (48), these directly arise from the canonical transformation and come from Bessel functions J −m which are symmetric for even m but anti-symmetric for odd m.  Figure 3a shows the typical behavior of J ex (E) while the lower panel shows behavior of B ex (E) for frequencies ω=9, 18 and 25.We observe that J ex (E, ω) can be controlled with the strength E and frequency ω of the electric field similarly as found in [12] for the single-orbital system i.e. J ex (E) can be reduced for frequencies above the Mott gap U + J H , enhanced for frequencies below the gap and reversed for stronger driving field E. The major contribution to J ex (E) comes from the second order contribution in the hopping. Figure 3b shows the contributions of the different biquadratic paths B  ex (P ν 0 ) in red and B [2] ex (P 0 2 ) in blue which are the strongest contributions to the biquadratic exchange. The middle panel displays B [1] ex (P 1 2 ) in red and B [2] ex (P 1 2 ) in blue. On the bottom panel, we plotted the weakest contributions to the biquadratic exchange: B [1] ex (P 0 2 ) in red and B [2] ex (P 0 ) in blue. By summing up all the B  3a. In equilibrium, B ex (E=0)<0 favors a collinear alignment of spins in the classical limit and |B ex (E=0)| is weak as compared to |J ex (E=0)|.
We observe that analogous to J ex (E, ω), also B ex (E, ω) can be controlled by the electric field strength and frequency. In the regime of low driving field strength E 1, |B ex (E, ω)| is reduced for frequencies above the Mott gap, ω=18 and 25, and enhanced for the frequency below the gap, shown here for ω=9. The enhancement of |B ex (E)| can be understood from Figure 3b where |B [1] ex (P 0 2 )| and |B [3] ex (P 0 )| are both enhanced at low driving field. In addition, the sum of B [1] ex (P 1 2 ) and B [2] ex (P 1 2 ) gives a reduction of B ex . Eventualy, the sum these four contributions dominates over the large enhancement of the B [1] ex (P 0 ) contribution leading to an enhancement of |B ex (E)|. The physical mechanism behind the increase/reduction of |B ex (E, ω)| for low driving field strength can be explained as follows: the virtual hopping to mω high energy states is enhanced or reduced as compare to equilibrium. This leads to an increase, decrease or change of sign of the C νν ,m dd products in different biquadratic paths B Eqs. (42)(43)(44)(45)(46). Summing all the biquadratic paths, the total |B ex (E, ω)| is enhanced or reduced as compare to its equilibrium value.
For larger driving field strenght E 1, the reduction of photo assisted hopping as well as the oscillatory origin of the Bessel function can lead to a change of sign of B ex (E, ω). Interestingly, for frequency ω=9, we identify a regime for which both J ex (E) and B ex (E) are positive for a range of driving field strenght E>1, this regime is diplayed with a gray area in Figure 3a. At ω=18, the rectangle box in Figure 3a shows a regime for which J ex (E)∼B ex (E)>0. Within this regime, both J ex (E) and B ex (E) are positive which leads to a competition between the exchange interactions. We note that our results suggest that in principle it is possible to change the ratio of J ex (E)∼B ex (E) over a large range, where the equilibrium phase diagram in 1D shows several distinct quantum phases. It will be very interesting to study the feasibility of dynamical transitions between such phases in future work. Analogously, it might be very interesting to study the emergence of non-collinear order in classical spin systems by perturbation of the ratio J ex (E)∼B ex (E). For resonant photo-excitation, this problem has been studied recently and it was indeed found that the non-collinear phase can emerge [22]. Next, we study the possible enhancement of B ex (E, ω)/J ex (E, ω) for a wider range of freqeuncies (9≤ω≤40) where ω=40 is larger than the highest energy of the undriven system E 0 2 −E 0 0 . The result is shown in Figure 4 as a color map as a function of E, ω. Positive values of the ratio are shown in yellow and negative values are shown in blue. Below the Mott gap (ω=12), accurate results can only be obtained in a frequency range 9≤ω≤9.5. Below and above this range until ω 17, the energy approximation Eq.(24) breaks down and orbital resolved spin dynamics [34] is required to have an accurate description of the exchange interactions. Therefore, the frequency range 10≤ω≤17 is not shown. Figure 4 clearly demonstrates that the exchange ratio can be enhanced as well as reduced depending on ω and E. The parameters for which B ex (E, ω)/J ex (E, ω) is strongly enhanced correspond to three types of situations: • Frequencies for which ω=E ν d −E ν d +mω, this corresponds to white lines at ω=10, 12, 16, 20 and 32. At these frequencies, B ex (E, ω) diverges, such that the canonical transformation breaks down. Note that ω=32 is the frequency that separates spin states from the doubly ionized state, such that a coupling appears close to this frequency. This coupling leads to charge dynamics and therefore, the spin-one model is not accurate in this region. This coupling to charge states is studied in detail in the next section.
• Field parameters ω and E for which B ex (E, ω)∼J ex (E, ω) are indicated in white curved areas. For these regime, the exchange ratio |B ex (E, ω)/J ex (E, ω)| is enhanced since J ex (E, ω) 0.
This leads to a regime where B ex (E, ω)>J ex (E, ω) is realised however, B ex (E, ω) itself remains small as compared to J ex (E=0).
• Parameters ω and E for which the relative sign of B ex (E, ω)/J ex (E, ω) is changed due to the change of sign of B ex (E, ω) leading to a slight enhancement of B ex (E, ω)/J ex (E, ω). This can be clearly seen for frequency ω 9 at E 2.
Summarizing, orbital degrees of freedom do not change the behavior of J ex (E, ω). Both sign and strength of J ex (E, ω) and B ex (E, ω) as well as their relative sign can be controlled by driving, while the regime for which B ex (E, ω)∼J ex (E, ω) is reached only for J ex (E, ω) J ex (E=0).

Beyond the spin-one model
Besides the additional term B ex , the inclusion of orbital degrees of freedom also gives rise to qualitatively new effects that go beyond a description in terms of a spin model alone. In particular, under driving there can be coupling to the doubly ionized charge sector (d=2), which is irrelevant in equilibrium due to the large energy difference between the sectors. Under non-equilibrium conditions, the spin charge coupling readŝ We now study the regime for whichP ν 0 andP ν 2 from different m sectors overlap. This overlap appears for frequencies ω close to E ν 0 =E ν 2 +mω. Although this seems a resonant condition, a direct optical transition is not possible since two hoppings are required to go from theP ν sc,∆m =0 . The first contribution represents the coupling within one Fourier sector m. This contribution remains weak sinceP ν 0 andP ν 2 states are strongly gapped when they belong to the same Fourier sector. We therefore focus on the second term that allows coupling between the spin sectorP ν 0 and the charge sectorP ν 2 with different m. For small E, the leading contribution to Eq. (49) arises from n=2 and ∆m=±1. Here we restrict to the coupling between betweenP ν 0 from m=0 andP ν 2 from m=−1 sector. This yieldŝ To illustrate the spin-charge coupling, we restrict ourselves to the space ν =0. These are two states that have all electrons either on site i or on site j. The full expression ofĤ (2) sc,|∆m|=1 (t) in terms of fermionic operators is given in Appendix E.
To show the effect of this coupling, we compute the low energy spectrum for driving frequencies ω=ω 0 +δω, where ω 0 =|E ν 0 −E 0 2 |. The spectrum is shown in Figure 5 for δω = 0.5 in the two-site system. In this case, the lowest energy state is a singlet state that couples to the two charge states ofP 0 2 . The latter are degenerate up to t 4 0 since four hoppings are needed in order to transfer the four electrons from one atom to the other one. Black lines, from top to bottom, show the quintet state (S=2) and the triplet state (S=1) that are not involved in  (2) eff that contains the spin-charge coupling terms, see Appendix E. The eigen-energies show an avoided crossing, which reveals a hybridization between the spin and charge states. For driving frequencies far away from ω 0 , the spin-charge coupling termsĤ (2) sc,|∆m|=1 (t) remain small and for these frequencies, the spin and charge states are gapped. Therefore, the hybridization is negligible and we recover the regime for which the effective spin model is valid. However, in the regime ω∼ω 0 , the hybridization between the spin and charge states cannot be neglected and the exchange interaction formula obtained in the previous section are no longer accurate.
To sketch the hybridization process, let us take the equilibrium ground state state of the system which is the singlet state (spin state). After switching on the electric field and by slowly changing the field amplitude (∂ t E 0 /E 0 ω), one anticipates that the system starts in a pure spin state and, approaching the avoided-crossing regime, charge states are mixed to the state. For strong E, the spin state becomes a pure charge state.
In the literature, a time-dependent traverse of an avoided crossing is widely studied. The basic example is the Landau-Zener (LZ) effect [38,39]. Also condensed matter systems can exhibit LZ physics. For instance, the Zener breakdown has been studied [40,41] in semiconductors and more recently in Mott insulators [42]. Nonetheless, distinct from these LZ effects which involve changes of the electrical conductivity, here we report coherent transfer of spin to charge degrees of freedom that keep the system in an insulating regime.
Summarizing, orbital dynamics gives access to charge states that are not accessible in equilibrium. The spin-charge coupling offers the possibility to induce coherent charge dynamics in the system. This phenomenon appears in the non-equilibrium low-energy spectrum as an avoided crossing. We stress that it is quite distinct from spin-orbit coupling since here we have an interplay with Coulomb and Hund interaction with the driving field. The spin-charge coupling is not present in single band systems, and we expect it to be universal for multi-orbital systems. Indeed, since multi-orbital systems offer the possibility of having multiple excited states (multiple doublons), multi-doublon excitation should be possible for multi-orbital systems in general. Note that here we did not studyP 1 2 states, they are nonetheless very interesting since they have a lower energy (at and below gap energy U +J H ) and therefore are reachable with lower frequencies ω.

Time-dependent numerical simulations
In the previous section, we showed that the generalized canonical transformation can capture a qualitatively new phenomenon that couples the spin and doubly ionized charge sector. To further support this finding, we perform an exact time propagation of a cluster of two sites described by the Mott-Hubbard Hamiltonian. We focus our attention on the coherent transfer of spin to charge degrees of freedom. In order to describe the charge dynamics, we define pseudo-spin one operatorsT 0 that are composed of Anderson pseudo-spin 1/2 operatorŝ [43,44]. The construction ofT 0 is inspired by the expression of the spin-one operatorŜ 0 . By using a similar procedure withT 0 and pseudo-spin 1/2 operators, we obtainT that characterizes fully occupied and completely empty sites from theP 0 2 sector. Operatorŝ T +1 andT −1 can be defined analogously, see Appendix C.
To solve the time-dependent Schrödinger equation, we use a second order commutator-free approximation of the time-propagator [45] and we compute the time-dependent wavefunction |Ψ(t) and evaluate observables as Ô = Ψ(t))|Ô|Ψ(t) . We focus on three different observables: First, the spin correlation S i · S j to show the spin dynamics during the laser pulse. Second, we characterize the charge statesP 0 2 with T 0 iT 0 j . Finally, to probe the states that possess one doublon d=1, we evaluate N d=1 , whereN d=1 =P 0 1dP 0 1 . Figure 5 shows simulated spin-charge dynamics for an electric field E(t)=E 0 cos(ωt) × exp −(t−t * ) 2 /τ 2 , where E 0 is the amplitude of the field, t * is the time at which E(t) peaks and τ is the pulse width.
We choose a Gaussian envelope with τ =4000π/ω, ω=ω 0 +δω, with δω=0.5, such that τ ×ω sc 1 where ω sc 0.1 is the energy splitting of the avoided crossing (see Figure 5). It has been shown that the effective Hamiltonian picture can break down for long-time dynamics in the thermodynamic limit, because the system heats up to infinite temperature [46]. Here we restrict ourselves to a two-site system to mimic the dynamics of a large system at relatively short timescales. However, for generic systems it is shown that heating can occur at short timescale since the adiabatic limit of Floquet does not exists [47]. Nevertheless, here the use of Floquet restricts to the derivation of an effective Hamiltonian which gives a qualitative picture of the avoided-crossing. We confirm the reversibility of the spin-charge coupling phenomenon within the two-site system with the time-dependent numerical simulations displayed in Figure  6.   Figure 6c shows S i · S j for different E. In equilibrium and for small E, S i · S j −1.9, which slightly deviates from the pure spin case ( S i · S j =−2) due to hybridization withP ν 0 ,P 0 1 andP 0 2 sectors. Figure 6c shows that, with increasing E the state has less spin characteristics. Moreover, it is observed in Figure 6a that with increasing E, the state has more charge characteristics T 0 iT 0 j . In addition, after the laser pulse, both charge excitation and the spin correlations return to their initial value, demonstrating that the coupling is reversible.
Further, we confirm that for frequencies away from ω 0 the spin-charge coupling dynamics is not present. This is shown in Figure 6a,c where the charge and spin dynamics are plotted in red lines for δω=3. In this case, no enhancement is observed. Similarly, in Figure 6c, the spin correlations are not diminished. Moreover, we show the single doublon states dynamics N d=1 in Figure 6b, for a field strength E=0.7. We observe that the laser pulse does not trigger any positive excitation ofP 0 1 states. This means that enhancement of charge dynamics is not due to resonant excitation of the intermediate excited statesP 0 1 . Interestingly, we actually observe depopulation of theP 0 1 states during the laser pulse. This means that moreP 0 1 states are virtually excited to the doubly inonized sector than spin states excited to theP 0 1 sector. Finally, the inset of Figure 6a shows values of the peak of the spin correlations S i · S j t peak for each E in dots and values of the spin correlations S i · S j as a function of E is obtained from the anaytical calculation of Section 3.2. Good agreement between analytical and numerical results is found, which confirms the predictions of the analytical theory. The slight discrepancies between S i · S j E and S i · S j t peak at zero field E can be reduced by taking into account higher order terms in the effective Hamiltonain as well as the spin-charge coupling dynamics to the full P ν 2 sector. In addition, we note that with careful tuning of δω, the reduction of | S i · S j | as a function of E at small E can be made even stronger e.g. by increasing δω, one can move the avoided-crossing closer to E=0 which enhances hybridization with the charge states at small E.

Conclusion
In summary, we obtain an analytical expression for the Heisenberg J ex (E, ω) and biquadratic B ex (E, ω) exchange interaction in a periodically driven two-orbital system. We show that J ex (E, ω) can be controlled analogous to the single-orbital case. We find that for low driving strength, |B ex (E, ω)| can be reduced for frequencies above the Mott gap and enhanced for frequencies below the gap. In addition, we show that B ex (E, ω) can even change sign for stronger driving field srength. We demonstrate that B ex (E, ω)/J ex (E, ω) can be controled by driving while the regime for which B ex (E, ω)∼J ex (E, ω) is reached only for J ex (E, ω) J ex (E = 0). Moreover, a new coupling between spin and charge states is demonstrated. While this coupling is negligible in equilibrium, it can be strongly enhanced and even dominate under driving. This coupling leads to a hybridization between spin and charge states for frequencies close to the spin-charge gap. In contrast to a common charge excitation by resonant photoabsorption, the spin-charge coupling allows non-resonant and reversible coupling to charge degrees of freedom. We have furthermore confirmed these results by simulating the electron dynamics on a two-site cluster.
Natural extension of this work are numerical studies for extended systems. This could be possible for example using multiband extensions of nonequilibrium Dynamical Mean Field Theory (DMFT) [48][49][50][51][52]. We emphasize that, besides the possibility to induce coherent charge dynamics, the presence of the spin-charge coupling should also be visible for short pulses, enabling the excitation of doubly ionized states which could remain coherent due to the gapping with the normal Mott-Hubbard gap. This suggests interesting perspectives for enhancing electronic coherence in correlated electron systems. In this context, it will be interesting to explore the applicability of the two-orbital model to experimental spinone systems such as KNiF 3 [8], which in the literature are considered as prototypical S=1 systems. We would like to point out that modification of the charge occupation is known to systematically influence phonon excitation. Since we did not take into account electronphonon interactions in our model, an outlook would be to study phonon excitation induced by the spin-charge coupling. In addition, interesting prospect of this work is the two-orbitals system under arbitrary fields similarly to what is done in [53] and, since our approach with the canonical transformation can also be applied for arbitrary time-dependent fields [54]. Its extention to more exotic forms of time-dependent fields seems feasible and is left for future work. Finally, we hope that our work can find applications in cold atoms systems, where multi-orbital systems can nowadays be engineered [55,56]. With an adiabatic ramping of the electric field strength, the fully reversible spin-to-charge conversion might be directly observed in double-well systems.

A Projection operators
In this section, we express the projection operatorsP ν d as well as the hopping operatorsT ±1 , T 0 in terms of the electron operatorsĉ † iασ (ĉ iασ ).
From Eqs. (11)(12)(13) of the main text, we can compute the following expressions for T ±1 and T 0 in terms of single electron operatorŝ iασĉ jασĥjασ (n iβσĥjβσ +ĥ iβσnjβσ ) and

B Effective Hamiltonian
In this section we provide explicit expressions for the effective Hamiltonian up to fourth order in the hopping in terms of electron operators. In addition, more details are given for the second canonical transformation that is used to obtain the mapping of the d = 0 sector. The second order contribution to the low energy effective Hamiltonian reads Note that for theP ν 1 sector, only ν=0 contributes since theP 1 1 sector is not connected toP ν d =1 for inter-orbital hopping t α =β =0. The fourth order contribution to the low energy effective Hamiltonian reads: where the first term,P ν 0Ĥ (4) eff (t)P ν 0 , is computed using the direct generalized canonical transformation, Eq.(27) while the second termP ν 0H (4) eff (t)P ν 0 is computed using the second canonical transformation with termsP ν

0Ĥ
(2) eff (t)P ν 2 +h.c. in order to obtain the fourth order contribution to the effective Hamiltonian in the d=0 sector. After developing the commutator in Eq. (27) of the main text and inserting identities whereT ±1 m =T +1 m +T −1 m . We introduce the shorthand notations: We express iŜ where we used C ν0,m d0 =C ν1,m d0 . Note that the last term of Eq. (67) describes the hopping process viaP ν 0 states and therefore, is simpler than the rest of the equation which describes hopping processes viaP 0 2 andP 1 2 . Performing the second canonical transformation we obtain: The effective hoppingsT ±1 m (t) in the second canonical transformation are determined by second order off-diagonal contributions toĤ (2) eff : Substitution of Eq.(25) yields andT −1 Note thatT qq kl =T qq lk , for k =l. Finally,P ν 0H (4) C Spin-one and pseudo spin-one operators Here we derive expressions for pseudo spin-one operatorsT q , where q=±1, 0. Starting from the spin-one operators in term of spin 1/2 operators acting on orbital α, β=a, b.

D Derivation of the effective spin-one Hamiltonian
Here we derive the effective Hamiltonian up to fourth order in the hopping in terms of spinone operators. First we describe the second-order effective Hamiltonian ν,ν P ν

0Ĥ
(2) effP ν 0 in terms of spin-one operators andR 1 ij . Second, we describe the contribution to the fourthorder effective Hamiltonian ν,ν P ν 0 Ĥ (4) eff +H (4) eff P ν 0 obtained with the first and the second canonical transformation in terms of spin-one operators,R 1 ij andR 2 ij . Third, we do the third canonical transformation which takes a state from P ν 0 sector as a high energy state. Since the low-energy subspace contains not only spin-one terms, we perform an additional downfolding.
The bilinear spin-one term reads The first two terms of S i · S j allows an exchange interaction process that transforms a state that does not violate Hund rule (ν=0) to a state which violates Hund rule (ν=1) and vice versa. The last term is a density term which stands forŜ 0 iŜ 0 j . One can derive the following relation and use it to writeP ν

0Ĥ
(2) eff (t)P ν 0 in terms of S i · S j . After time averagingH

0Ĥ
(2) eff (t)dt, we obtain: The termR 1 ij contains a description for the non spin-one states fromP ν 0 , the coupling between a S=1 and a S =1 state and a constant term. All these features can be easily seen after the rotation ofR 1 ij into the spin-one basis. Here we restrict toR 1 ij written in terms of single electron operators: Similarly, we mapP 0H (4) effP 0 onto the spin-one model. The biquadratic spin-one term reads The first summation in Eq.(83) contains terms which are very similar to S i · S j i.e. it contains density terms which describe the states in theP 0 0 sector and terms which connect theP 0 0 states to theP 1 0 states. The second summation contains density terms which describe the states in theP 1 0 and terms wich operate an internal mixing of theP 0 0 states as well as an internal mixing of theP 1 0 states. From Eq.(67) we obtain the following equalities k,l,m,n ν,ν ,ν for the hopping process viaP ν 0 sector, for the hopping process viaP 0 2 sector, k,l,m,n ν,ν for the hopping process viaP 1 2 sector. We obtain the following expression forR 2 R 2 ij rotated into the spin-one basis contains a fourth order contribution to the energy of the non spin-one states, a contribution to the coupling between the spin-one and the non spin-one state and the same constant as inR 1 ij .
We use these equalities to write the time averaged effective Hamiltonian in terms of spinone operators where K 1 (E, ω)=J ex (E, ω) is the second order Heisenberg exchange interaction and reads K 2 (E, ω) is responsible for an energy contribution to both the fourth order Heisenberg exchange as well as the biquadratic exchange interaction, see Eqs. (38) and (40) of the main text.

21
− C 00,m 10 ), We now interest ourselves to the coupling between the spin-one and the non spin-one state of P ν 0 that is described byR 1 ij . The basis transformation which allows one to go from a electron occupation number basis to the angular momentum basis is the following where C SM S S i M i ,S j M j are Clebsch-Gordan coefficients. From this basis transformation, we obtain three spin-one states, namely a singlet (S=0), a triplet (S=1) and a quintet state (S=2) and |2, 0, 1, 1 = 1 √ 6 | ↑, ↑ i | ↓, ↓ j + | ↓, ↓ i | ↑, ↑ j + | ↑, ↓ i | ↑, ↓ j + | ↓, ↑ i | ↓, ↑ j + | ↑, ↓ i | ↓, ↑ j + | ↓, ↑ i | ↑, ↓ j , and three states which are non spin-one states. We define spin-one projection operators such that νP ν 0 =P S +P R , where subscripts S and R refere to S=1 and S =1 states respectively [34]. The third canonical transformation leads to the following fourth order contribution to the effective HamiltonianH where, E S and E R are energies of the S= 1 and S =1 states, respectively. The coupling between theP S andP R subspaces only involves the singlet state |0, 0, 1, 1 and the S =1 state |0, 0, 0, 0 |0, 0, 0, 0 = 1 2 | ↑, ↓ i | ↑, ↓ j + | ↓, ↑ i | ↓, ↑ j − | ↑, ↓ i | ↓, ↑ j − | ↓, ↑ i | ↑, ↓ j , (98) Note that, the energy approximation of Eq. (24) is only for d =d . Here, within the d=0 sector, we take the exact value for the energies of |0, 0, 1, 1 and |0, 0, 0, 0 , 0 and 4J H respectively. After time averaging and projection onto the singlet state, yields This is the additional fourth order energy contribution to the singlet state due to the additional downfolding ofP ν 0 . Using the spin-one Hamiltonian one can obtain a relation between J ex and B ex and the spin-one state energies Since the additional downfolding ofP ν 0 gives rise to an energy contribution for the singlet state only, we have the additional energy contribution to the biquadratic exchange interaction B [3] ex (P 0 )=Ẽ (4) Singlet /6, which leads to B [3] ex (P 0 ) = − k+l+m+n=0 A 1 klmn (E) (C 00,k 01 −C 00,l 10 )(C 00,m 01 −C 00,n 10 ) 4(4J H + (k + l)ω) , where Note that B [3] ex (P 0 )∝J ex (E, ω)/J H , which means that the canonical transformation gives an accurate description of B [3] ex (P 0 ) as long as J ex (E, ω) J H . This inequality is always fulfilled for the regime of frequencies studied here but breaks down when orbital resolved spin dynamics is needed, see the discussion in Section 3.1.

E Effective Hamiltonian with Spin-Charge coupling
In this section we discuss in more detail the effective Hamiltonian which describes the spincharge coupling phenomenon. We study the effective Hamiltonian responsible for the nonequilibrium spin-charge coupling betweenP ν 0 states from Fourier sector m=0 andP 0 2 states from the m=−1 sector as shown in Figure 5. This effective Hamiltonian forms a 8×8 matrix in the occupation number basis states |φ k of theP ν 0 +P 0 2 sector, which yields the following matrix structure where F = m J |m| (E) 2 2(U + J H + mω) , | ↓, ↑ i | ↓, ↑ j , | ↑, ↓ i | ↓, ↑ j , | ↓, ↑ i | ↑, ↓ j , where |σ a , σ b i =ĉ † ibσ ĉ † iaσ |0 . The lower right block of the matrix corresponds to the effective Hamiltonian inP 0 2 sector withm=−1 where the basis is the following We diagonalize the matrix, Eq. (105) and, after time averaging, obtain the spectrum shown in Figure 5. We emphasize that the diagonalization and time averaging do not commute e.g.

T
T 0 I(t)dt=0 which washes out the coupling.