Semi-Empirical Haken-Strobl Model for Molecular Spin Qubits

Understanding the physical processes that determine the relaxation $T_{1}$ and dephasing $T_2$ times of molecular spin qubits is critical for envisioned applications in quantum metrology and information processing. Recent spin-echo $T_1$ measurements of solid-state molecular spin qubits have stimulated the development of quantum mechanical models for predicting intrinsic spin qubit timescales using first-principles electronic structure methods. We develop an alternative semi-empirical approach to construct Redfield quantum master equations for molecular spin qubits using a stochastic Haken-Strobl model for a central spin with a fluctuating gyromagnetic tensor due to spin-lattice interaction and a fluctuating local magnetic field due to interactions with other lattice spins. Using a vanadium-based spin qubit as a case study, we compute qubit population and decoherence timescales as a function of temperature and magnetic field using a bath spectral density parametrized with a small number of $T_{1}$ measurements. The theory quantitatively agrees with experimental data over a range of conditions beyond those used to parametrize the model, demonstrating the generalization potential of the method. The ability of the model to describe the temperature dependence of the ratio $T_2/T_1$ is discussed and possible applications for designing novel molecule-based quantum magnetometers are suggested.


Introduction
Electron spins in molecules have emerged as promising qubit systems [1][2][3][4][5] due to their relatively long intrinsic relaxation and decoherence timescales, comparable with traditional solid-state spin systems such as NV-centers [6,7].Molecular spin qubits can be assembled into highly stable crystalline structures with tunable qubit densities [8][9][10] and can be integrated into other solid-state platforms such as superconducting resonators for controlling two-qubit interactions [11][12][13][14][15][16].The vibrational and spin environment of molecular qubits can also be engineered using synthetic chemistry strategies [17][18][19].Future advances in spin qubit coherence and controllability would require a precise characterization of the mechanisms that lead to spin decoherence in molecular materials [20].
Multi-scale ab-initio modeling techniques have been proposed and successfully used to predict population (T 1 ) and decoherence (T 2 ) timescales of solid-state molecular spin qubits, only using the chemical composition and geometry of the underlying crystal lattice as input [21][22][23][24][25][26][27].First-principles techniques offer valuable atomistic insights that can be used for designing chemical strategies aimed at improving qubit performance.However, the large computational overhead of atomistic simulations is a challenge for the implementation of large-scale computational discovery strategies for molecular spin qubit materials.
Recently, the relaxation dynamics of high-spin single-molecule magnets was theoretically studied with a semi-empirical approach [28].The dynamics problem was reduced in complexity by partitioning the crystal Hamiltonian into a part that is straightforward to compute using quantum chemistry packages and a parametrized term that models the interaction of the central spin with lattice phonons.The free model parameters were then determined by comparing theoretical predictions with magnetic relaxation measurements.
In this work, we further develop this semi-empirical approach by constructing parametrized Redfield quantum master equations to describe the interaction of molecular spin qubits with lattice phonons and electron spin baths.The method is based on Haken-Strobl theory [29][30][31], which treats system-reservoir interactions as stochastic fluctuations of the system Hamiltonian.Haken-Strobl theory has been successfully used to study exciton transport and spectroscopy in molecular aggregates [30] and light-harvesting complexes [32].We use it here to model spin-lattice interaction as a random fluctuation of the molecular gyromagnetic tensor of the spin qubit, and interaction with a reservoir of electronic and nuclear spins embedded in the crystal lattice, which result in random fluctuations of the local magnetic field.
Semi-empirical Redfield tensors can be used to compute relaxation and dephasing times, T 1 and T 2 , for solid-state molecular spin qubits, by parametrizing the bath spectral density through a fitting procedure that matches the model prediction to a small number of T 1 measurements.We demonstrate the procedure using experimental data for two different vanadium (IV) spin qubits that have been recently developed [4,5].In both cases, the vanadium ion has an energetically isolated 3d 1  xy valence electronic configuration and thus represents an ideal S = 1/2 spin system (see figure 1).The theory gives quantitatively accurate predictions for T 1 far beyond the experimental conditions used to parametrize the model, thus demonstrating the potential for generality of the method.
The rest of the article is organized as follows: in section 2 we construct the proposed semi-empirical quantum Redfield equation.In section 3 we discuss the T 1 and T 2 results for two vanadium-based spin qubits as a function of temperature and magnetic field.In section 4 we conclude and discuss perspectives for future work.

Theoretical framework
We use Haken-Strobl theory [29,33] to describe the interaction of a central electron spin with magnetic and thermal environments.Following [34], the total spin Hamiltonian is written as Ĥ(t) = ĤS + ĤSB (t), where ĤS describes the stationary system and ĤSB (t) describes the fluctuations of the system due to system-bath interaction.To model the anisotropy of the electronic Zeeman interaction and capture the dependence of the spin dynamics on the magnetic field orientation [35], we write the total Hamiltonian as where µ B is the Bohr magneton, h the reduced Planck constant, ⃗ B the external magnetic field vector, ↔ g the gyromagnetic tensor, and ⃗ S the spin angular momentum of the molecular qubit.
We consider random fluctuations of the gyromagnetic tensor ← → δg (t) as well as local magnetic field fluctuations δB i (t).The latter result from the interaction of the central spin with other electronic and nuclear atomic spins the lattice structure.By neglecting higher order fluctuations, the total Hamiltonian can be partitioned into a central spin Zeeman term where g ij are gyromagnetic tensor elements and σj are Pauli matrices, and a system-bath Hamiltonian of the form where describes lattice-induced fluctuations of the g-tensor, and describe spin-bath fluctuations.Contraction over Cartesian indices i = (x, y, z) is implied throughout.We focus on the dynamics of the reduced density matrix ρs (t) = Tr B [ρ(t)] in the Born approximation for the total state ρ(t) = ρs (t) ⊗ ρB , where ρB is the stationary bath state.In the Markov approximation, the vibrational and spin degrees of freedom of the reservoir relax faster than the system timescales.In the interaction picture with respect to ĤS , the equation of motion for the matrix elements ρ s ab of the reduced spin density matrix in the basis of system eigenstates |a⟩, can generally be written as where R ab,cd are elements of the Redfield tensor which encodes the relaxation and decoherence dynamics of the central spin system.The tensor elements are given by [22] R ab,cd = 1 h2 with Γ ab,cd (ω dc ) being decay rate functions evaluated at the system transition frequencies ω ab = ω a − ω b .These rates can in turn be written in terms of spin selection rules and system-bath coupling strengths as where J jj ′ (ω) is the bath spectral density, which describes the coupling strength of the system with its environment at a given transition frequency and temperature [36].To model molecular spin qubits, we partition the total spectral density J jj ′ (ω) in equation ( 8) as here J δg jj ′ (ω) is the contribution from g-tensor fluctuations and J δB jj ′ (ω) describes field fluctuations due to couplings with the spin bath.
We construct J δg jj ′ by writing the bath autocorrelation function as with g-tensor fluctuations of the form The amplitude matrix A ij is assumed to be isotropic for simplicity, i.e.A ij = A g , α is a temperature scaling power, γ g and Ω g give the bandwidth and the resonance frequency of the spectral density.These four parameters can be calibrated using experimental measurements of T 1 , as explained below with an example.Equation (11) implies that the lattice-induced g-tensor fluctuations grows with temperature as a power law and have a finite correlation time τ g ≡ 1/γ g determined by the phonon spectrum, which is consistent with first-principles theory [22].The quantum mechanical spectral density J δg jj ′ (ω dc ) that defines the Redfield tensor is constructed by taking the Fourier transform of the bath autocorrelation in equation (10) multiplied by the semiclassical harmonic correction factor hω/2k B T, which is needed to satisfy detailed balance relations [37,38].We obtain ( This expression vanishes at ω = 0 and therefore does not contribute to pure dephasing [39].
We follow a similar procedure to construct an effective spectral density that describes coupling to the low-frequency spin-bath J δB jj ′ , introducing the bath autocorrelation function with the local magnetic field noise dynamics given by This expression describes pure dephasing with a finite amplitude at zero frequency given by where B is the magnetic field magnitude.The quadratic field dependence of A is compatible with a direct spin relaxation process [35].The dephasing correlation time is τ pd = 1/γ pd .The parameters a and c determine the magnetic field dependence of the bath fluctuations and can be determined from T 2 measurements.Taking the Fourier transform of equation ( 13), we obtain the spin bath spectral density To constrain the parameters {A g , α, γ g , Ω g , a, c} that determine the open quantum system model, we solve the Redfield equation for the reduced density matrix elements ρ s ab = ⟨a|ρ s |b⟩ in the energy eigenbasis, using the total spectral density J jj ′ (ω) from equation ( 9) with trial parameters.T 1 and T 2 times are obtained from the calculated decay dynamics of the populations ρ aa and coherences ρ ab , and the parameters are iteratively improved by comparing the theoretical T 1 and T 2 times with suitable experimental measurements.Once the spectral density parameters are fixed to match a small set of target measurements, the model predictions can be extrapolated over a broader set of conditions, as illustrated in what follows with two examples.

Results and discussion
We demonstrate the parametrization of the spectral densities using selected T 1 measurements as a function of temperature and magnetic field for VO(acac) 2 [qubit (1)] taken from [4] and [(Ph [5].For these molecular complexes, we used the electronic structure methods described in [40,41] to obtain the diagonal g-tensor elements g xx = 1.920 102, g yy = 1.978 044 and g zz = 1.981 320, for qubit (1), and g xx = 1.956 584, g yy = 1.988 596 and g zz = 1.989 734, for qubit (2).The calculated qubit spectra agree well with measured Zeeman splittings for B > 0.3 T. The near-zero field spectrum is not well reproduced because explicit hyperfine interactions are ignored in equation (1).
The spectral density parameters used to construct the Redfield tensor are chosen to match the reported value T exp 1 = 0.940 ms at B = 5 T and T = 10 K, for qubit (1) [4], and T exp 1 = 0.198 ms at B = 6.8 T and T = 5 K for qubit (2) [5].The parameters obtained for each model and qubit are given in table 1. Different model parameters give the same T exp 1 value because the parametrization procedure corresponds to an over-determined multi-valued optimization problem.Further constrains can be placed on the parameters using additional T 1 and T 2 measurements.For qubit (1) we obtain Ω g = 11 cm −1 , γ g = 2.387 cm −1 , and γ pd = 0.1 cm −1 .The temperature scaling power α = 4.6 was extracted from the measured temperature dependence of the T 1 time at B = 5 T [4].At the magnetic fields and temperatures used, the fitting procedure Table 1.Spectral density parameters used in the T1 and T2 calculations in the text.The amplitude of the fluctuations of the local magnetic field are given by the polynomial A(B) = a + cB 2 .These sets of parameters are chosen to match the experimental measurement T exp 1 = 0.940 ms from [4], for qubit (1) at B = 5 T and T = 10 K, and T exp 1 = 0.198 ms from [5] for qubit (2) at B = 6.8 T and T = 5 K.In all models the temperature scaling parameter is α = 4.6 for (1) and α = 2.2 for (2).
In table 1, we specify three sets of model parameters that give the same T exp 1 .Model I assumes that only spin-lattice relaxation occurs.Models II assumes that the magnetic field fluctuations do not depend on the magnetic field magnitude and Model III assumes a quadratic dependence with B. The magnitude of the bath-induced Hamiltonian fluctuations are defined as δg ≡ √ ⟨δg 2 ii (0)⟩ and δB ≡ √ ⟨δB 2 i (0)⟩.

Vanadyl qubit relaxation
Figure 2 shows the spectral density J xx (ω) for qubit (1) corresponding to Model III parameters in table 1.The individual contributions of J δg xx and J δB xx are also shown.At higher frequencies (large magnetic fields), J δg jj ′ dominates the system-bath interaction by orders of magnitude, but its contribution to pure dephasing below 0.3 T is negligibly small.The finite zero-frequency amplitude of the spin bath spectral density ensures that the molecular spin qubit is coupled to the reservoir over the entire frequency range.Figure 3(a) shows the calculated T 1 times for qubit (1) as a function of magnetic field at T = 10 K.The predictions of different model parametrizations are compared to experimental data from [4].As expected and in agreement with ab-initio calculations [22], the reservoir model without magnetic field fluctuations (Model I) is in good agreement with T exp 1 at higher magnetic fields, where direct spin-phonon coupling dominates the relaxation dynamics.However, Model I is unable to capture the crossover around B ≈ 1 T, below which T 1 decreases with decreasing magnetic field, as the influence of the low-frequency spin reservoir is stronger.
By including magnetic field fluctuations (Models II and III), the theory continues to reproduce the T exp 1 values at higher magnetic fields, but now also capture the reduction of T 1 below the crossover.The spin-bath model gives T 1 values of the same order as the low-field measurements (B < 0.3 T), but quantitative agreement is limited.For T = 10 K, we find that using a field-dependent magnetic noise model (Model III) does not significantly modify the calculated crossover position and the T 1 values relative a theory with field-independent spin noise (Model II). Figure 3(b) shows the calculated T 1 times for qubit (2) as a function of magnetic field at T = 5 K. Theoretical predictions are compared with measurements from [5].All models parametrizations give T 1 values in agreement with experiments at higher fields (B > 3 T) but Models I and II deviate significantly from the experimental trend at magnetic fields below the crossover.

Effective spin dephasing models
The Redfield equation in the eigenbasis of ĤS describes population transfer and decoherence of spin sublevels.The tensor in equation ( 7) also contains off-diagonal non-secular contributions that couple populations and coherences, depending on the details of the spectral density [42][43][44].For secular open quantum systems, populations and coherences evolve independently [45].In that case, dephasing and relaxation times are related by 1/T 2 = 1/2T 1 + 1/T * 2 , where T * 2 is the contribution from pure dephasing.Since pure dephasing does not involve transitions between system eigenstates, it is only present when the spectral density has finite amplitude at zero frequency.Therefore, only J δB contributes to pure dephasing.The phonon bath spectral density scales as J δg ∼ ω in the low frequency limit, and thus can only contribute to T 2 via spin relaxation (e.g.1/2T 1 ) and non-secular processes.
For systems where only T 1 measurements are available to parametrize the total spectral density, the spin bath model parameters a, c, and γ pd are not fully constrained.There can be magnetic fields and temperatures for which different bath models give the same prediction for T 1 , but differ in the estimation of T 2 .This is clear in table 1, where three sets of model parameters match the same T exp 1 , but predict T 2 times that can differ by up to an order magnitude.
Figure 4(a) shows predictions for the dephasing-to-relaxation ratio T 2 /T 1 as a function of magnetic field for qubit (1).The three sets of model bath parameters from table 1 are compared.As expected, when the spin bath is neglected (Model I), the pure relaxation limit T 2 = 2T 1 is obtained for a broad magnetic field interval.Non-secular terms in the Redfield tensor increase the ratio beyond the pure-relaxation limit at higher fields (B > 3 T).By including the contribution to dephasing from the spin bath spectral density (Models II and III), T 2 drops below T 1 at low fields, with a zero-field ratio T 2 /T 1 given √ A(0) (see equation (15)).The ratio further decreases with increasing magnetic field, mostly due to the increase of T 1 up to the crossover field B ≈ 1 T.For magnetic fields beyond the crossover, the qubit spin dynamics becomes increasingly dominated by phonon-induced relaxation, making T 2 /T 1 tend towards the pure relaxation limit at high fields.The spin bath model with field-dependent fluctuation amplitude (Model III) has a slower approach to the pure relaxation limit, as in this case T 2 ∼ 1/B 2 .
Figure 4(b) shows T 2 /T 1 as a function of magnetic field for qubit (2).The results are qualitatively similar to qubit (1).Model I captures the pure relaxation limit over a broad range of magnetic fields with non-secular deviations occurring from B ≈ 1 T.For Model III, we have T 2 < T 1 throughout the interval, because the local field fluctuations δB are relatively strong in this case, even at high magnetic fields.
Figure 5(a) shows the temperature dependence of the ratio T 2 /T 1 for qubit (1), as predicted by Model III.For low magnetic fields (B = 0.15 T), the ratio is largely insensitive to temperature, with only a modest increase above T ∼ 100 K due to the decrease of T 1 .In this regime, the system dynamics is mainly governed by local magnetic field noise, which is temperature independent in the model.At intermediate fields = 1 T), the dissipative dynamics has comparable contributions from magnetic field noise and g-tensor fluctuations, which scale as δg ∼ T α/2 with temperature.Since δg dominates at higher fields (B = 5 T), the temperature dependence becomes stronger.The T → 0 limit of Model III is strongly field dependent, as the relative contributions of the fluctuations δB and δg can vary across different field regimes (see figure 2). Figure 5(b) shows the corresponding temperature dependence of the ratio T 2 /T 1 predicted by Model III for qubit (2).The behaviour at different magnetic fields is qualitatively similar to qubit (1), with the difference that T 2 < T 1 over the entire temperature range.
The temperature dependence of T 2 is directly related with the relative strength of the qubit coupling to the phonon and spin baths, characterized by the spectral densities J δg and J δB , respectively.The model only gives T 2 < T 1 from cryogenic to room temperature when magnetic field noise is not negligible.In our model this implies that the magnetic field should remain below or near the T 1 crossover in figure 3. Another possibility would be to introduce a temperature dependence to the zero-field fluctuation parameter, i.e. a → a 0 T β , with a scaling parameter β to be determined from experiments or ab-initio simulations.

Conclusions and outlook
We introduced a semi-empirical Haken-Strobl method to construct Redfield quantum master equations that describe the population and coherence dynamics of solid-state molecular spin qubits in an external static magnetic field, as a function of temperature.In Haken-Strobl theory, the interaction of the molecular spin with lattice phonons is treated as a stochastic fluctuation of the gyromagnetic tensor and the interaction with other atomic spins in the crystal as fluctuations of the local magnetic field experienced by the qubit.Analytical expressions for the associated bath spectral autocorrelation functions can be parametrized using a small set of experimental measurements.We demonstrate the parametrization procedure using experimental T 1 data for two recently reported vanadium (IV) spin qubit implementations [4,5].Quantitative agreement is obtained with respect to spin relaxation measurements at high magnetic fields, where spin-lattice relaxation dominates.Qualitative agreement is found at low magnetic fields, where spin-spin interactions become dominant.The model captures the experimentally observed optimal T 1 time as a function of magnetic field strength, which marks the crossover from a regime dominated by the spin-bath interaction at low magnetic fields, from the regime where qubit relaxation is dominated by spin-lattice interaction.
The generalization ability of the proposed semi-empirical approach beyond the set of conditions used to parametrize the Redfield tensor could improve significantly by comparing the model prediction with simultaneous measurements of the dephasing and relaxation times, T 2 and T 1 , over a broad range of temperatures T ∼ 1 − 10 K and magnetic fields B ∼ 0.1 − 10 T. Complementary ab-initio calculations may also be carried out to reduce the number of free parameters and qubit measurements needed to specify the relaxation model, thus facilitating the analysis of hypothetical crystal structures of potential interest for novel qubit implementations.Explicitly modeling the hyperfine interaction of the central molecular spin with the local spin-bath would enable the treatment of spin clock transitions [46] and a more accurate analysis of the system dynamics at low magnetic fields.Such model improvements would enable a better understanding of the factors that limit the spin dephasing time T 2 at room temperature, a key feature to optimize in the development of quantum magnetometers based in solid-state molecular spins [19].

Figure
Figure (a) Spin relaxation time T1 as a function of magnetic field B for qubit (1).Predictions using three spectral density models parametrized as in table 1 are compared with experimental data at T = 10 K from[4].(b) Same as panel (a) for qubit (2), with experimental data at T = 5 K from[5].

Figure 4 .
Figure 4. (a) Dephasing-to-relaxation ratio T2/T1 as a function of magnetic field for qubit (1) at T = 10 K. Model parameters are given in table 1.The dashed horizontal line shows the pure relaxation limit T2 = 2T1.(b) Same as panel (a) for qubit (2) models at T = 5 K.

Figure 5 .
Figure 5. (a) Dephasing-to-relaxation ratio T2/T1 as a function of temperature for qubit (1) at low (solid line), intermediate (dashed line) and high magnetic fields (dot-dashed line).Model III predictions are shown, with parameters given in table 1.The dashed horizontal line shows the pure relaxation limit T2 = 2T1.(b) Same as panel (a) for qubit (2).