Electronic spin-spin decoherence contribution in molecular qubits by quantum unitary spin dynamics

TL;DR

This paper introduces a computational method to analyze spin decoherence in molecular qubits, highlighting the impact of long-range interactions on coherence times and spin-echo signals.

Contribution

It presents a novel quantum unitary spin dynamics approach to understand and quantify spin decoherence mechanisms in molecular qubits.

Findings

Long-range electronic spin-spin interactions significantly affect decoherence.

Magnetic dilution influences spin dephasing behavior.

The methodology enables detailed analysis of dephasing pathways.

Abstract

The realisation of quantum computers based on molecular electronic spins requires the design of qubits with very long coherence times, T2. Dephasing can proceed over several different microscopic pathways, active at the same time and in different regimes. This makes the rationalisation of the dephasing process not straightforward. Here we present a computational methodology able to address spin decoherence processes for a general ensemble of spins. The method consists in the propagation of the unitary quantum spin dynamics on a reduced Hilbert space. Then we study the dependence of spin dephasing over the magnetic dilution for a crystal of Vanadyl-based molecular qubits. Our results show the importance of long-range electronic spin-spin interactions and their effect on the shape of the spin-echo signal.

Figures

Equations

$${\hat{H}}_{\mathrm{S}}=\sum _i^{N_{\mathrm{S}}}{\mu }_{\mathrm{B}}\vec{\mathbf{B}}\cdot \mathbf{g}(i)\cdot \vec{\mathbf{S}}(i)+\frac{1}{2}\sum _{ij}^{N_{\mathrm{S}}}\vec{\mathbf{S}}(i)\cdot {\mathbf{D}}^{\mathrm{dip}}(ij)\cdot \vec{\mathbf{S}}(j),$$

$$\frac{d\hat{\rho }}{dt}=-\frac{i}{\hslash }[{\hat{H}}_s,\hat{\rho }].$$

$$\{|S_z(1),S_z(2),...,S_z(N)⟩\}=\underset{i=1}{\overset{N}{⨂}}\{|S_z(i)⟩\}.$$

$${\rho }_{ab}(t)=U_{ab,ab}(t-t_0){\rho }_{ab}(t_0)=e^{-i{\omega }_{ab}(t-t_0)}{\rho }_{ab}(t_0),$$

$$\rho (0)\stackrel{R(\frac{\pi }{2})}{\to }{\rho }^{\frac{\pi }{2}}(0)\stackrel{U(\tau )}{\to }{\rho }^{\frac{\pi }{2}}(\tau )\stackrel{R(\pi )}{\to }{\rho }^{\frac{3\pi }{2}}(\tau )\stackrel{U(\tau )}{\to }{\rho }^{\frac{3\pi }{2}}(2\tau ),$$

$$M_x(t)=(M_x(0)-M_x(\infty ))e^{-(\frac{t}{T_2}{)}^2}+M_x(\infty ).$$

$$M_x(t)=(M_x(0)-M_x(\infty ))e^{-(\frac{t}{T_2}{)}^{\beta }}+M_x(\infty ).$$

Full text

I Introduction

In recent years, the possibility of using quantum properties of materials as active elements of new technologies has emerged, setting the foundation for a quantum revolution. Quantum computing, even though theorised several years ago, has only recently been implemented in several physical systems, ranging from Cooper’s pairs in superconductors Nakamura et al. (1999) to electronic spins Godfrin et al. (2017). In a nut-shell, the vast dimension of a quantum system’s Hilbert space and the quantum mechanical way to operate over it, allows one to perform operations at a speed much faster than in classical computers. The basic working principle of this technology requires the ability to initialise, manipulate and detect a many-body quantum state in its Hilbert space DiVincenzo (2000). The coherence over time of the information “stored” in the quantum Hilbert space, generally referred as the characteristic time, $T_{2}$, is an essential property of a given system and depends on both intrinsic Atzori et al. (2016a) and extrinsic factors Zadrozny et al. (2015). Quantum coherence has to be preserved long enough for the computation to be carried out and it represents the basic requirement for any class of potential qubits.

Among the most promising paradigms to build quantum computers, transition-metal-based magnetic molecules have been brought forward as potential qubits, thanks to their long spin-coherence times and easy tunabilityLeuenberger and Loss (2001); Ardavan et al. (2007); Winpenny (2008); Lehmann et al. (2009); Wedge et al. (2012); Ferrando-Soria et al. (2016); Shiddiq et al. (2016); Gaita-Ariño et al. (2019). Indeed, molecular compounds represent a rather rich materials platform, where several synthetic strategies can be exploited to selectively tune specific inter- or intra-molecular interactions, with the goal of optimising the compound and extracting the basic physical picture of the dephasing. In recent years, an extensive amount of work has been performed in order to extend both the spin-lattice and spin-coherence times. Several guidelines to build more robust molecular qubits had been proposed, highlighting the importance of both the nature of the first coordination shell Atzori et al. (2016a) and of surrounding nuclear spins density Zadrozny et al. (2015).

From a rational design point of view, the optimisation process has to be led by a detailed knowledge of those microscopic mechanisms leading to spin dephasing and this represents the focus of this work. Restricting our study to the sole spin-spin interactions, i.e. assuming the temperature to be low enough to exclude spin-lattice relaxation processes, several sources of spin dephasing might be operative. In a very general fashion, every spin in the system, other than the one addressed, cause a dynamical fluctuation of the spin local magnetic field, generating decoherence. In particular, three different types of spin-spin interactions are generally present in a crystal of molecular magnets: electronic-electronic dipolar interactions, nuclei hyperfine interaction and electronic-nuclear dipolar interactions. From a computational point of view, several approaches has been developed to study the central spin problem Yao et al. (2006); Witzel and Das Sarma (2006); Yang and Liu (2008); Maze et al. (2008), where a single electronic spin loses its correlation due to the interaction with a bath of nuclear spins, finding application in both solid-state defect quantum bits Seo et al. (2016) and molecular qubits Bader et al. (2016). In order to complement the existing picture of spin dephasing, we present here a computational strategy to address the role of electronic spin-spin interactions over the dephasing time for an ensemble of generic spins. This approach, able to address up to hundreds of spins at the same time, is used to simulate the spin-echo decay experiment for a crystal of Vanadyl qubits, revealing a linear dependence of $T_{2}$ with respect to the magnetic dilution. Moreover, insights regarding the shape of the spin-echo envelop are reported, showing a transition from a Gaussian profile, for a perfectly ordered crystal, to a stretched exponential, for the random-distributed magnetically diluted phase.

II Unitary Quantum Dynamics and Spin Echo

A spin initialised in a specific state is always subject to dephasing, unless it represents a perfect close system, namely it is isolated from the rest of the universe. Such a strict condition is never met in practice, so that a quantum system is never fully isolated and thus is subject to some non-trivial dynamics. The study of this process, in principles, requires the knowledge of all the spin-“external world” interactions and the solution of the Schrödinger equation for the coupled system.

In order to formalise the problem, let us introduce the total spin Hamiltonian for a system of $N_{\mathrm{S}}$ coupled spin 1/2, $\mathbf{S}(i)$,

[TABLE]

where the interaction of the spins with an external field, $\vec{\mathbf{B}}$, is gauged by the Bohr magneton, $\mu_{\mathrm{B}}$, and the $g$-tensor, $\mathbf{g}(i)$, while the spin-spin interaction, $\mathbf{D}^{\mathrm{dip}}$, is assumed to be only dipolar in nature. The state of the spin system can be described in term of the spin density matrix, $\rho$, whose dynamics is regulated by the Liouville equation

[TABLE]

The matrix elements of $\hat{H}_{\mathrm{S}}$ are usually computed on a basis set corresponding to the tensor product of all single spins $\hat{S}_{z}(i)$ operators’ eigenkets

[TABLE]

By diagonalising $H_{\mathrm{S}}$ one can solve numerically Eq. (2) to obtain an expression for the time evolution of $\rho_{ab}$ in terms of the eigenvalues of $H_{\mathrm{S}}$, $E_{a}$,

[TABLE]

where $\omega_{ab}=(E_{a}-E_{b})/\hbar$ and $U_{ab,ab}(t-t_{0})$ is the unitary evolution operator.

Eq. 4 describes the exact dynamics of the $N_{\mathrm{S}}$ interacting spins, but its use is limited by the size of the matrix that one can actually diagonalise. The complexity of the problem arises mainly from the exponential dependence of the size of the spin Hilbert space on $N_{\mathrm{S}}$. This limits the solution to $\sim 16$ spins on common parallel machines, a general problem in quantum many-body physics. Nevertheless, several approaches to overcome such limitation exist. The specific problem concerning the correlation loss of a central spin embedded in a spins bath has been successfully tackled with a range of methods falling under the name of cluster approximations Yao et al. (2006); Witzel and Das Sarma (2006); Yang and Liu (2008); Maze et al. (2008), where the total contribution to dephasing is decomposed into the contribution of small spin clusters. Another possible approach, somehow based on a similar physical principle, consists in working on a truncated Liouvillian space, where the states considered are chosen among the ones deemed relevant. The basic idea behind these methods consists in assuming that only a small part of the total Liouville space will be actually explored during the time evolution. Moreover, the states explored will be those close enough to the states at time $t=0$. Coupled clusters Roger and Hetherington (1990), density-matrix renormalisation group Daley et al. (2004) and adaptive state-space restriction Kuprov et al. (2007) approaches all fall into this category.

In this paper we will take advantage of a non-perturbative Hilbert space size reduction technique, belonging to the last class of methods and inspired by Kuprov’s works Kuprov et al. (2007). In the presence of an external magnetic field pointing along the $z$ direction and in the limit $B_{z}\rightarrow\infty$, the ground state will tend to a state with all the spins aligned towards the field. In this situation, the ground state can be represented with only one vector in the $S_{z}$ basis set and all the excited states are represented by single, double, …, $N$-th fold spin flips. In the presence of spin-spin interactions a mixing between these eigenstates is introduced, generating a non-trivial dynamics, leading to decoherence. However, as long as the Zeeman term remains the leading interaction, the nature of the eigenstates will be only slightly affected. Therefore, by assuming that the starting configuration can be generated close to the ground state of the system, we can restrict the size of the Hilbert space to those states within a few spin flips from the initial state. Such a restriction makes the total number of spin states needed to investigate the dynamics to scale polynomially with $N_{\mathrm{S}}$ and the problem expressed by Eq. (2) and Eq. (4) becomes easily tractable even for $10^{2}$ spins.

The method just outlined is here employed to study the effect of electronic spin-spin interactions in the decoherence processes through the simulation of the spin Echo experiment. The spin Echo experiment is often use as a mesure of T2, as it averages out inhomogeneous effects, and it consists in a two pulses sequence. The first pulse brings the magnetization along the x-axis, then, after a free evolution time $\tau$, a $\pi$ pulse is applied. After the second pulse, another interval $\tau$ of free evolution time is waited before reading out the x component of the residual magnetization. This can be described by the path,

[TABLE]

where $\hat{\rho}(0)$ is the starting density matrix, $U(t)$ is the time propagator operator as expressed in Eq. (4) and $R(\theta)$ is a rotation operator that simulates the effect of a micro-wave electron paramagnetic resonance (EPR) pulse.

III Computational Details

The implementation of the method just outlined follows three steps: 1) the construction of the total Hamiltonian matrix, 2) the diagonalization and 3) the propagation of an initial density matrix according to the scheme introduced in Eq. (5). The truncated basis set is generated by flipping a given number of spins with respect to a state where they are all aligned along the field direction. When possible, the spin Hamiltonian is rotated in a block-diagonal form by forming a translation-invariant basis set. The spin Hamiltonian diagonalisation is performed with scalapack libraries Choi et al. (1996). The initial spin-density matrix is generated as a pure state where all the spins are aligned along the external magnetic field. Tests have also been performed with an initial spin-density matrix equal to the canonical ensemble distribution, but no differences in the extracted T2 values have been observed. The rotation operator is constructed by converging a Taylor expansion of $R(\theta)=\mathrm{exp}[i(\vec{\mathbf{S}}\cdot\vec{\mathbf{n}})\theta/2\pi]$, where $\vec{\mathbf{S}}$ is the spin-vector operator of the spin to be rotated, $\vec{\mathbf{n}}$ is the normal rotation axis and $\theta$ is the angle of rotation. The time evolution of the expectation value for a general operator $\hat{A}$ is $\langle A(t)\rangle=\mathrm{Tr}\{\hat{A}\hat{\rho}(t)\}$.

We have also tested an alternative to the diagonalisation of the spin Hamiltonian matrix that involves the use of the original $S_{z}$ basis set. In this basis set almost all the operators have a very sparse matrix representation that can be conveniently exploited for storing purposes. However, the time propagator $\hat{U}(t-t_{0})$, being an exponential operator, has a dense matrix representation and must be multiplied to the density matrix twice at each time-step. Moreover, the generation of the exponential operator in this basis set requires the use of a very short time step $\Delta t=t-t_{0}$, usually requiring the calculation of more than 106 steps to reach the desired microsecond time-scale. Consequently, here we employ the direct diagonalisation approach that makes it possible to set the desired time step in the dynamics propagation. The $\mathbf{g}$ operators, when not available from experiments, are computed by CASSCF calculations with the software Orca Neese (2012). An active space including the $d$ electrons and five 3$d$ orbitals is used. The spin-spin distance used to compute the dipole-dipole interaction is implemented using periodic boundary conditions in order to correctly reproduce the solid state environment.

IV Results

In this work we study two molecular systems belonging to the Vanadyl class of molecular quibits: VO(acac)2 Tesi et al. (2016), being acac=acetylacetonate, and VO(dmit)2 Atzori et al. (2016a), with dmit = 1,3-dithiole-2-thione-4,5-dithiolate. These magnetic molecules have recently become subjects of intense research as they naturally show rather long coherence times, up to several $\mu s$ even at room temperature Atzori et al. (2016b), suggesting the possibility of further performance increase if specifically tailored.

The first system investigated is the VO(acac)2 crystal. The phase-coherence properties of this compound have not been investigated experimentally, and here we use it to show the consistency of our computational method. The $\mathbf{g}$ tensor for this specie is calculated as explained in the “Computational Details” section and it reproduces the typical values for this class of molecules: $g_{xx}=1.9430$, $g_{xy}=0.0179$, $g_{xz}=-0.0170$, $g_{yy}=1.9749$, $g_{yz}=0.0084$ and $g_{zz}=1.9735$. The orientation of the molecule in space is as found in the crystal’s unit-cell aligned with the crystallographic $\vec{\mathbf{a}}$ vector along $x$ and the $\vec{\mathbf{b}}$ vector lying in the $xy$ plane. The lattice parameters can be found in ref. Tesi et al. (2016). Finally the external field was chosen aligned along the $z$ direction. Fig. 1 reports the time evolution of the expectation value of the transverse magnetisation, $M_{x}$, when the multiplicity of the spin excitations included in the basis set of Eq. (3) is changed. The study is performed on a 2$\times$2$\times$2 super-cell of the VO(acac)2 crystal, containing a total of 16 independent and interacting spins. From the figure it is possible to see how $\langle M_{x}\rangle$ for a given spin initialised at its maximum value after the $\pi/2$ pulse, decays with a Gaussian profile

[TABLE]

The calculated time constant is in the nanosecond range. The short-time spin dynamics converges very rapidly with respect to the inclusion of multiple spin-flips in the basis set. In fact, already for the minimal basis set, including only double spin-flips, the decay part of the dynamics is well converged and only long-time features need to be perfected. This is an important aspect as it allows us to use the minimal basis-set to study the dephasing process, being the short-time dynamics the important one.

Strong of these results we use the two spin-flips basis-set to study the dependence of the $\langle M_{x}\rangle$ dynamics on the total number of spins included in the simulation. This is reported in Fig. 2. The spin dynamics features converge fast with respect to the number of interacting spins, and virtually no difference is observed on the short-time profile for 3$\times$3$\times$3 and 4$\times$4$\times$4 VO(acac)2 super-cells, containing 54 and 128 spins, respectively. It is also interesting to note how the long-time fluctuations get damped with the inclusion of more spins in the simulation. This happens as a result of an increased size of the Hilbert space, which makes the chance of the system to visit states with a large magnetisation less probable, in analogy to the classical recurrent time in the phase space.

Being the limits of our method understood, we are now ready to study the effects of the magnetic dilution on the transverse magnetisation decaying profile. A simple model for addressing this problem consists in designing a lattice of spins with the same symmetry of the crystal unit-cell and rescale the size of the lattice vectors in order to modulate the spin-spin distance as in a dilution experiment. By applying this methodology to 27 VO(acac)2 spins, we scan magnetic dilution values ranging from 50% to 0.78% and the results are reported in Fig. (3). The dilution percentage is calculated as the number of spins with respect the total number of Vanadium sites availble in the volume considered. The echo signal presents the same features regardless of the dilution, since the symmetry of the interactions is preserved. However, the time-scale changes, leading to larger $T_{2}$ values for diluted systems. These reach values of 0.1 $\mu s$ for $\sim 1\%$ dilution.

Fig. 4 shows the $T_{2}$, extracted by fitting the echo profiles to Eq. (6), as a function of the dilution, revealing a linear dependence. This fact has also been previously noted in the context of nuclear magnetic impurities interacting with an electronic spin Maze et al. (2008) and in diluted radicals Brown (1973). It originates from the fact that, even though the dipolar interactions goes as $R^{-3}$, the average distance between spins does not scale linearly with the dilution percentage.

Interestingly, even though the predicted order of magnitude of the spin dephasing is approaching the correct value for the typical dilutions used in experiments, the shape of the spin-echo envelope is qualitatively different. Spin-echo signals are usually interpreted with generalised exponentials as

[TABLE]

For this class of molecular complexes an experimental value $\beta<1$ is found Atzori et al. (2016b), while here the profile is Gaussian and $\beta=2$.

In order to investigate this discrepancy we turn our attention to the study of another molecule in the Vanadyl family, namely VO(dmit)2. This compound has been extensively investigated experimentally Atzori et al. (2016a) and the spin-echo decay profile is available. Furthermore, an experimental $\mathbf{g}$ tensor has been extracted and it is used here: $g_{xx}=g_{yy}=1.9843$ and $g_{zz}=1.9735$. In this second example we overcome the limitations of the previous diluted study by designing several super-cells, where the Vanadium sites are randomly occupied by a number of spins consistent with the desired level of dilution. A first set of 80 supercells contains 48 spins located on the sites of a 6$\times$8$\times$6 VO(dmit)2 supercell, and a second set of 80 contains 70 spins located on the sites of a 7$\times$9$\times$7 supercell. All the supercells correspond to a dilution of 4%, a value comparable with experimental conditions. The spin echoes for all the different simulations performed on the same number of spins are then averaged. This protocol offers a more realistic representation of the experimental conditions, where spins experience different environments.

Our results are summarised in Fig. 5, where we show a non-Gaussian decaying profile of $M_{x}$ for both the sets of supercells, with the only noticeable difference among them being the value of $M_{x}(\infty)$. The fit to Eq. (7) of the curve corresponding to the simulation with 70 spins returns $\beta=0.70$ and $T_{2}=0.20~{}\mu s$. In comparison, VO(acac)2 relaxes with $T_{2}\sim 0.10~{}\mu s$ for the same level of dilution, as extrapolated from Fig.4. Notably, the stretch factor $\beta$ is now in excellent agreement with the experimental phenomenology. At the same time we predict a $T_{2}$ slightly shorter than the experimental one, which is estimated to be $\sim 5\mu s$. This is an interesting result, since usually an incomplete model overestimates the relaxation time, while our predictions surpass the experimental rate of decoherence. This suggests that there is an active mechanism, not cover in our simulations, slowing down the dephasing process. However, other sources of error might may also come from the simulated pulses being instantaneous and from the fact that in experiments multiple spins resonate at the same frequency and would flip at once under micro-wave pulses. Our model instead assumes that one single electronic spin is affected by the pulse.

V Discussion

The study of spin dephasing in molecular magnets has been mainly driven by experimental investigation and over the years several relevant figures of merit for the optimisation of $T_{2}$ have been proposed. However, the main focus of these investigations concern the role of the H1 spins in the solid-state matrix. The physics of this interaction is relatively well understood and the basic guideline to enhance $T_{2}$ consists in synthesising organic ligands and crystal’s matrices without elements with abundant spin-active nuclear isotopes. The role of the residual electronic spins interactions is usually less investigated and dilution is generally thought to be sufficient to eliminate it.

Here, however, we have shown that $T_{2}$ is only linearly dependent on the dilution and extremely low concentrations are needed to bring it beyond the $\mu s$ boundary. According to Fig. 3, a 0.01% dilution (or 0.3 mM) is needed to obtain $T_{2}>10\mu s$. This estimate also sets the minimal distance between qubits inside a device to 15 nm. For levels of dilution lower than this threshold, the electronic spin-spin contribution is comparable to that of H1. Therefore, the method proposed here, is expected to be of help in the discrimination processes originating from different interactions. A less investigated aspect, due to the synthetic complexity in gauging such interaction, concerns the role of the hyperfine coupling of the nuclear spin borne by the same atom hosting the electronic one. The method presented in this work can also be applied to this situation, as well as to the central electronic spin interacting with an H1 spin bath, and it will be argument of future applications. Another important ingredient that can be added to our model, is the inclusion of temperature effects and spin-phonon relaxation. As we approach ambient temperature, spin-phonon coupling becomes the bottle-neck for the coherence time Bader et al. (2014); Atzori et al. (2016a). In particular the understanding of the correlation between the chemical structure and the spin-lattice relaxation might lead to new, unforeseen synthetic designing rules Atzori et al. (2017); Lunghi et al. (2017a, b) and it must be considered a research priority.

VI Conclusions

In conclusion, we have here presented a general methodology to generate a truncated Hilbert space, designed to reproduce the dynamics of interacting spins in external magnetic fields. The method has been tested for VO(acac)2 and VO(dmit)2 molecular qubits, showing that even simple two-spin flips and second neighbour shells are able to capture the relevant physics leading to spin dephasing. The study of VO(acac)2 has also been extended to reproduce magnetic dilution experiments and has shown a linear dependence of $T_{2}$ with respect to the dilution. Our results therefore highlight the important role played by the residual electronic spin dipolar interaction in crystals of molecular qubits in setting a upper limit to $T_{2}$, reachable at not-extreme dilutions. Finally, our simulations also shed some light on the origin of the stretched exponential-type of decay usually encountered in EPR investigations of transition metal complexes, pointing to the disorder of magnetic centres inside the solid-state host as its origin.

**Data Availability

**All the relevant data discussed in the present paper are available from the authors upon request.

**Acknowledgments

**This work has been sponsored by Science Foundation Ireland (grant 14/IA/2624). We acknowledge the MOLSPIN COST action CA15128. Computational resources were provided by the Trinity Centre for High Performance Computing (TCHPC) and the Irish Centre for High-End Computing (ICHEC).

**Author contributions

**All the authors discussed the results and contributed to the manuscript.

**Competing financial interests

**The authors declare no competing financial interests.