1 Introduction

The implementation of quantum processors requires precise control of quantum states exhibiting long decoherence times \(T_2\) in order to perform quantum algorithms or construct quantum memories. However, interactions with environmental degrees of freedom limit the coherence time and the qubit operability. In the case of spin qubits in solid state materials, examples of decoherence sources are spin–spin interactions, spin diffusion, fluctuations in the magnetic field or charge/electric noise [1, 2]. These interactions lead to a pure dephasing time \(T_\phi \) of the qubit much shorter than its energy relaxation time \(T_1\).

A number of techniques have been developed recently aiming to replace the pure dephasing with an actively controlled decoherence time by means of Electron Spin Resonance (ESR) pulses such as Dynamical Decoupling (DD) [3,4,5,6,7] and Concatenated Dynamical Decoupling (CDD) [8, 9]. The former uses a sequence of spin flips to improve spin coherence via a spin-echo technique, while the latter is using additional layer(s) of DD to solve for the imperfections and fluctuations of the first layer. The methods make use of strong \(\pi \)-pulses which complicates the operations due to gate imperfections; similar issues are noted in the case of a proposal using strong continuous waves [10, 11]. In general, DD techniques applied to electronic spin qubits have been used on nitrogen vacancy (NV) centers and with relative success [8, 12,13,14].

However, other spin qubit implementations with different anisotropic/isotropic properties [15,16,17], spin levels [18], electro-nuclear configurations [19], and host materials need to be studied for their potential use in quantum computing and quantum memories. Rare earth ions such as CaWO\(_4\):Er\(^{3+}\) [20], Y\(_2\)SiO\(_5\):Er\(^{3+}\) [21] and Y\(_2\)SiO\(_5\):Yb\(^{3+}\) [22] or \(^{28}\)Si:Bi [23] show long relaxation times as well as significant coherence times. Therefore, a universal method is needed, applicable to any type of qubit (spin or superconducting qubit for instance), using weak pulses which do not alter a quantum algorithm.

Generally speaking, a Rabi nutation opens up a gap in the spectrum of quasi-energies (Hamiltonian in the rotating frame) creating a or the ? so-called sweet spot where the Rabi rotation is insensitive to fluctuations of the static field detuning. In other words, the Rabi rotation protects itself, given that the operation is done with no detuning with one [24] or more photons [25]. In Ref. [26] we demonstrated a protocol based on Floquet resonances [27] where such a dynamical sweet spot is created for any detuning of the qubit away from resonance (see also a related approach implemented on NV centers [28, 29]). The method uses two microwave drives acting in parallel on the quantum state of a qubit and able to increase the coherence time under drive up to the relaxation time within the experimental conditions of Ref. [26]. The main drive induces imperfect Rabi oscillations while a circularly polarized second drive, of about two orders of magnitude smaller in amplitude, serves to sustain the oscillations. Rather than decoupling the qubit from the bath using a strong excitation, we use very weak pulses and alter the dynamics of the entire system. In practice, the time interval between two regular pulsed gates could contain an integer number of Rabi periods protected by our protocol and therefore the quantum information would be preserved from one gate to the next. This method was demonstrated in Ref. [26] on any initial state of spin systems with different anisotropies, spin sizes and spin-orbit couplings. The protocol is universal and it should be applicable to other qubit implementations such as nuclear spins or superconducting qubits.

In this paper, the experimental protocol described in Ref. [26] is studied analytically and by computer simulation. The aim of this paper is to scrutinize the microscopic physical mechanisms that may cause the sustained, long-time Rabi oscillations to appear and to give a consistent microscopic description of the experimental data displayed in Fig. 1.

Fig. 1
figure 1

Experimental data of the Rabi oscillations in CaWO\(_4\):Gd\(^{3+}\) at \(T=40\,\mathrm {K}\) for different values of the detuning \(\Delta \), phase differences \(\phi =0\) and \(\phi =45^\circ \), and fixed amplitudes \(h_d=15\,\hbox {MHz}\) and \(h_i=0.12h_d=1.8\,\hbox {MHz}\), respectively

Full size image

We consider three different models for the source of decoherence. The first model describes a spin qubit, that is a magnetic moment, subject to external magnetic fields and interacting with a bath of two-level systems. These two-level systems are not directly affected by the external magnetic fields. This kind of bath, which we refer to as bath-I, mimics an environment consisting of two-level systems representing e.g., nuclear spins, electronic spins of different species, defects, etc.,

The second model describes an ensemble of identical magnetic moments, all driven by the same microwave fields and interacting with each other via all-to-all, dipolar-like weak couplings. The effect of these couplings on any of the magnetic moments is to cause decoherence. Thus, from the viewpoint of an arbitrarily chosen magnetic moment, all others may be regarded as representing a kind of “bath”, which we refer to as bath-II. The third model extends the second one by taking the inhomogeneity of the microwave fields into account.

The paper is organized as follows. In Sect. 2, we present and discuss new electron spin resonance results for CaWO\(_4\):Gd\(^{3+}\). Sect. 3 gives an overview of the three different models that are at the basis of our simulation work. In Sect. 4, we specify the bath-I model in detail and present the simulation results. Section 5 introduces the Hamiltonian of the bath-II model and discusses the simulation results for the case without and with the inhomogeneity of the microwave fields. In Sect. 7, we present our conclusions. Technical details and additional results can be found in the Supplemental Material [30].

2 Electron spin resonance experiments

Experiments probing the magnetization dynamics under external driving fields such as the ones performed in Ref. [26], pose interesting theoretical problems. In these electron spin resonance (ESR) experiments, one measures a signal that is proportional to the expectation value of the z-component of the magnetization where the z direction is defined by the direction of the static magnetic field \({\mathbf {B}}_0\). In symbols:

$$\begin{aligned} \mathrm {Signal}\propto & {} \sum _{n=1}^{N}\langle S^z_n(t) \rangle , \end{aligned}$$
(1)

where N is the number of magnetic moments in the sample. The presence of a high frequency resonant microwave field perpendicular to the static field induces Rabi oscillations. A characteristic feature of Rabi oscillations is that their frequency changes linearly with the applied microwave power. These oscillations decay due to various dissipation effects.

In the ESR experiments under scrutiny, Rabi oscillations are induced by applying an electromagnetic drive pulse of strength \(h_d\) and frequency \(f=f_0+\Delta \) where \(f_0\) is the Larmor frequency and \(\Delta \) is a (small) detuning. Simultaneously with the drive pulse, an image pulse of strength \(h_i\) and of frequency \(f=f_0-\Delta \), coherent with the drive pulse, is applied. The power of the image drive is much less than that of the drive pulse itself, i.e. \(h_i\ll h_d\). The phase difference \(\phi \) between the drive and image pulse can be used to manipulate the Rabi oscillations. Under optimal image drive conditions, experiments show sustained Rabi oscillations up to the maximum amplifier gate length of \(15\,\mu \hbox {s}\) [26].

Beside data previously reported in Ref. [26], we present results of new experiments in Fig. 1, performed under the same conditions as for the experiments earlier. These data have been obtained for Gd\(^{3+}\) moments in a CaWO\(_4\) host lattice. Although these moments have spin \(S=7/2\), the orientation of the static magnetic field and the frequency and power of the microwave pulses are chosen to select transitions between one pair of Zeeman levels only [26]. Therefore, the spin system may be viewed as an effective two-level system performing Rabi oscillations. In this paper, we take the CaWO\(_4\):Gd\(^{3+}\) system for our case study and use the results presented in Fig. 1 as reference for comparison with the model calculations. In this paper, we use Fig. 1 as a reference for comparing with our simulation results for the three different models.

The main conclusion drawn on the basis of the experimental data [26] and new experimental data presented in Fig. 1 may be summarized as follows:

  1. 1.

    Regarding the effect of the image pulse on the decay of the Rabi oscillations, the experimental findings are generic, meaning that the main features do not significantly depend on the kind of impurity that consitutes the qubit, the presence or absence of nuclear spins, etc. See Ref. [26] for more details.

  2. 2.

    In the absence of the image pulse, the amplitudes of the Rabi oscillations vanish rapidly, on a time scale of \(1\,\upmu \mathrm {s}\) [26].

  3. 3.

    For a detuning \(\Delta =5\,\mathrm {MHz}\), the Rabi oscillations decay rapidly, on a time scale of \(1\,\upmu \mathrm {s}\) see Figs. 1a, f.

  4. 4.

    For a detuning \(\Delta =7\,\mathrm {MHz}\), we observe either Rabi oscillations that decay rapidly, on a time scale of \(1\mu \mathrm {s}\) if \(\phi =0\) (see Fig. 1b) or sustained Rabi oscillations if \(\phi =45^\circ \) (see Fig. 1g).

  5. 5.

    For a detuning equal of \(\Delta =8.66\,\hbox {MHz}\) corresponding to the first Floquet resonance frequency (see below), we observe sustained Rabi oscillations with some signatures of beating, Fig. 1c, h, depending on the value of the phase \(\phi \) between the drive and image pulse, in qualitative agreement with earlier experiments on MnO:Mn\(^{2+}\) [26].

  6. 6.

    Increasing the detuning further to \(\Delta =10,12\,\hbox {MHz}\), see Fig. 1d, e, i, j, we observe sustained Rabi oscillations with small amplitudes.

  7. 7.

    Figure 1c, d, g show that the appearance of sustained Rabi oscillations depends on both \(\Delta \) and \(\phi \). Therefore this appearance cannot be solely attributed to the existence of a Floquet resonance, the frequency of which depends on \(\Delta \) and \(h_d\) and is independent of \(\phi \) (see below).

In summary, from the experimental results shown in Fig. 1, we conclude that there are two distinct features. Depending on the value of the detuning further to \(\Delta \) (for constant drive and image drive amplitudes \(h_d\) and \(h_i\), respectively), sustained long-time Rabi oscillations appear. This is the main feature. In contrast to the usual decaying Rabi oscillations observed in the absence of the detuning (\(\Delta =0\)), in some situations the Rabi oscillations may persist for a long time. The second main feature is the additional structure of these sustained oscillations which depends on the detuning \(\Delta \) and the phase shift \(\phi \).

3 Overview of the microscopic models

The magnetic moments being probed in the ESR experiments are distributed in the host lattice in a highly diluted manner [26]. Thus, the most basic model is that of non-interacting spins subject to a time-dependent magnetic field, see Sect. 4. In this paper, we study a model in which the contributions to the magnetic field are (1) a strong static magnetic field, (2) a microwavedrive field and (3) the novel aspect, the microwave image field [26]. The pure quantum dynamics of this model is straightforwardly studied by solving the time-dependent Schrödinger equation for a single spin, see Sect. 4. To account for the explicit time-dependence of the system + bath Hamiltonian, we use a second-order decomposition formula for ordered matrix exponentials [31, 32] to solve the TDSE.

Obviously, there is no decoherence or dissipation in this most basic model. In the presence of the image drive, the magnetization dynamics exhibits modulated Rabi oscillations. The amplitude of the modulation depends on the detuning \(\Delta \) and the phase shift \(\phi \).

In the present paper, we study the decay of the Rabi oscillations by solving for the pure quantum dynamics of models with many spins in which decoherence and dissipation appear automatically. Decoherence and dissipation effects may also be studied through quantum master equations [33]. This approach makes several implicit assumptions and often involves many parameters [33]. We found it difficult to find a set of parameters which reproduce, even qualitatively, the dependence on the detuning \(\Delta \) and the phase shift \(\phi \) observed in experiments. In section IV of the Supplemental Material [30], it is shown that, even when the inhomogeneity of the microwave field is included, this method does not describe all the qualitative features of the experimental data shown in Fig. 1. Therefore, in this paper we focus on microscopic many-body models and explore the conditions under which these models reproduce the systematic trends of the features observed experimentally.

As alluded to earlier, in real samples, the microscopic mechanisms accounted for separately in the present paper are most likely simultaneously active. Solving the time-dependent Schrödinger equation for 28 mutually interacting spins, the largest systems we have studied systematically, already requires substantianal computational resources. Adding another 28 two-level systems to represent a bath requires computational resources which, at present, are prohibitive. Performing such simulations may be a challenging project for the future.

In the next three sections, we scrutinize theoretical models by visually comparing the outcomes of numerical simulations with the experimental results shown in Fig. 1.

4 Single-spin system

This section introduces the model for the magnetic moment of a single Gd\(^{3+}\) ion subject to external magnetic fields. To good approximation (for our purposes), this model describes a two-level system driven by a periodic field. Although too simple to serve as a model for the system studied experimentally, its dynamics is already complicated. For instance, because of the presence of a driving fields \(h_d\) and \(h_i\), the model can exhibit resonances and beating oscillations. For pedagogical purposes, this section presents a few results of the dynamics of the single-spin model in time-dependent magnetic fields.

The system (S), contains one magnetic moment in an external time-dependent magnetic field and is defined by the Hamiltonian

$$\begin{aligned} H_{\mathrm {S}}= & {} -\omega _0 S^z - 2\omega _d S^x\sin [2\pi (f_0+\Delta ) t + \phi ] \nonumber \\&- 2\omega _i S^x\sin [2\pi (f_0-\Delta ) t - \phi ] , \end{aligned}$$
(2)

where \({\mathbf {S}}=(S^x,S^y,S^z)\) denotes the three components of magnetic moment. We use units such that \(\hbar =1\). To facilitate the comparison with experimental data, we express frequences such as \(f_0=\omega _0/2\pi \), \(\Delta \), \(h_d=\omega _d/2\pi \), and \(h_i=\omega _i/2\pi \) in MHz. The first term in Eq. (2) describes the Larmor rotation (frequency \(f_0\)) of the spin due the static magnetic field \({\mathbf {B}}_0\). The second and third term in Eq. (2) are called the microwave drive and image field, respectively.

In ESR experiments, \(\omega _0\) typically is several orders of magnitude larger than \(\omega _d\) and \(\omega _i\). Instead of adopting the standard resonance condition, we choose \(\omega =\omega _0+2\pi \Delta =2\pi (f_0+\Delta )\), see Ref. [26], and use \(U=\exp \left( i\omega t S^z\right) \) to find (for details on the calculation, see section I of the Supplemental Material [30])

$$\begin{aligned} H_\mathrm {S,RF}(t)\equiv & {} -iU^\dagger \frac{\partial }{\partial t}U+U^\dagger H_\mathrm {S} U \nonumber \\= & {} 2\pi \big \{ \Delta S^z - h_d \left( S^x \sin \phi + S^y\cos \phi \right) \nonumber \\&+ h_i \left[ S^x \sin (4\pi t\Delta +\phi ) \right. \nonumber \\&\left. - S^y \cos (4\pi t \Delta + \phi ) \right] \big \} \nonumber \\= & {} -{\mathbf {B}}(t)\cdot {\mathbf {S}}\;, \end{aligned}$$
(3)

where we have omitted terms that oscillate with the (very large) frequency \(2f_0\) and introduced the time-dependent magnetic field

$$\begin{aligned} {\mathbf {B}}(t)=2\pi \left( \begin{array}{c} h_d\sin \phi -h_i\sin (4\pi t\Delta +\phi )\\ h_d\cos \phi +h_i\cos (4\pi t\Delta +\phi )\\ \Delta \\ \end{array} \right) \;. \end{aligned}$$
(4)

If \(h_i=0\), \({\mathbf {B}}(t)\) is constant in time and it follows immediately that the spin \({\mathbf {S}}\) will perform oscillations with the Rabi frequency \(F_{\mathrm {R}}=(\Delta ^2+h_d^2)^{1/2}\). Furthermore, if \(\Delta =0\), the magnetization \(\langle S^z(t)\rangle \) displays the usual Rabi oscillations, that is the spin performs rotations about the vector \((h_d\sin \phi ,h_d\cos \phi ,0)^{{\mathbf {T}}}\).

In sect. 1 of the Supplemental Material [30], we show that if \(2\Delta =F_\mathrm{R}\) (or equivalently \(\Delta =h_d/\sqrt{3}\)) and \(h_i\not =0\), the spin performs a second rotation with a frequency \(F_\mathrm{R}^{(2)}=3h_i/4\). In other words, the condition \(2\Delta =F_\mathrm{R}\) defines a special point in the parameter space of the model defined by Eq. (3).

Fig. 2
figure 2

Simulation results obtained by solving the TDSE for the Hamiltonian Eq. (3), describing a single spin subject to the time-dependent magnetic field \({\mathbf {B}}(t)=2\pi (h_d\sin \phi -h_i\sin (4\pi t\Delta +\phi ), h_d\cos \phi +h_i\cos (4\pi t\Delta +\phi ),\Delta )^{\mathrm {T}}\), with the amplitudes \(h_d=15\,\hbox {MHz}\) and \(h_i=0.12 h_d=1.8\,\hbox {MHz}\). For additional information, see sections I–II of the Supplemental Material [30]

Full size image
Table 1 The Rabi frequency and the six smallest (in absolute value) quasi-energies obtained by solving the Floquet problem Eq. (S18) for the subspace spanned by \(\{\psi _{-100},\ldots ,\psi _{-1},\psi _{0},\psi _{1},\) \(\ldots ,\psi _{100}\}\) for different values of the offset \(\Delta \), \(h_d=15\,\hbox {MHz}\) and \(h_i=1.8\,\hbox {MHz}\)
Full size table

As \({\mathbf {B}}(t)={\mathbf {B}}(t+1/2\Delta )\), the Hamiltonian Eq. (3) is a periodic function of time. Some insight into the dynamics of the periodically driven system described by Eq. (3) can be obtained by resorting to Floquet theory [34, 35]. It is easy to show (for details see section III of the Supplemental Material [30]) that in the perturbation regime \(\vert h_i\vert \ll h_d\), Floquet theory predicts the existence of resonances if the parameters \(\Delta \) and \(h_d\) satisfy the condition

$$\begin{aligned} \Delta =\frac{h_d}{\sqrt{4n^2-1}}, \quad n=1,2,\ldots \;. \end{aligned}$$
(5)

For \(n=1\), Eq. (5) yields \(\Delta =h_d/\sqrt{3}\) and therefore we refer to \(2\Delta =F_\mathrm{R}\) as the condition for the first Floquet resonance. It is essential to note that the Floquet resonance opens up a gap or level repulsion for the quasi-energies in resonance which can be seen as a dynamical sweet spot: the levels are flat and thus insensitive to noise in \(h_d\) or \(\Delta \). Our numerical simulations uncover the microscopic, physical mechanisms developing inside the bath (resonant or non-resonant to the drive).

Experimentally reasonable values for the CaWO\(_4\):Gd\(^{+3}\) system are \(h_d=15\,\hbox {MHz}\) and \(h_i=0.12h_d=1.8\,\hbox {MHz}\) [26] and the condition for the first Floquet resonance reads \(\Delta =h_d/\sqrt{3}=8.66\,\hbox {MHz}\) with \(F_{\mathrm {R}}=2\Delta =17.33\,\hbox {MHz}\) [26].

The dynamics of the single-spin system Eq. (3) is studied by solving the TDSE with Hamiltonian Eq. (3) numerically. In Fig. 2, we present simulation results for the Rabi oscillations for different values of the detuning \(\Delta \) and two values of the phase shift \(\phi \) (see section II of the Supplemental Material [30] for the corresponding pictures in the frequency domain).

As expected, the effect of the image drive is most pronounced if the parameters match the condition for the first Floquet resonance. Then, the magnetization exhibits considerable beating, a manifestation of the presence of the second Rabi frequency \(F_\mathrm{R}^{(2)}=3h_i/4\ =1.35\,\hbox {MHz}\).

The Fourier transform of the data depicted in Fig. 2c (see Fig. S1(c) of the Supplemental Material [30]) shows a peak at the Rabi frequency \(F_\mathrm{R}=17.33\,\hbox {MHz}\) and a lower/upper sideband at \(F_\mathrm{R}\mp \delta f\) with \(\delta f=1.34\,\hbox {MHz}\). The spectral weight of these sidebands is about 8% of the main signal. In addition, there is a weak signal (spectral weight \(\approx 0.8\,\%\)) at zero frequency with side bands (spectral weight \(\approx 1\,\%\)) at \(\delta f=1.34\,\hbox {MHz}\;\approx F_\mathrm{R}^{(2)}=1.35\,\hbox {MHz}\).

From Table 1, third row, if follows that \(F_\mathrm{R}=\kappa _2+\kappa _1 = 17.32 \,\hbox {MHz}\) and \(\delta f\approx \kappa _2-\kappa _1 = 1.34 \,\hbox {MHz}\). The frequency of the modulation in the time dependent signal Fig. 2c is approximately \(1.34\,\hbox {MHz}\), the same as to the difference between the Rabi frequency and the frequency the sidebands. Thus, we can relate the Rabi frequency and the sideband frequencies to the quasi-energies obtained from Floquet theory and the modulation of the time dependent signal. Moreover, the sideband frequencies are, to a very good approximation, given by \(F_\mathrm{R}^{(2)}=1.35\,\hbox {MHz}\).

Fig. 3
figure 3

Bloch sphere representation of the spin motion for the cases shown in Fig. 2a–e. The brightness of the trajectory changes from dark to light as time proceeds

Full size image
Fig. 4
figure 4

Same as Fig. 3 except that the pictures corresponds to the cases shown in Fig. 2h–j

Full size image

From Fig. 2 and the pictures of the motion of the magnetization on the Bloch sphere shown in Figs. 3 and 4, it is clear that the presence of an image drive can affect the Rabi oscillations in a nontrivial manner. The salient features of the magnetization dynamics, that is the Rabi oscillations and the frequency of the amplitude modulation (if present) of these oscillations caused by the image drive can be related to quasi-energies obtained from Floquet theory, see section II of the Supplemental Material [30]). The key point, advanced in Ref. [26] and illustrated by Fig. 1, is that the image drive can be used to significantly reduce the decay time of the Rabi oscillations.

A simple, phenomenological approach to include effects of dissipation and decoherence is to resort to a description in terms of the Markovian quantum master equation. Section IV of the Supplemental Material [30], shows that also this approach, even when the inhomogeneity of the microwave field is included, does not describe all the qualitative features of the experimental data shown in Fig. 1.

5 Single-spin system interacting with a bath of two-level systems

In this section, we study a microscopic model in which the spin qubit, the magnetic moment of the Gd\(^{+3}\) ion referred to as system (S), interacts with a bath (B) of pseudo-spins. In this section, the spins other than that of the Gd\(^{+3}\) ion are the bath spins. We analyze the spin dynamics of this many-body \(S=1/2\) system by solving the TDSE numerically. Here, \({\mathbf {S}}=\varvec{\sigma }/2\) where \(\varvec{\sigma }=(\sigma _x,\sigma _y,\sigma _z)=(\sigma _1,\sigma _2,\sigma _3)\) are the Pauli matrices.

The Hamiltonian of the system (S) + bath (B) takes the generic form

$$\begin{aligned} H(t)= & {} H_\mathrm {S}(t) + H_\mathrm {B} + \lambda H_\mathrm {SB} , \end{aligned}$$
(6)

where \(H_\mathrm {B}\) and \(H_\mathrm {SB}\) are the bath and system-bath Hamiltonians, respectively. The overall strength of the system + bath-I interaction is controlled by the parameter \(\lambda \). The Hamiltonian for the system-bath interaction is chosen to be

$$\begin{aligned} H_\mathrm {SB}= & {} -\sum _{n=1}^{N_{\mathrm {B}}} \;\;\sum _{\alpha =x,y,z}J_{n}^\alpha I_{n}^\alpha S^\alpha , \end{aligned}$$
(7)

where \(N_{\mathrm {B}}\) is the number of spins in the bath, the \(J^\alpha _{n}\) are uniform random frequencies in the range \([-J,+J]\) and \(I_{n,k}\) is the k-th component if the bath spin \({\mathbf {I}}_n\). As the system-bath interaction strength is controlled by \(\lambda \), we may set \(J=1\,\hbox {MHz}\). As the bath-I Hamiltonian, we take

$$\begin{aligned} H_\mathrm {B}= & {} -\sum _{n=1}^{N_{\mathrm {B}}}\;\;\sum _{\alpha =x,y,z} K_{n}^\alpha I_{n}^\alpha I_{n+1}^\alpha , \end{aligned}$$
(8)

where the \(K_{n,k}\)’s are uniform random frequencies in the range \([-K,K]\). In our simulation work, we use periodic boundary conditions \(I_{n}^\alpha =I_{n+N}^\alpha \). The Hamiltonian Eq. (8) describes a collection of magnetic moments, located on a ring of lattice sites, and interacting with their nearest neighbors. Because of the random couplings, it is unlikely that Eq. (8) is integrable (in the Bethe-ansatz sense) or has any other special features such as a conserved magnetization etc.

Expressing the motion of the spin qubit in the rotating frame has no effect on the bath Hamiltonian Eq. (8) because the transformation to the rotating frame only affects the spin \({\mathbf {S}}\), not the operators describing bath-I. However, expressing the motion of the spin qubit in the rotating frame changes the system-bath Hamiltonian Eq. (7) to

$$\begin{aligned} H_\mathrm {SB}= & {} -\sum _{n=1}^{N_{\mathrm {B}}} J_{n}^z I_{n}^zS^z -\sum _{n=1}^{N_{\mathrm {B}}} J_{n}^x I_{n}^x\left( S^x \cos \omega _0 t\right. \nonumber \\&\left. +S^y \sin \omega _0 t\right) \nonumber \\&-\sum _{n=1}^{N_{\mathrm {B}}} J_{n}^y I_{n}^y\left( S^y \cos \omega _0 t-S^x \sin \omega _0 t\right) \;. \end{aligned}$$
(9)

Discarding the terms that oscillate with the high angular frequency \(\omega _0\), only the coupling between the z-component of the system and bath spins survives, yielding

$$\begin{aligned} H_\mathrm {SB,RF}= & {} -\sum _{n=1}^{N_{\mathrm {B}}} J_{n}^z I_{n}^zS^z . \end{aligned}$$
(10)

Therefore, within this approximation, the Hamiltonian of the system (S) + bath (B) in the rotating frame reads

$$\begin{aligned} H(t)= & {} H_\mathrm {S,RF}(t) + H_\mathrm {B} + \lambda H_\mathrm {SB,RF} . \end{aligned}$$
(11)

As \(H_\mathrm {S,RF}(t)\) and \(H_\mathrm {SB,RF}\) do not commute, the system and bath-I can exchange energy if \(\lambda \not =0\).

As \(J^\alpha _n\) is chosen uniformly random in the interval \([-1\,\hbox {MHz},+1\,\hbox {MHz}]\), the average of strength of the system-bath interaction is \(\lambda \langle \vert J^\alpha _n\vert \rangle =\lambda /2 \,\hbox {MHz}\). On the other hand, the characteristic energy scale of the spin qubit is \(2\pi h_d\approx 94\,\hbox {MHz}\). Therefore, for \(\lambda <10\), this strength is at least an order of magnitude smaller than the characteristic energy scale of the spin qubit. Similarly, for \(K=10 \,\hbox {MHz}\) and \(\lambda =6\) (see Fig. 5)), \(\langle \vert K^\alpha _n\vert \rangle =5\,\hbox {MHz}\), of the same order of magnitude as \(\lambda \langle \vert J^\alpha _n\vert \rangle =3\,\hbox {MHz}\). In the absence of the image drive, the characteristic decay time of the Rabi oscillations is \(1\,\upmu \mathrm {s}\) or less [26]. Performing a simulation with \(h_i=0\) and \(\Delta =8.66\,\hbox {MHz}\) (data not shown) we find that the decay time of the Rabi oscillations is about \(0.3\,\upmu \mathrm {s}\). Thus, our choice for the model parameters yields a decay time of the Rabi oscillations which, in the absence of the image drive, is short as observed experimentally for very different samples [26].

As explained in more detail in section V of the Supplemental Material [30], to study the Rabi oscillations, the initial state of the whole is taken to be

$$\begin{aligned} \vert \Psi (t=0)\rangle= & {} \vert {}\uparrow \rangle \otimes \vert \Phi \rangle \;, \end{aligned}$$
(12)

where \(\Phi \) is a random vector in the \(N_{\mathrm {B}}\)-dimensional Hilbert space.

Evidently, the energy of the whole system, which is closed, is a conserved quantity. However, the system S and the bath B are in contact with each other and can exchange energy. Thus, the system S can show decoherence and relaxation, equilibration, thermalization, all depending on how the simulation is performed [36, 37].

Fig. 5
figure 5

Simulation results for the magnetization, as obtained by solving the TDSE for Hamiltonian Eq. (11), describing the system S coupled to bath-I, both subject to a time-dependent magnetic field. Bath-I contains \(N_{\mathrm {B}}=28\) interacting spins, the spin-bath coupling \(\lambda =6\), \(J=1\,\hbox {MHz}\), \(K=10\,\hbox {MHz}\), \(h_d=15\,\hbox {MHz}\), and \(h_i=0.12\;h_d=1.8\,\hbox {MHz}\)

Full size image

Figure 5 shows the simulation data obtained by solving the TDSE for Hamiltonian Eq. (11) with a bath of \(N_{\mathrm {B}}=28\) spins. Qualitatively, model Eq. (11) reproduces all essential features of the experimental data shown in Fig. 1.

In model Eq. (11), there are only two adjustable parameters, namely the system-bath interaction strength \(\lambda \) and the scale of the inter-bath interactions K. In view of the fact that solving the TDSE for a system with 29 spins is quite expensive in terms of computer resources, we did not try to improve the qualitative agreement with the experimental results by adjusting these parameters.

It should be noted that being exactly at the Floquet resonance frequency \(\Delta =8.66\,\hbox {MHz}\) is not a necessary condition to observe sustained Rabi oscillations: they also appear in Fig. 5b, g (and Fig. 1b, g).

Fig. 6
figure 6

Simulation results for the energy of the spin bath (dashed red line), the system (dotted blue line) and total energy (solid black line), as obtained by solving the TDSE for Hamiltonian Eq. (11), describing the single-spin system S coupled to a spin bath, both subject to a time-dependent magnetic field. The bath consists of \(N_{\mathrm {B}}=28\) spins, the spin-bath coupling \(\lambda =6\), \(J=1\,\hbox {MHz}\), \(K=10\,\hbox {MHz}\), \(h_d=15\,\hbox {MHz}\), \(h_i=0.12h_d=1.8\,\hbox {MHz}\), and \(\Delta =8.66\,\hbox {MHz}\). The lines for the system and total energy nearly overlap. For \(t\ge 4\,\upmu \mathrm {s}\), the energies oscillate on a scale that is hardly visible in the plots. Results for non-resonant values of \(\Delta =8.66\,\hbox {MHz}\) are presented in Fig. S9

Full size image

In Fig. 6 we present TDSE simulation data of the energy of the system \(E_\mathrm {S}=\langle \Psi (t)\vert H_\mathrm {S,RF}\vert \Psi (t)\rangle \), of bath-I \(E_\mathrm {B}=\langle \Psi (t)\vert H_\mathrm {B,RF}\vert \Psi (t)\rangle \), and the total energy \(E_\mathrm {S}+E_\mathrm {B}+\lambda E_\mathrm {SB}\langle \Psi (t)|H(t)|\Psi (t)\rangle \). From a more general, many-body physics viewpoint, these data show interesting behavior.

It is not evident that the qubit will evolve to a stationary state if it is driven by a time dependent field (of strength \(h_i\)) [38, 39]. Recall that for \(\lambda =6\), the qubit-bath interaction is rather weak. Nevertheless, as Fig. 6 shows, in the case under study, it is clear that as time progresses, both the qubit and bath-I end up in a stationary state, the details of which depend on the offset \(\Delta \) and the phase shift \(\phi \). Furthermore, if the value of \(\Delta \) matches the resonance condition \(\Delta =h_d/\sqrt{3}\), the system and total energy decay to the stationary state on a much longer time scale than for \(\Delta =5,7,10,12\,\hbox {MHz}\) (see Fig. S9 of the Supplemental Material [30]).

The period of the large-amplitude oscillations in Fig. 6 is approximately \(0.7\,\upmu \hbox {s}\approx P^{(2)}_\mathrm {R}=0.74\,\upmu \hbox {s}\). Thus, at resonance (\(\Delta =h_d/\sqrt{3}\)), the system and total energy exhibit a synchronized, damped oscillatory behavior with a frequency that is approximately given by the second Rabi frequency \(F^{(2)}_\mathrm {R}=3h_i/4\).

In Fig. 7 we show the correlation function

$$\begin{aligned}&B(t)\nonumber \\&=\frac{1}{N_{\mathrm {B}}} \sum _{\alpha =x,y,z}\sum _{n=1}^{N_{\mathrm {B}}}\; \langle \Psi (t=0)\vert I_{n}^\alpha (t) I_{n+1}^\alpha (t)\vert \Psi (t=0) \rangle \;,\nonumber \\ \end{aligned}$$
(13)

that is, the sum of all equal-time, nearest-neighbor correlations of all bath-I pseudo spins components at the Floquet resonance \(\Delta =8.66\,\mathrm {MHz}\). A first observation is that these correlations are small (relative to the maximum of B(t) which is equal to one). A second observation, less clear than in the case \(h_i=1.8\,\mathrm {MHz}\) and \(\phi =0\), is that as time proceeds, the correlations seem to saturate. The transient dynamics includes signatures of the Floquet and Rabi dynamics, because the spin qubit is driving bath-I via the \(H_{SB}\) interaction. Also, it’s important to note that a saturation of the bath-I internal dynamics can lead to an increase of the coherence time, as it was observed in the case of molecular magnets placed in high magnetic fields [40]. The behavior of the bath-I correlation B(t) is very different from the that of the correlations in the model of interacting system spins, discussed in the next section.

Fig. 7
figure 7

Simulation results of correlation Eq. (13) obtained by solving the TDSE for Hamiltonian Eq. (11) describing the single-spin system S coupled to a spin bath, both subject to a time-dependent magnetic field. The bath consists of \(N_{\mathrm {B}}=28\) interacting spins

Full size image

In summary: the model Eq. (11) reproduces the main features (see Sect. 1) of the spin dynamics observed in the ESR experiments on CaWO\(_4\):Gd\(^{3+}\).

6 System of interacting magnetic moments

In the previous section, we studied effects of environments on the spin dynamics and found decaying and sustained Rabi oscillations. Regarding mechanisms for the magnetization decay, as a alternative to the coupling to environment, we may also consider interactions among the Gd\(^{3+}\) moments and the inhomogeneity of the microwave field. In this section, we study these two aspects by considering a collection of spins, each one representing the magnetic moment of a Gd\(^{3+}\) ion (and without additional spins representing a bath). The model describes \(N_\mathrm {S}\) spins that interact with the external magnetic fields and with each other.

The Hamiltonian of the whole system takes the form

$$\begin{aligned} H(t)= & {} \sum _{n=1}^{N_\mathrm {S}} \big \{-\omega _0 S^z_n - 2\omega _d S^x_n\sin [2\pi (f_0+\Delta ) t + \phi ] \nonumber \\&\text{ t }o 1cm{} -2\omega _i S^x_n\sin [2\pi (f_0-\Delta ) t - \phi ] \big \} \nonumber \\&-\lambda \sum _{1=m<n}^{N_{\mathrm {S}}}\;\;\sum _{\alpha =x,y,z} C_{m,n,\alpha } S_{m}^\alpha S_{n}^\alpha \;, \end{aligned}$$
(14)

where the parameter \(\lambda \) is used to control the overall strength of the coupling between the spins. The first term in Eq. (14) describes the spins in the external time-dependent magnetic field. The interaction between the spins is described by the second term in Eq. (14). The coefficients \(C_{m,n,k}\)’s are taken to be uniform random frequencies in the range \([-1\,\hbox {MHz},1\,\hbox {MHz}]\). These randomly chosen \(C_{m,n,k}\)’ s are assumed to mimic, in a simple manner, the essential features of the dipolar interactions between the Gd spins.

As before, it is advantageous to eliminate the very fast motion associated with the Larmor angular frequency \(\omega _0\) by changing to the rotating frame. The unitary transformation that accomplishes this is \(U=\exp \left[ i\omega _0(S^z_1+\ldots +S^z_{N_{\mathrm {S}}})\right] \). Keeping only the secular terms, i.e. the terms that do not depend on \(\omega _0\), the interaction Hamiltonian does not change and Eq. (14) becomes

$$\begin{aligned} H(t)= & {} 2\pi \sum _{n=1}^{N_{\mathrm {S}}} \big \{ \Delta S^z_n - h_d \left( S^x_n \sin \phi + S^y_n\cos \phi \right) \nonumber \\&+ h_i \left[ S^x_n \sin (4\pi t\Delta +\phi ) - S^y_n \cos (4\pi t \Delta + \phi ) \right] \big \} \nonumber \\&-\lambda \sum _{1=m<n}^{N_{\mathrm {S}}}\;\;\sum _{\alpha =x,y,z} C_{m,n,\alpha } S_{m}^\alpha S_{n}^\alpha \;. \end{aligned}$$
(15)

As the C’s in Eq. (15) are chosen uniformly random in the interval \([-1\,\hbox {MHz},+1\,\hbox {MHz}]\), the average strength of the spin–spin interactions is \(\lambda \langle \vert C_{m,n,\alpha }\vert \rangle =\lambda /2 \,\hbox {MHz}\). On the other hand, the characteristic energy scale of the spins is \(2\pi h_d\approx 94\,\hbox {MHz}\). Therefore, for \(\lambda =5\) (the value adopted in our simulations), this strength is at least an order of magnitude smaller than the characteristic energy scale of the qubit. In other words, the coupling between the spins may be considered to be weak.

If we naively consider a uniform distribution of impurities over the host lattice with the concentration that we believe applies to the sample, we estimate the magnitude of \(\lambda C\) (see Eq. (17)) to be approximately 0.3 MHz. This is a factor of 10 smaller than the average interaction strength (2.5 MHz) used in our simulations. On the other hand, in our simulations, only a small (not more than 28 spins) cluster of spins suffices to create the “sustained Rabi oscillations”. Assuming that average distance between the Gd\(^{3+}\) ions in this cluster is only about two times smaller than for a uniform distribution Gd\(^{3+}\) ions, changes the magnitude of magnitude of \(\lambda C\) to approximately 3 MHz. A slightly different, clustered distribution of impurities would suffice. Clearly, there is a lot of uncertainty in these estimates but at least they are not off by orders of magnitude.

Another aspect, not explicitly included in our simulation models is the following. As a spin system, Gd\(^{3+}\) has \(2S+1=8\) levels, two of which are effectively used in the ESR experiment. It is well-known, e.g. in the field of superconducting qubits [41, 42], that multilevel structures contribute to decoherence of the two-level system regarded as the qubit. The two-level systems of the bath are a very simple model that can mimic this mechanism. Of course, two-level systems can mimic “almost anything”, as we also point out in the text. Also in this case, it is hard to put a reliable number on the interaction strength.

In the present case, to study the Rabi oscillations, we solve the TDSE with Hamiltonian Eq. (15) and take as the initial state the product state of all spins up (along the z-axis), that is

$$\begin{aligned} \vert \Psi (t=0)\rangle= & {} \vert {}\uparrow \rangle _1 \otimes \cdots \otimes \vert {}\uparrow \rangle _{N_{\mathrm {S}}} \;. \end{aligned}$$
(16)

Instead of the plotting the expectation value of the z-component of each spin, we now plot \(\langle M^z(t)\rangle =\sum _{n=1}^{N_{\mathrm {S}}} \langle \sigma ^z_n(t)\rangle \), that is the z-component of the total magnetization.

Figure 8 shows the simulation data obtained by solving the TDSE for Hamiltonian Eq. (15) with a system of \(N_{\mathrm {S}}=28\) interacting spins. Clearly, there is a significant qualitative difference between \(\phi =0\) and \(\phi =45^\circ \) data. Furthermore, the latter does not compare well with the data shown in Fig. 1(h–j). For completeness, Figs. S11 and S12 of the Supplemental Material [30] show the simulation data for the time evolution of the energy of the model Eq. (15).

6.1 Inhomogeneity of the microwave fields

A trivial modification to the Hamiltonian of the system of interacting spins allows us to account for the inhomogeneity of the microwave fields. As we now show, the combination of spin–spin interactions and the inhomogeneity of the microwave fields yields satisfactory results for the main features observed experimentally (for details see section IV C of the Supplemental Material [30]).

Instead of Eq. (15), we now take

$$\begin{aligned} H(t)= & {} 2\pi \sum _{n=1}^{N_{\mathrm {S}}}\big \{ \Delta S^z_n - h_{d,n} \left( S^x_n \sin \phi + S^y_n\cos \phi \right) \nonumber \\&+ h_{i,n} \left[ S^x_n \sin (4\pi t\Delta +\phi ) - S^y_n \cos (4\pi t \Delta + \phi ) \right] \big \} \nonumber \\&-\lambda \sum _{1=m<n}^{N_{\mathrm {S}}}\;\;\sum _{\alpha =x,y,z} C_{m,n,\alpha } S_{m}^\alpha S_{n}^\alpha \;. \end{aligned}$$
(17)

In words, we place each spin n in its own microwave fields \((h_{d,n},h_{i,n})\) where \(h_{d,n}/h_d\) and \(h_{i,n}/h_i\) are uniform random numbers in the range \([1-\delta ,1+\delta ]\), \(\delta \) being a dimensional measure of the microwave inhomogeneity. This modification has no impact on the computer resources it takes to solve the TDSE. For the simulations reported in this paper, the number of different \((h_{d,n},h_{i,n})\) pairs is equal to \(N_{\mathrm {S}}=28\).

Figure 9 shows the simulation data obtained by solving the TDSE with Hamiltonian Eq. (17). Although the microwave inhomogeneity is weak (\(\delta =0.1\)), it has considerable impact on the magnetization dynamics.

Fig. 8
figure 8

Simulation results for the total magnetization obtained by solving the TDSE with Hamiltonian Eq. (15) for a system of \(N_{\mathrm {s}}=28\) interacting spins, subject to a time-dependent magnetic field. The spin–spin coupling strength \(\lambda =5\), \(h_d=15\,\hbox {MHz}\), and \(h_i=0.12h_d=1.8\,\hbox {MHz}\)

Full size image
Fig. 9
figure 9

Simulation results obtained by solving the TDSE with Hamiltonian Eq. (17) for a system of \(N_{\mathrm {S}}=28\) interacting spins, spin-bath coupling \(\lambda =5\) and microwave field inhomogeneity \(\delta =0.1\)

Full size image

In the absence of spin–spin interactions (\(\lambda =0\)) and inhomogeneities, all the spins perform the same motion, as discussed in Sect. (4). If \(\lambda =0\), small differences in the microwave amplitudes \(h_{d,n}\) and \(h_{i,n}\), see Eq. (17), lead to quantitative but no qualitative differences in the motion of the spins. However, if \(\lambda \not =0\), the random interactions among the spins act against the coherent motion of the spins and the Rabi oscillations decay on a short time scale, as illustrated in Fig. 9a, f, except when the value of \(\Delta \) satisfies the condition for the first Floquet resonance, see Fig. 9c, h.

In the case of the bath-I model, see Sect. 5, the spin and the bath-I spins establish some kind of stationary state.. In the case of model Eq. (17), the interactions among all the spins may destroy the correlations between the individual spins as time proceeds. In the present case, we take as a measure for these correlations

$$\begin{aligned} C(t)=C_0\sum _{1=m<n}^{N_{\mathrm {S}}}\; \langle \Psi (t=0)\vert S_{m}^z(t) S_{n}^z(t)\vert \Psi (t=0) \rangle \;,\nonumber \\ \end{aligned}$$
(18)

where \(C_0={2}/{N_{\mathrm {S}}(N_{\mathrm {S}}-1)}\). Equation (18) is the sum of all equal-time correlations of the z-components of all spins. Here we do not present real-time data for C(t) (see Fig. S13)) but rather show its Fourier transform, see Fig. 10. As our main interest is in the long-time behavior of C(t), we discard the transient regime by considering real-time data for \(1\,\upmu \mathrm {s}\le t \le 5\,\upmu \mathrm {s}\) only.

In the absence of the image drive \(h_i=0\) (see Fig. 10a, d), the Fourier transformed C(t) exhibits almost no structure and has low intensity in comparison to the maximum intensities observed if the image drive is active (\(h_i=1.8\,\mathrm {MHz}\), see Fig. 10b, c, e, f).

If the image drive is present (\(h_i=1.8\,\mathrm {MHz}\), \(\Delta =5\,\mathrm {MHz}\) and \(\phi =0\)), Fig. 10b suggests a correlated motion with a frequency of \(\approx 10\,\mathrm {MHz} \approx 2 \Delta \) (and some weak higher harmonics). The spin system is outside the Floquet resonance regime and therefore the only source of correlation is imposed by the image field drive, similarly to the case of bath-I presented in the previous section. Since the image drive runs at \(2\Delta \), the correlation is driven at the same frequency.

This picture changes drastically if the image drive is present and \(\Delta \) satisfies the condition for the first Floquet resonance \(2\Delta =(h_d^2+\Delta )^{1/2}=F_\mathrm{R}\approx 17.3\,\mathrm {MHz}\). Then Fig. 10e suggests a correlated motion with a frequency of \(\approx 35\,\mathrm {MHz} \approx 4 \Delta =2F_\mathrm{R}\), that is a frequency of four (not two as in the case \(\Delta =5\,\mathrm {MHz}\)) times the Rabi frequency. Assuming that the z-component of each spin is oscillating with some frequency, then multiplying two such components (see \(S_{m}^z(t) S_{n}^z(t)\) in Eq. (18)) gives rise to an oscillation with twice that frequency. Thus, the clear signal at \(\approx 4 \Delta \) suggests that z-component of each spin is oscillating with a frequency of \(2\Delta \) imposed by the image drive and that there is some collective motion of these components. In other words, at the Floquet resonance, the image drive induces long-time correlations between the spins seen in the Fourier transform at a frequency of \(2F_\mathrm{R}\). For \(\phi =45^\circ \) the magnetization dynamics shows a much richer structure (see Fig. 9h) which is also reflected in a more complicated spectrum, see Fig. 10g, but the signal at \(2F_\mathrm{R}\) is also clearly visible.

Fig. 10
figure 10

Simulation results of the Fourier transformed correlation Eq. (18) obtained by solving the TDSE for Hamiltonian Eq. (17), describing a system of \(N_{\mathrm {S}}=28\) interacting spins, The microwave field inhomogeneity \(\delta =0.1\). The positions of \(2\Delta \) and \(4\Delta \) are indicated by the two vertical dashed lines. The plots b, c, e and f correspond to plots a, c, f, and h in Fig. 9, respectively. At the Floquet resonance, the spins tend to develop a synchronized motion seen as a signal at \(4\Delta \)

Full size image

The results the TDSE simulations for the model of interacting spins can be summarized as follows:

  1. 1.

    Qualitatively, model Eq. (15) does not reproduce the features of the Rabi oscillations shown in Fig. 1. In particular, there is a clear qualitative difference between Figs. 1h–j and 8h–j, that is if \(\phi =45^\circ \). Therefore, model Hamiltonian Eq. (15) alone does not seem to describe well such oscillations.

  2. 2.

    However, accounting for the inhomogeneity of the microwave field as in model Eq. (17) considerably improves the qualitative agreement with the experimental data (compare Fig. 1 with Fig. 9).

  3. 3.

    At the Floquet resonance, all interacting spins are developing a synchronised motion. The effect disappears quickly when the resonance condition is not met.

7 Conclusion

Following its experimental implementation presented in Ref. [26], a spin qubit subject to two different microwave drives, one that induces Rabi oscillations and a much weaker second one aiming to preserve coherence was shown to extend the coherence time of a spin qubit. In this paper, this protocol is analyzed theoretically by numerically solving the time-dependent Schrödinger equation of two different microscopic many-body systems. The models considered are: (i) a spin qubit interacting with a collection of two-level systems representing a bath (bath-I) and (ii) weakly interacting spin qubits in which the spin qubits themselves act as a bath (bath-II).

We discuss the conditions for properly describing the measured Rabi oscillations [26] and analyze the microscopic internal dynamics simulated for each type of bath. We find a small amount of saturation building up in the spin–spin correlations of bath-I at long time scales while the transitory regime shows signatures of the Floquet dynamics. For the second model with bath-II, one clearly observes the build up of a synchronised motion of the bath spins, sustained for very long times. Both microscopic models are shown to be capable of reproducing the main feature of the spin dynamics observed in the ESR experiments on CaWO\(_4\):Gd\(^{3+}\).

The results of our simulation work suggest that in order to capture the essence of the real-time dynamics of a spin as observed in the ESR experiments described above, the effect of dynamics of this spin on the dynamics of the bath degrees of freedom has to be taken into account. The dynamics of the bath degrees must be treated on the same footing as the dynamics of the spin itself. The build up of correlations between the spin and the bath spin result in a synchronous motion of all the spins, having the effect of partial preservation of the phases.

In the standard situation the bath determines the temperature. However, if the whole system is driven by the external field, as is case for the systems studied in this paper, it is not so clear what “temperature” means. Still, the effect of the energy exchange between the spin bath and the system is important, as our results also show. Under these circumstances the role of the bath is not to determine the temperature but to promote relaxation by energy exchange. Of course, if the Hamiltonian is time independent and the whole system is closed, its dynamics is governed by the TDSE and there is no dissipation. However, a spin that is subject to external, time-dependent microwave fields and in contact with an environment can, and also does, exchange energy with the spin bath. Thus, the total energy is not conserved and the exchange of energy between the spin and the spin bath is the source of dissipation and decoherence (in this case).”

Even if initially the bath energy is zero, the energy of the total system is not conserved and the final energy does not need to be zero because the whole system is driven. Nevertheless, the fact that also in this driven case, the energy of the single spin and the total energy approach stationary values is an important characteristic of the systems studied in the present paper.

We have demonstrated that with the computational resources that are available today, it is possible to simulate the real-time dynamics of many-body quantum spin systems over a time span that is experimentally relevant. Our numerical simulations show how correlations are building up in the baths when the Floquet protocol is implemented, leading to a preservation of the coherence of the qubit. Both microscopic models are shown to be capable of reproducing the main feature of the spin dynamics observed in the ESR experiments on CaWO\(_4\):Gd\(^{3+}\).

In view of the various physical mechanisms that may be at play, it would be best to consider a model that simultaneously includes several of them. Unfortunately, limitations imposed by computer resources force us to consider the effects of different mechanisms separately.