Floquet engineering of Kitaev quantum magnets

In recent years, there has been an intense search for materials realizing the Kitaev quantum spin liquid model. A number of edge-shared compounds with strong spin-orbit coupling, such as RuCl3 and iridates, have been proposed to realize this model. Nevertheless, an effective spin Hamiltonian derived from the microscopic model relevant to these compounds generally contains terms that are antagonistic toward the quantum spin liquid. This is consistent with the fact that the zero magnetic field ground state of these materials is generally magnetically ordered. It is a pressing issue to identify protocols to drive the system to the limit of the Kitaev quantum spin model. In this work, we propose Floquet engineering of these Kitaev quantum magnets by coupling materials to a circularly polarized laser. We demonstrate that all the magnetic interactions can be tuned in situ by the amplitude and frequency of the laser, hence providing a route to stabilize the Kitaev quantum spin liquid phase. Kitaev materials are prime candidates to study quantum spin liquids but despite many known examples it is still a challenge to realize systems that have a total absence of magnetic ordering. Here, the authors propose Floquet engineering and the use of light-matter interactions to achieve a Kitaev quantum spin liquid phase.

L ight-matter interaction not only provides an important way to probe materials' properties but also offers exciting opportunities to control the physical properties of quantum materials [1][2][3] . It has been demonstrated experimentally that light can induce superconductivity 4-6 , magnetism 7,8 , topological states [9][10][11] , and other novel quantum states of matter 12 . One particular fruitful protocol is the so-called Floquet engineering by driving materials periodically with a continuous laser. In this Floquet approach, the coupling between light and matter can tune the microscopic parameters, generate new interactions, and even stabilize entirely new states and excitations that do not exist in equilibrium 13,14 . Inspired by the success of light control of quantum states in the past, we investigate the light-driven transition into the quantum spin liquid in strongly correlated magnetic materials.
Quantum spin liquid (QSL) is one class of highly entangled quantum states of matter, which hosts fractionalized excitations 15,16 . The elegant exact solution of the Kitaev spin model unambiguously shows the existence of a quantum spin liquid 17 . The presence of an external magnetic field can open a gap and the model supports the Majorana fermions which have the potential for applications in robust topological quantum computation 18,19 . Kitaev's work has motivated tremendous efforts both theoretically and experimentally to materialize the Kitaev model in quantum magnets [20][21][22] . Jackeli and Khaliullin showed that one can realize the Kitaev spin interaction in edgeshared octahedral structure with strong spin-orbit coupling 23 . Several candidate materials, including RuCl 3 24 and iridates, have been identified and have been extensively studied. Encouraging signs of QSL has been reported 25 , but unambiguous identification of the QSL remains a challenge. One obstacle from a theoretical perspective is that there exist competing magnetic interactions originated from the different correlation induced exchange processes in materials, which drives the model Hamiltonian away from the ideal Kitaev QSL model 26,27 . Numerical study shows that the Kitaev QSL is stable only in a small region of the model parameter space 26 . Therefore, it is highly desirable to design a protocol for continuously in situ tuning the magnetic interactions in materials to drive the system into QSL phases.
In this work, we propose to use light to tune magnetic interactions in spin-orbit coupled Mott insulators to favor QSL. We consider the multi-orbital strongly spin-orbit coupled Hubbard model, which describes several candidate materials including RuCl 3 and iridates. In a simple picture, the coupling of electrons to light has two effects. First, the periodic modulation of the hopping of electrons between different sites effectively normalizes the hopping parameters. Second, electrons can absorb or emit photons during the virtual hopping process, and as a consequence, the energy barrier due to Coulomb repulsion also gets modified. Therefore, the strength and even the sign of magnetic interactions can be controlled by light. In addition, an effective Zeeman field can be generated by circularly polarized light through the process called inverse Faraday effect, which provides a new handle to control the quantum states 28,29 . Here, we derive a low-energy spin Hamiltonian from the Hubbard model, both using the numerical exact diagonalization (ED) method and the analytical perturbation theory, and explore in detail how the magnetic interactions can be tuned by the frequency and amplitude of the light.

Results
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 zbond, 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 Here d y iασ (p y iσ ) creates a hole on TM (ligand) in the ith-cell, which includes a TM and a ligand site with spin σ and orbital α = {xy, yz, zx}, n iασ ¼ d y iασ d iασ , U (U 0 ) 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, λ=2∑ i ∑αβ; σσ 0 d y iασ ðL αβ Á S σσ 0 Þd iβσ 0 , assuming λ ( U; U 0 ; J H ; Δ 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 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 The time-dependence of the drive appears as a Peierls phase in these multi-orbital Hamiltonian as noted in Eq. 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 0 ¼ U À 2J 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 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 j 2Sþ1 L; L z ; S z ; 0i, 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: 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 αβ ; t pd ( U; U 0 ; J H ; Δ], 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 Since the trigonal distortion is ignored in our model, we have Γ 0 ¼ 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, Here, H m and jψ α;m i are Hamiltonian and wavefunction in the frequency basis for mth 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 j 2Sþ1 L; L z ; S z ; 0i, where L, S[=(0, 0), (1, 1), (2, 0)], as mentioned earlier.
The full Floquet Hamiltonian (H) can be written in the eigendecomposed form as, H ¼ ∑ n E n kϕ n ihϕ n j, where, E n and ϕ n ð¼∑ α;m a n α;m jψ α;m iÞ are the eigenvalues and eigenvectors, respectively. Notice that jψ α;m i 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 HP s,0 . The H PG can be used to estimate the various parameters in (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 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 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 2 pd =Δ 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 multiorbital 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 spinexchange 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. 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 modelbased 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 Fig. 4 Schematics for evaluating spin Hamiltonian parameters at the third-order in perturbation theory. We start with one particle each at the transition metal (TM) site. In the intermediate states, the particle moves to a ligand and then to the other TM site creating a doublet. Finally, it returns to the original site, with or without a spin-flip. We also consider a timereversal partner of this hopping path.
the various couplings defined in Eq. (5) as where J m ðζÞ 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 timedependent 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 h ¼ 4t 2 pd 9 ∑ 1 n;m¼À1 J m;n ðζÞ sin½ðm À nÞψ 0 Δ þ mω where J m;n ðζÞ ¼ J mþn ðζÞJ Àm ðζr ij =R ij ÞJ Àn ðζ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 superexchange 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 4 pd =Δ 2 U. Within our parameter choice, this term is much smaller than the second or third-order contributions and can safely be neglected.

Discussion
In our discussion above, we have not taken into account the role of heating due to laser irritation. When the laser frequency is close to resonances in the system, electrons can be excited efficiently resulting in severe heating. Under the drive, the model Hamiltonian can be divided into different Floquet Hamiltonian sectors according to the number of the photon, m, the system absorbs/emits. Initially, electrons occupy the m = 0 sector. The laser drives the electrons to other Floquet sectors, which causes heating. Meanwhile, the effective Floquet Hamiltonian at m = 0 sector is modified through hybridization with other sectors. There are resonances within the same Floquet sector. In our system, at the single-ion level, we have crystal field splitting between t 2g and e g orbitals, and splitting between J eff = 3/2 and J eff = 1/2 manifold due to the spinorbit coupling. The frequency of the laser should be tuned away from these resonances. In the Mott insulator, one dominant mechanism of heating is the excitation of the doublons by laser. This resonance occurs when the photon energy is close to the energy difference between different Hubbard bands. The population of the doublons increases continuously with time in the presence of laser irritation, which eventually invalidates the effective spin Hamiltonian description. It is calculated in ref. 38 that the rate of doublon generation can be very low when the laser frequency is away from resonance set by the lower and upper Hubbard band. The effective spin Hamiltonian we derived here is valid in a long time window, where the effect of doublons on magnetic interactions is also negligible. This long time window is known as the Floquet prethermal region, which can be exponentially long in time before the system enters into a featureless infinite temperature region due to heating [38][39][40][41][42][43][44][45] . In addition, the heating can be mitigated by coupling the driven system to dissipative bath, such as phonon bath 46,47 .
We remark that not all the magnetic interactions can be tuned independently because we only have two tuning parameters, i.e., amplitude and frequency, of the incident circularly polarized light. As demonstrated explicitly in Fig. 3, we can achieve the region where K is dominant over all other magnetic interactions, thus realizing a favorite condition for the Kitaev QSL. Furthermore, the Kitaev interaction K can be tuned to be ferromagnetic or antiferromagnetic by light in a single material. It is believed that RuCl 3 and iridates are described by a ferromagnetic Kitaev interaction without a light drive. The antiferromagnetic Kitaev model has attracted considerable interest recently because it may host a new gapless QSL in the intermediate magnetic field region [48][49][50][51][52][53] . Both the antiferromagnetic Kitaev interaction and magnetic field can be induced by light, thus the driven system allows us to access this gapless QSL. We can also tune into a parameter region where Γ interaction is dominant, where a new type of QSL has been suggested recently 54 .
So far, we have focused on the circularly polarized light, our results can be readily generalized to the case with a linearly polarized light. In this case, the effective Zeeman field due to the inverse Faraday effect is absent because the time-reversal symmetry is preserved. It is natural to expect one can tune the system into an anisotropic spin Hamiltonian by linearly polarized light, which can also support QSL 33 .
While we are working on the current project, we noticed a few recent theoretical works 29,55,56 on the tuning of magnetic interactions in Kitaev quantum magnets by light, which have some overlap with our current work. In refs. 55,56 , the authors treat the ligands by an effective hopping between two transition metal sites, and then employ Floquet theory for the two-site problem. In this treatment, the ratio between K and Γ interaction is fixed and cannot be tuned by laser. The effects of time-reversal symmetry breaking and hence the inverse Faraday effect is absent because the direct hopping of electrons between two sites does not produce a net hopping phase associated with the circularly polarized light. Also, ref. 29 investigated the model with all t 2g orbitals and ligand on a few site cluster using brute force exact diagonalization and investigated the spin Hamiltonian parameters numerically. In our work, instead, we use the formalism developed in pioneering work by Rau et al. 26 and solve two-site problem, in which e g orbitals and J eff = 3/2-manifold of t 2g are integrated out. This model is usually used as a starting point to understand the Iridate and Ruthenates. We have added ligands to that formalism, the details of which are presented in Supplementary Note 4. Integrating out the J eff = 3/2-manifold degrees of freedom makes our problem relatively tractable and allows us to provide analytical forms for the various interactions in the spin Hamiltonian, which are valuable to see the dependence of magnetic interactions on drive frequency and amplitude explicitly.
In conclusion, we show that light can tune magnetic interactions in Kitaev quantum magnets by explicitly calculating an effective spin Hamiltonian from the multi-orbit spin-orbit coupled Hubbard model in the presence of a circularly polarized light. We employ both the exact diagonalization of the Floquet Hamiltonian and analytical perturbation theory. We demonstrate that magnetic interactions favoring the QSL phase can be achieved by tuning the frequency and amplitude of the light. Our work points to a promising route to stabilize quantum spin liquid by coupling quantum magnets to light.

Methods
Floquet theory. We use Bloch wave theory, jψðtÞi ¼ e Àiϵ α t jψ α ðtÞi and solve the Schrodinger equation given by HðtÞjψðtÞi ¼ i ∂ ∂t jψðtÞi 7 . Using the Fourier transform given by H m ¼ 1 T R T 0 e imωt HðtÞdt and jψ α;m i ¼ 1 T R T 0 e imωt jψ α ðtÞi, we can solve the Schrodinger equation, where the solution is given by, ϵ α þ mω À Á jψ α;m i ¼ ∑ m 0 H mÀm 0 jψ α;m 0 i. Here, α are the basis states for the singly occupied, doubly occupied states on the TM and the singly occupied states on the ligand. The results for the quasi-energies converge quickly for off-resonance and in our case |m| < 8, is sufficient to get converged results. We use m = 0-Floquet sector to estimate the parameters in the spin model.
For the system with ligands, the second-order results are the same as Eq. (8). The third-order corrections are evaluated using Here, L j i (E L ) are the eigen-states (energies) for the occupancy on the ligand sites. T pd is the hopping between the TM and the ligand sites. The parameters, J, K, Γ, and h evaluated using this formalism are given in Eq. (7a-d).

Data availability
The data that support the findings of this study are available from the corresponding authors on reasonable request.

Code availability
Codes used to produce the findings of this study are available from the corresponding author upon reasonable request.