## Floquet engineering of Kitaev quantum magnets

### Multi-orbital system

We focus on the edge-shared iridates (5d) and ruthenates (4d) as our target materials. In the octahedral crystal field, these materials have been proposed to be described by the *J*–*K*–Γ Hamiltonian^{26} and as possible candidates to realize Kitaev QSL phase. Here, the *t*_{2g} manifold of the transition metal (TM) is known to host a hole in the *J*_{eff} = 1/2 sector, which is separated from *J*_{eff} = 3/2 states due to a large spin-orbit coupling (*λ*) in the system. The *J*_{eff} = 1/2 multiplet is composed of *d*_{xy}, *d*_{yz}, and *d*_{zx} orbitals of the TM atom. On the other hand, the TM atoms form a two-dimensional (2D) honeycomb lattice perpendicular to [111] direction composed of the standard *x*, *y*, *z-*bonds as shown in Fig. 1a. For subsequent discussion, we constrain our analysis to the *z*-bond which is oriented along the [110] direction. The *x* and *y* bonds can be recovered by *C*_{3}-rotation of this bond. Along the *z*-bond, only *p*_{z}-orbital of the oxygen/chloride ligand is active and has a finite hopping between the *d*_{yz} and *d*_{zx} orbitals as shown in panels (b) and (c) of Fig. 1, respectively. These multi-orbital systems are well captured by Hubbard-Kanamori Hamiltonian and can be mapped to effective spin Hamiltonians using exact diagonalization (ED)^{7} or perturbation theory^{26,30}. The Hamiltonian of the system is given by \({{{{{{{\mathcal{H}}}}}}}}(t)={{{{{{{{\mathcal{H}}}}}}}}}_{K}(t)+{\sum }_{i = 1,2}{{{{{{{{\mathcal{H}}}}}}}}}_{i}^{C}\) where

$${{{{{{{{\mathcal{H}}}}}}}}}_{i}^{C}=\; \frac{U}{2}\mathop{\sum}\limits_{\alpha ,\sigma }{n}_{i\alpha \sigma }{n}_{i\alpha \bar{\sigma }}+\frac{U^{\prime} }{2}\mathop{\sum}\limits_{\alpha \ne \beta \atop\sigma ,\sigma ^{\prime}}{n}_{i,\alpha ,\sigma }{n}_{i\beta \sigma ^{\prime} }\\ -\frac{{J}_{{{{{{{{\rm{H}}}}}}}}}}{2}\mathop{\sum}\limits_{\alpha \ne \beta \atop\sigma ,\sigma ^{\prime}}{d}_{i\alpha \sigma }^{{{{\dagger}}} }{d}_{i\alpha \sigma ^{\prime} }{d}_{i\beta ,\sigma ^{\prime} }^{{{{\dagger}}} }{d}_{i\beta \sigma }\\ +\frac{{J}_{{{{{{{{\rm{P}}}}}}}}}}{2}\mathop{\sum}\limits_{\alpha \ne \beta \atop\sigma}{d}_{i\alpha \sigma }^{{{{\dagger}}} }{d}_{i\alpha \bar{\sigma }}^{{{{\dagger}}} }{d}_{i\beta \sigma }{d}_{i\beta \bar{\sigma }}+{{\Delta }}\mathop{\sum }\limits_{\sigma }{p}_{i\sigma }^{{{{\dagger}}} }{p}_{i\sigma }.$$

(1)

Here \({d}_{i\alpha \sigma }^{{{{\dagger}}} }\) (\({p}_{i\sigma }^{{{{\dagger}}} }\)) creates a hole on TM (ligand) in the *i*th-cell, which includes a TM and a ligand site with spin *σ* and orbital *α* = {*x**y*, *y**z*, *z**x*}, \({n}_{i\alpha\sigma }={d}_{i\alpha \sigma }^{{{{\dagger}}} }{d}_{i\alpha \sigma }^{{{\phantom{{\dagger}}}} }\), *U* (\(U^{\prime}\)) is the strength of the intra (inter)-orbital Coulomb repulsion, *J*_{H} is the Hund’s coupling for the orbitals *α*, *β* ∈ {*d*_{xy}, *d*_{yz}, *d*_{zx}}, and *J*_{P}(=*J*_{H}) is the pair hoppings. Δ parametrizes the ligand charge-transfer energy. Note that we have ignored the atomic spin-orbit coupling (SOC) term, \(\lambda /2{\sum }_{i}{\sum }_{{{\alpha \beta ,}\atop {\sigma \sigma ^{\prime}}}}{d}_{i\alpha \sigma }^{{{{\dagger}}} }({{{{{{{{\bf{L}}}}}}}}}_{\alpha \beta }\cdot {{{{{{{{\bf{S}}}}}}}}}_{\sigma \sigma ^{\prime} }){d}_{i\beta \sigma ^{\prime} }\), assuming \(\lambda \ll U,U^{\prime} ,{J}_{{{{{{{{\rm{H}}}}}}}}},{{\Delta }}\) ^{26,31}.

Focusing on the four-site cluster [two TM atoms oriented along the *z*-bond with their two edge-shared ligand sites], we now write down the tight-binding Hamiltonian in presence of a circularly polarized light (CPL) as

$${{{{{{{{\mathcal{H}}}}}}}}}_{K}(t)= \, \mathop{\sum}\limits_{\alpha ,\beta ,\sigma}{e}^{i\delta F(t)}{t}_{\alpha \beta }{d}_{1\alpha \sigma }^{{{{\dagger}}}}{d}_{2\beta \sigma} \\ +{t}_{pd}\mathop{\sum}\limits_{\sigma}\bigg({e}^{i\delta {F}_{x}(t)}{d}_{1zx\sigma }^{{{{\dagger}}}}{p}_{1\sigma}-{e}^{i{F}_{y}(t)}{p}_{1\sigma}^{{{{\dagger}}}}{d}_{2yz\sigma}\\ +\,{e}^{i\delta {F}_{y}(t)}{d}_{1yz}^{{{{\dagger}}}} {p}_{2\sigma}-{e}^{i{F}_{x}(t)}{p}_{2\sigma}^{{{{\dagger}}}}{d}_{2zx\sigma}\bigg)+{{{{{{{\rm{h.c.}}}}}}}},$$

(2)

where *t*_{pd} is the hopping amplitude between the *p*– and *d*-orbitals. The hopping matrix between the three *d*-orbitals *t*_{αβ} is obtained through Slater–Koster^{32} interatomic matrix elements as

$${t}_{\alpha \beta }{d}_{1\alpha \sigma }^{{{{\dagger}}}}{d}_{2\beta \sigma}=\left[\begin{array}{ccc}{d}_{1yz\sigma }^{{{{\dagger}}} }&{d}_{1zx\sigma }^{{{{\dagger}}} }&{d}_{1xy\sigma }^{{{{\dagger}}} }\end{array}\right]\left[\begin{array}{ccc}{t}_{3}&{t}_{2}&0\\ {t}_{2}&{t}_{3}&0\\ 0&0&{t}_{1}\\ \end{array}\right]\left[\begin{array}{c}{d}_{2yz\sigma }\\ {d}_{2zx\sigma }\\ {d}_{2xy\sigma }\end{array}\right].$$

(3)

The time-dependence of the drive appears as a Peierls phase in these multi-orbital Hamiltonian as noted in Eq. (2)^{7,33}. For the sake of simplicity, we have considered the light along [001] direction with circular polarization, \({{{{{{{\bf{F}}}}}}}}(t)={E}_{0}(\hat{{{{{{{{\bf{x}}}}}}}}}\sin \omega t+\hat{{{{{{{{\bf{y}}}}}}}}}\cos \omega t)/\omega\). The phase along the TM–TM bond \(\left(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\right)\) is given by \(\delta F(t)=\zeta \sin (\omega t-\pi /4)\) and the phases along TM–ligand bonds are given by \(\delta {F}_{x}(t)=\zeta \sin \omega t/\sqrt{2}\) and \(\delta {F}_{y}(t)=\zeta \cos \omega t/\sqrt{2}\). Here, the dimensionless parameter *ζ* is defined as *ζ* = *R*_{ij}*E*_{0}/*ω* [*R*_{ij} is the bond-length between the sites for the associated hopping]. The values of the parameters entering Eqs. (1)–(3) are adapted from the recent ab-initio^{34} and photoemission^{35} studies relevant for *α*-RuCl_{3}: *t*_{1} = 0.044, *t*_{2} = 0.08, *t*_{3} = 0.109, *t*_{pd} = −0.8, *U* = 3.0, *J*_{H} = 0.45, and Δ = 5.

Here, we assumed a reduced value of *t*_{2} so as to account for finite *t*_{pd}. All the values are in units of eV, unless stated otherwise. We also use relations, \(U^{\prime} =U-2{J}_{{{{{{{{\rm{H}}}}}}}}}\) and *J*_{P} = *J*_{H} in the Hubbard–Kanamori Hamiltonian, which preserves the rotational invariant form. Consequently, the Hubbard–Kanamori term in Eq. (1) can be mapped to

$${{{{{{{{\mathcal{H}}}}}}}}}_{i}^{C}=\frac{U-3{J}_{{{{{{{{\rm{H}}}}}}}}}}{2}{\left({N}_{i}-5\right)}^{2}-2{J}_{{{{{{{{\rm{H}}}}}}}}}{{{{{{{{\bf{S}}}}}}}}}_{i}^{2}-\frac{{J}_{{{{{{{{\rm{H}}}}}}}}}}{2}{{{{{{{{\bf{L}}}}}}}}}_{i}^{2}+{{\Delta }}\mathop{\sum}\limits_{\sigma }{n}_{i\sigma }^{p},$$

(4)

where *N*_{i} is electron number, **S**_{i} is total spin, and **L**_{i} is total orbital angular momentum at site *i*^{26,31}. The above simplified form of the Hubbard-Kanamori Hamiltonian [Eq. (4)] allows us to evaluate the energy of the doubly occupied states at the TM sites. Since, the SOC, *λ* is ignored, we label the doubly occupied sites in terms of the orbital and spin angular momentum **L**, **S**, respectively, as \({\left|\right.}^{2S+1}L,{L}_{z},{S}_{z},0\left.\right\rangle\), where *L*, *S *[=(0, 0), (1, 1), (2, 0)] are the orbital and the spin angular momentum of the doublet, respectively, with *L*_{z} and *S*_{z} being their respective *z*-components. The allowed values of *L* and *S* provides three different energies for the doublets (two particles on TM site) given by: *E*_{S} = *U* + 2*J*_{H}, *E*_{P} = *U* − 3*J*_{H}, *E*_{D} = *U* − *J*_{H} and singly occupied case (one particle in each TM site), *E*_{0} = 0.

### Effective spin-exchange model

In this section, we present the effective spin-exchange model derived from the multi-orbital model in presence of a drive discussed above. We evaluate the various parameters in the spin-exchange model for two cases; without ligand and with ligands in the multi-orbital model. In each case, the results are evaluated numerically using ED and analytically using pertubation theory. We further present a comparision between the numerical and analytical results.

Since the tight-binding parameters are small compared to interaction strengths [\({t}_{\alpha \beta },{t}_{pd}\ll U,U^{\prime} ,{J}_{{{{{{{{\rm{H}}}}}}}}},{{\Delta }}\)], the model can be analyzed by a low-energy effective spin-exchange Hamiltonian by ignoring the high-energy charge degrees of freedom. The static limit of our model has been extensively studied previously leading to the well-known *J*−*K*−Γ spin Hamiltonian^{26,36} (also see Supplementary Note 1). Here, we utilize the Floquet approach to derive a time-independent effective Hamiltonian from Eqs. (1, 2) and subsequently perform ED calculations to estimate the various spin-exchange parameters of the concomitant spin Hamiltonian. On the other hand, the time-independent version of the latter can be naturally captured by a minimal magnetic model [based on the mirror reflection, inversion, and broken time-reversal symmetry of TM–ligand–TM atomic cluster associated with the *z*-bond] as

$${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{eff}}}}}}}}}=\left[{S}_{1}^{x}\,{S}_{1}^{y}\,{S}_{1}^{z}\right]\,\left[\begin{array}{lll}J&{{\Gamma }}&{{\Gamma }}^{\prime} \\ {{\Gamma }}&J&{{\Gamma }}^{\prime} \\ {{\Gamma }}^{\prime} &{{\Gamma }}^{\prime} &J+K\end{array}\right]\,\left[\begin{array}{l}{S}_{2}^{x}\\ {S}_{2}^{y}\\ {S}_{2}^{z}\end{array}\right]+{{{{{{{\bf{h}}}}}}}}\cdot ({{{{{{{{\bf{S}}}}}}}}}_{1}+{{{{{{{{\bf{S}}}}}}}}}_{2}).$$

(5)

Since the trigonal distortion is ignored in our model, we have \({{\Gamma }}^{\prime} =0\)^{26}. The emergent Zeeman magnetic field **h** is an effect of the broken time-reversal symmetry due to the applied CPL and is not accounted for in the static version of the associated model.

We now briefly discuss how the mapping of the multi-orbital model to the spin model is carried out. In the method section, we discuss how the time-dependent Hamiltonian can be mapped to a time-independent form. The time-independent matrix in Bloch structure form is given by, \(({\epsilon }_{\alpha }+m\omega )|{\psi }_{\alpha ,m}\rangle ={\sum }_{m^{\prime} }{H}_{m-m^{\prime} }|{\psi }_{\alpha ,m^{\prime} }\rangle\). Here, *H*_{m} and \(|{\psi }_{\alpha ,m}\rangle\) are Hamiltonian and wavefunction in the frequency basis for *m*th Floquet sector index [see “Methods” section for details]. Here index *α* consists of (i) a single hole in each TM site and an empty ligand site, (ii) two holes [*double occupancy*] on a TM site while the ligand and the other TM site is empty and (iii) one particle each on a ligand and a TM site (see Supplementary Note 2 for the exact forms). The singly occupied states are a product of *J*_{eff} = 1/2 states on the two TM sites, the doublets consists of states \({\left|\right.}^{2S+1}L,{L}_{z},{S}_{z},0\left.\right\rangle\), where *L*, *S*[=(0, 0), (1, 1), (2, 0)], as mentioned earlier.

The full Floquet Hamiltonian (*H*) can be written in the eigen-decomposed form as, \(H={\sum }_{{\mathfrak{n}}}{E}_{{\mathfrak{n}}}\|{\phi }_{{\mathfrak{n}}}\rangle \langle {\phi }_{{\mathfrak{n}}}|\), where, \({E}_{{\mathfrak{n}}}\) and \(\left|{\phi }_{{\mathfrak{n}}}\right\rangle \,(={\sum }_{\alpha ,m}{a}_{\alpha ,m}^{{\mathfrak{n}}}|{\psi }_{\alpha ,m}\rangle )\) are the eigenvalues and eigenvectors, respectively. Notice that \(|{\psi }_{\alpha ,m}\rangle\) contains all the configurations for the two-site problem. In order to restrict the basis to *m* = 0-Floquet sector and singly occupied basis, one can use a projector; *P*_{s,0} = *P*_{α∈s}*P*_{m=0}. We then have the projected spin Hamiltonian given by *H*_{PG} = *P*_{s,0}*H**P*_{s,0}. The *H*_{PG} can be used to estimate the various parameters in \({{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{eff}}}}}}}}}\). For examples, \(J=2\left\langle -\frac{1}{2},+\frac{1}{2}\right|{H}_{PG}\left|+\frac{1}{2},-\frac{1}{2}\right\rangle\) (see Supplementary Note 3 for additional details). The spin model is valid only when the Floquet band for the singly occupied is well separated from the other Floquet band and the upper Hubbard band. We define \({{{\Delta }}}_{l}=({E}_{\min }(m=0,s)-{E}_{{{{{{{{\rm{min-1}}}}}}}}})/W\) and \({{{\Delta }}}_{u}=({E}_{{{{{{{{\rm{max+1}}}}}}}}}(m=0,s)-{E}_{\max }(m=0,s))/W\) for each *ζ*, i.e., Δ_{u}, Δ_{l} ≪ 1. Here \({E}_{\min }(m=0,s)\) and \({E}_{\max }(m=0,s)\) are the minimum and maximum energy for the singly occupied states in the *m* = 0 Floquet sector, and \(W={E}_{\max }(m=0,s)-{E}_{\min }(m=0,s)\) for a given *ζ*.

Figures 2 and 3 plot the estimates using ED for the system without and with ligand, respectively. Inset in the panels of these figures plots the energy states for the *m* = 0 Floquet states with singly occupied configuration (in black) and its nearby states. The singly occupied states (in black) are separated from other states, which validates our calculation for the drive frequencies presented in our work.

#### Spin-exchange model: ED without ligands

The ligand degrees of freedom (oxide/chloride *p* orbitals), in the edge-shared ruthenates or iridates, are usually integrated out leading to an effective description of the model in terms of the TM *d*-orbitals only^{26}. In this section, we perform ED on the Floquet Hamiltonian obtained from \({{{{{{{\mathcal{H}}}}}}}}(t)\) and analyze the results by integrating out the ligand. The system without the ligand is simulated by turning off the direct hopping *t*_{pd} between TM and ligand sites and rescaling the TM–TM hopping *t*_{2} to \({t}_{2}+{t}_{pd}^{2}/{{\Delta }}\) in Eq. (2). Consequently, the four-site cluster effectively becomes a two-site system. In this case, the emergent Zeeman magnetic field **h** term is absent as there is no residual orbital current along the TM–TM bond. However, the multi-orbital nature leads to finite non-vanishing spin-exchange couplings *J*, *K*, and Γ.

The corresponding results are shown in Fig. 2. We choose four distinct light frequencies avoiding the three critical resonant frequencies *ω* = {*E*_{P}, *E*_{D}, *E*_{S}}[={1.65, 2.45, 3.9}] of the doublet states and illustrate the relative variation of *J*, *K*, and Γ with applied laser strength in the four panels [Fig. 2]. Inset in each panel shows the separation of the *m* = 0 Floquet singly occupied states from the adjacent *m* ≠ 0 states. Panel (a) shows the variation of the couplings at drive frequency *ω* = 1.2, below the doublet *E*_{P} energy. In this case, we notice that *J*, *K*, Γ can all change their sign and enhance the magnitude, but the relative tunability is missing. Panel (b) shows the result for *ω* = 2.1, in between doublet *E*_{P} and *E*_{D} energies and retains the similar feature as in the prior case. Panel (c) and (d) show the results for *ω* = 3.2 [*E*_{D} < *ω* < *E*_{S}] and *ω* = 7 [*ω* > *E*_{S}], respectively. In these two cases, however, we notice that the sign of the parameters cannot be changed, and also the magnitudes of the couplings can be enhanced only mildly. As will be shown by perturbation calculations below, the ratio of Γ/*K* is fixed, independent of the amplitude and frequency of the light.

#### Spin-exchange model: ED with ligands

We now come to the main part of our paper. Here, we present the results using the full Hamiltonian given by Eqs. (1, 2) for the parameters discussed previously and consider the effect of ligands in full glory. Recently, the inclusion of ligands has been investigated in a few studies^{37}. The inclusion of ligand degree of freedom in the multi-orbital description is crucial primarily for its two-fold significance. First, in real materials [iridates/ruthenates] the ligand effect is unavoidably present and a thorough quantitative material-specific analysis needs to incorporate it properly. Second and most importantly, the inclusion of ligand *p*-orbital along with the TM *d*-orbital offers the desired relative tunability of the spin-exchange couplings—*J*, *K*, and Γ [see Fig. 3], which was absent as we saw in the previous section. Finally, the inclusion of ligand sites constitutes the full TM–ligand–TM atomic cluster as shown in Fig. 4, with its associated orbital current which leads to a finite Zeeman magnetic field through inverse Faraday effect^{28,29}. In Fig. 3 we show the variation of all the parameters defined in Eq. (5) with the driving strength (*ζ*) evaluated using ED calculations.

As before, we again choose six distinct light frequencies avoiding the four critical resonant frequencies *ω* = {*E*_{P}, *E*_{D}, *E*_{S}, Δ} [={1.65, 2.45, 3.9, 5}] of the doublet states and illustrate the relative variation of *J*, *K,* and Γ with applied laser strength in the six panels (a)–(f) [Fig. 3]. The ED analysis is always performed away from the Floquet resonances in order to avoid the breakdown of the spin picture. Inset in each panel shows the separation of the *m* = 0 Floquet singly occupied states from the adjacent *m* ≠ 0 Floquet states. For field strength *ζ* = 0 [*static case*], we have *J* = 0.0007, *K* = 0.005, and Γ = 0.003 in each panel, whereas the Zeeman magnetic field **h** is zero. The vanishing of the magnetic field is due to the preserved time-reversal symmetry of the multi-orbital system in the absence of CPL.

Panel (a) shows the variation of the various parameters defined in Eq. (5) with the driving amplitude *ζ* at the drive frequency *ω* = 1.2 eV, below the lowest doublet energy *E*_{P}. Indeed for the reasons mentioned earlier, we find a finite Zeeman magnetic field **h** in the driven system, in contrast to vanishing **h** in the case without ligand. In this case, the parameters *J*, *K*, Γ, and **h** can all change their sign and their magnitudes can be enhanced. In addition, these parameters vanish at different values of *ζ* [compared to the situation in Fig. 2] allowing for their relative tuning. We also note that the magnitude of the Kitaev term *K* and anisotropy Γ can be enhanced by a few orders of magnitude compared to the Heisenberg exchange term *J*. Panel (b) shows the variation of the parameters with *ζ* for *ω* = 2.3 in between the *E*_{P} and *E*_{D} doublet. In addition to our findings in panel (a), we notice that *K* at *ζ* = 3.7 can be tuned to a few orders of magnitudes larger compared to all the other parameters, making it ideal for realizing the Kitaev QSL phase. The results for *ω* = 3.2 [*E*_{D} < *ω* < *E*_{S}] are shown in panel (c). In this case, the magnetic field **h** can achieve the largest magnitude compared to the other parameters. Moreover, we notice limited tunability at larger drive strength *ζ*. Panel (d) shows the result for *ω* = 4.5, in the regime between *E*_{S} and Δ. In this case, the system can be tuned to large *J* and *K*, leading to a realization of the Kitaev-Heisenberg model. In addition, the almost vanishing Zeeman field provides a promising route to realize gapless Kitaev QSL. For *ω* = 7 above the largest energy scale Δ, in our model [see panel (e)], we notice a significant enhancement of *K* at low drive strength, making it another frequency regime suited for realizing the Kitaev phase. Panel (f) shows the result for *ω* = 12, large frequency, well above Δ. In this case, we observe that all the parameter values decrease with the driving strength similar to the case without ligand [see Fig. 2d]. In contrast, the tunability of the signs of these parameters persists.

### Spin-exchange model: perturbation theory

So far we analyzed our model in the Floquet approximation utilizing the ED calculations and discussed the results in various regimes of the model-based parameter values. In this section, we derive the low-energy spin-exchange Hamiltonian based on perturbation theory and compare the analytical expressions of the couplings *J*, *K*, Γ, and **h** with the numerical estimates from the ED.

In the absence of ligand degree of freedom, we perform a second-order perturbation theory (see Supplementary Note 2) and time-dependent Schrieffer–Wolff transformation^{28} to obtain the various couplings defined in Eq. (5) as

$$J=\frac{4}{27}\mathop{\sum }\limits_{m=-\infty }^{\infty }{{{{{{{{\mathcal{J}}}}}}}}}_{m}^{2}(\zeta )\left[\frac{6{t}_{1}({t}_{1}+2{t}_{3})}{{E}_{P}+m\omega }+\frac{2{({t}_{1}-{t}_{3})}^{2}}{{E}_{D}+m\omega }+\frac{{(2{t}_{1}+{t}_{3})}^{2}}{{E}_{S}+m\omega }\right],$$

(6a)

$$K=\frac{8}{9}\mathop{\sum }\limits_{m=-\infty }^{\infty }{{{{{{{{\mathcal{J}}}}}}}}}_{m}^{2}(\zeta )\frac{{J}_{{{{{{{{\rm{H}}}}}}}}}[{({t}_{1}-{t}_{3})}^{2}-3{t}_{2}^{2}]}{({E}_{P}+m\omega )({E}_{D}+m\omega )},$$

(6b)

$${{\Gamma }}=\frac{16}{9}\mathop{\sum }\limits_{m=-\infty }^{\infty }{{{{{{{{\mathcal{J}}}}}}}}}_{m}^{2}(\zeta )\frac{{J}_{{{{{{{{\rm{H}}}}}}}}}{t}_{2}({t}_{1}-{t}_{3})}{({E}_{P}+m\omega )({E}_{D}+m\omega )},$$

(6c)

where \({{{{{{{{\mathcal{J}}}}}}}}}_{m}(\zeta )\) is the Bessel function of the first kind of order *m*. Notice that the expressions for *K* and Γ have nodes at the identical values of the drive strength *ζ* dictated by the combination of the Bessel functions and the denominator only and also are found to follow a constant relative strength as illustrated in Fig. 2. Therefore, it forbids the relative tunability of these parameters in the system. The comparison of our ED results with the analytical expressions for the couplings *J*, *K*, and Γ is shown in Fig. 5a. We find a very good agreement between the two in this case.

We now focus on the final part of our work and obtain the spin-exchange model using perturbation theory and time-dependent Schrieffer–Wolff transformation^{28}. In this case, the ligand degrees of freedom are taken into consideration and we have to rely on a third-order perturbation theory to properly capture the ligand effects. The perturbation theory is performed utilizing the processes illustrated in Fig. 4 and we obtain the couplings defined in Eq. (5) as

$$J=\frac{8{t}_{pd}^{2}}{27}\mathop{\sum }\limits_{n,m=-\infty }^{\infty }{{\mathfrak{J}}}_{m,n}(\zeta )\left[\frac{\sin [(m-n){\psi }_{0}]}{{{\Delta }}+m\omega }\left(\frac{2{t}_{1}+{t}_{3}}{{E}_{S}+l\omega }+\frac{{t}_{1}-{t}_{3}}{{E}_{D}+l\omega }+3\frac{{t}_{1}+{t}_{3}}{{E}_{P}+l\omega }\right)\right],$$

(7a)

$$K=\frac{8{t}_{pd}^{2}}{9}\mathop{\sum }\limits_{n,m=-\infty }^{\infty }{{\mathfrak{J}}}_{m,n}(\zeta )\left[\frac{{J}_{{{{{{{{\rm{H}}}}}}}}}}{({{\Delta }}+m\omega )}\frac{\sin [(m-n){\psi }_{0}]({t}_{1}-{t}_{3})-3\cos [(m-n){\psi }_{0}]{t}_{2}}{({E}_{P}+l\omega )({E}_{D}+l\omega )}\right],$$

(7b)

$${{\Gamma }}=\frac{8{t}_{pd}^{2}}{9}\mathop{\sum }\limits_{n,m=-\infty }^{\infty }{{\mathfrak{J}}}_{m,n}(\zeta )\left[\frac{{J}_{{{{{{{{\rm{H}}}}}}}}}}{({{\Delta }}+m\omega )}\frac{\cos [(m-n){\psi }_{0}]({t}_{1}-{t}_{3})+\sin [(m-n){\psi }_{0}]{t}_{2}}{({E}_{P}+l\omega )({E}_{D}+l\omega )}\right],$$

(7c)

$$h=\frac{4{t}_{pd}^{2}}{9}\mathop{\sum }\limits_{n,m=-\infty }^{\infty }{{\mathfrak{J}}}_{m,n}(\zeta )\frac{\sin [(m-n){\psi }_{0}]}{{{\Delta }}+m\omega }\left[\frac{{t}_{1}-{t}_{3}}{{E}_{P}+l\omega }+\frac{{t}_{1}-{t}_{3}}{{E}_{D}+l\omega }\right],$$

(7d)

where \({{\mathfrak{J}}}_{m,n}(\zeta )={{{{{{{{\mathcal{J}}}}}}}}}_{m+n}(\zeta ){{{{{{{{\mathcal{J}}}}}}}}}_{-m}(\zeta {r}_{ij}/{R}_{ij}){{{{{{{{\mathcal{J}}}}}}}}}_{-n}(\zeta {r}_{ij}/{R}_{ij})\), *r*_{ij} is the bond-length between the TM and ligand sites, *l* = *m* + *n* and *ψ*_{0} is the angle between the TM–TM and TM–ligand bond. Note that previously *ψ*_{0} was considered to be equal to *π*/4. We observe that in the static case *ζ* = 0 the super-exchange coupling *J* and the Zeeman magnetic field **h** vanish. This vanishing of super-exchange coupling is consistent with the Jackeli–Khaliluin formalism^{23} and the absence of **h** at *ζ* = 0 is consistent with the intact time-reversal symmetry of the multi-orbital system.

We further observe the distinct analytical structure of the Kitaev term *K* and the anisotropy Γ [7b, c)] which dictates that the zeros of these parameters occur at different drive strength *ζ*. This is in stark contrast to the results discussed without ligand in Eq. (6b, c) and is consistent with the variation of these parameters as shown in Fig. 3. Figure 5b shows the results for the relative comparison between the perturbation [adding the contributions from Eq. (6a–c) to Eq. (7a–d)] and ED calculations. In this case, we find an excellent agreement for the case of the magnetic field **h**, but have a relatively poor agreement in the case of the other parameters. Note that with the inclusion of ligand, there are many more paths possible for the super-exchange in the high order processes. Therefore, this deviation can possibly be accounted for by higher-order contributions, which cannot be neglected due to larger values of *t*_{pd}, whereas our perturbation term accounts only till the third order. On the other hand, the better agreement of the perturbation results and the ED calculations for the Zeeman magnetic field **h** is justified. In this case, there is only one process, in the next order perturbation, which contributes and is roughly proportional to \({t}_{pd}^{4}/{{{\Delta }}}^{2}U\). Within our parameter choice, this term is much smaller than the second or third-order contributions and can safely be neglected.