Hyperfine Interaction in a MoS₂ Quantum Dot: Decoherence of a Spin-Valley Qubit

Modern life is powered by information, which leads to an increase in the demand for computational power so much that new generations of semiconductor technologies are employed consistently. Since then, there have been some reports headed for the development and improvement of valleytronic devices [1, 2], and TMD qubits, including valley qubits, spin qubits, spin-valley qubits and even impurity based qubits [3, 4, 5, 6, 7, 8, 9]. Quantum dots in a monolayer transition metal dichalcogenide (TMD-QD) such as Molybdenum Disulfide MoS₂ is a new pattern for qubit and holds the promises for new quantum devices but it still suffer from decoherence. Precisely, the major obstacles that are addressed pose a challenge to researchers, noise hinders the transmission of quantum signals. One of the essential tools in the study of spin-valley qubits in quantum dots is stability, which is delicate and very sensitive to the disturbances of a noisy environment. Indeed, the inevitable coupling and that each tiny interaction of these objects with their environment very quickly destroys the phase relations (superposition of incompatible states) between the quantum states until they become classical states. The electron spin in a quantum dot has two main decoherence channels, a (Markovian) phonon-assisted relaxation channel, due to the presence of spin-orbit interaction, and a (non-Markovian) spin bath constituted by the spins of the nuclei in the quantum dot that interact with the electron spin via the hyperfine interaction [10, 11, 12, 13] where the number of nuclear spins ranges from $\sim$ $10^{2}$ up to $10^{6}$ and full polarized baths are employed to facilitate qubit operations and extend coherence times.
A powerful tool for dealing with such systems is provided by an alternative way of deriving an exact master equation is the second-order time-convolutionless (TCL) master equation [14, 15, 16] of initial bath states, including the full polarized bath. Although we will focus on the example of spin qubits, our results are potentially applicable to any central spin problem.
In this work, we solve the central spin problem using the TCL master equation. This equation enables us to study the dynamics of the central spin, and more interestingly showing the fidelity of the spin qubit to investigate the non-Markovian character in the dynamical decoherence of open quantum systems. Motivated by this consideration, in this paper, we consider a fluctated environment with a full polarized nuclear spins to study a decoherence effect at the qubit system which can be solved exactly.
The paper is organized as follows. In Sec.(II) we discuss the best candidate monolayer for a spin qubit. To obtain realistic values of the parameters appearing in the theory we have performed density functional theory (DFT) calculations. Next, the TMD-QD Hamiltonian is given with an external magnetic field perpendicular to the quantum dot is considered, and numerical solutions to the necessary external field strengths at a given quantum dot radius are shown at which a spin-degenerate state within a given valley is expected. In Sec.(III) we derive an effective Hamiltonian describing this single electron spin and $\mathcal{N}$ spin-nuclei in interaction. Our approach is based on deriving an appropriate exact non-Markovian time-convolutionless (TCL) master equation (ME) describing the evolution of the qubit system. We further apply this result to quantify the fidelity loss that the noise induces, and Sec.(IV) concludes with a discussion about our finding.

A huge effort is underway to develop semiconductor nanostructures as low noise hosts for qubits. To provide the possibility for Transition Metal Dichalcogenide (TMD) monolayer candidates to use as a host for a spin qubit, the two spin states ($\ket{\uparrow},\ket{\downarrow}$) selected for the desired qubit need to be degenerate, or tuneable by some external influence (Magnetic field) into a degeneracy. Understanding the band structure and external field replies is a requirement for achieving qubits with carriers in TMD ML. However, the strong intrinsic spin–orbit coupling (SOC) within this material (TMD) exhibits the opposite influence and makes this a non-trivial task. Indeed, the large splitting occurred of the spin states, $\sim meV$, within the same valley ($\mathcal{K}(\mathcal{K}^{\prime}$)) in the conduction band (CB), means that a single electron within a quantum dot (QD) in TMD will not naturally explain the required degeneracies wanted for a spin qubit. Favorably, the band crossing seen in the spin-resolved CB structures in MoS₂-ML which submit that it is reasonable to accomplish spin degeneracy localized within a given valley ($\mathcal{K}(\mathcal{K}^{\prime}$)), see Fig. (1). Consequently, Such spin-degenerate regimes allow the possibility of realizing the desired spin qubits in the ML-MoS₂ quantum dot [3].

Electron spins can be manipulated via external controls by applying an external magnetic field and have been used as qubits, as shown in Sec. (II). For this case, the electron wavefunction is localized inside a quantum dot with $R=26$ nm. The effective Hamiltonian for this single electron spin interacts with a bath of a $\mathcal{N}$ spin-I₀ nuclei, through the contact hyperfine (HF) interaction, in a magnetic field b_z along the z-axis is, ( setting $\hbar=1$),

$$\mathcal{H}_{\mathrm{tot}}=\omega_{s}S_{z}+\omega_{k}I_{z}+\mathbf{h}.\mathbf{S}$$

(5)

where $\mathbf{S}=\left(S_{x},S_{y},S_{z}\right)$ is the electron spin operator. $\omega_{s}=g_{sp}\mu_{B}b_{z}$ $\left(\omega_{k}=g_{I_{k}}\mu_{N}b_{z}\right)$ is the electron (nuclear) Zeeman splitting in a magnetic field $b_{z}$, with effective $g$-factor $g_{sp}\left(g_{I}\right)$ for the electron (nuclei) and Bohr (nuclear) magneton $\mu_{B}\sim 2000\,\mu_{N}\,\left(\mu_{N}\right)$, see Table (2). For magnetic field $b_{z}\sim 1.5\times 10^{-2}$T, $\omega_{s}=1.7\,\mu eV$.

In Eq. (5) we have neglected the anisotropic hyperfine interaction, electron-electron interaction, dipole-dipole interaction between nuclear spins, and nuclear quadrupolar splitting, which may be present for nuclear spin I${}_{0}>1/2$ [20]. The total action by an environment of the $\mathcal{N}$ nuclear spins can be interpreted as a nuclear magnetic field, is the so-called Overhauser field [21]

$$\mathbf{h}=\left(h_{x},h_{y},h_{z}\right)=\sum_{k=0}^{\mathcal{N}-1}A_{k}\mathbf{I}_{k}$$

(6)

where $\mathbf{I}_{k}$ $=\left(I_{k}^{x},I_{k}^{y},I_{k}^{z}\right)$ is the nuclear spin operator at lattice site $k$ at position $\textbf{r}_{k}$. $I_{z}=\Sigma_{k}I_{k}^{z}$ is the total $z$ component of nuclear spin and $A_{k}$ is the associated hyperfine coupling constant. We have also introduced raising and lowering operators $S_{\pm}=S_{x}\pm iS_{y}$, $I_{k}^{\pm}=I_{k}^{x}\pm iI_{k}^{y}$ and the nuclear magnetic field operators $h^{\pm}=h_{x}\pm ih_{y}$. We can rewrite the Eq. (5)

	$\displaystyle\mathcal{H}_{\mathrm{tot}}=\omega_{s}S_{z}+$	$\displaystyle\sum_{k}\omega_{k}I_{k}^{z}+\sum_{k}\frac{A_{k}}{2}\left(S_{+}I_{k}^{-}+S_{-}I_{k}^{+}\right)$		(7)
		$\displaystyle+\sum_{k}A_{k}S_{z}I_{k}^{z}$

The third and the fourth terms in Eq. (7) is the hyperfine contact interaction between the spin electron and the nuclei in the quantum dot, which describe the flip-flop interaction and (longitudinal) Overhauser’s field, giving rise to inhomogeneous broadening and dephasing, indeed the last terme, $\sum_{k}A_{k}S_{z}I_{z}^{k}$, produces an effective magnetic field for the electron $B_{\text{eff }}=b_{z}-\mathcal{B}_{N}$, where [22]

$$\mathcal{B}_{N}=\sum_{k}A_{k}I_{z}^{k}/(g_{sp}\mu_{B})$$

(8)

which results in the well known Overhauser shift. However, when $g_{sp}\mu_{B}B_{\text{eff }}\ll(\sum_{k}A_{k}/2\left(S_{+}I_{k}^{-}+S_{-}I_{k}^{+}\right))$, spin exchange becomes the dominate effect [23]. The strength $A_{k}$ is determined by the electron density at the site of nuclei [24]

$$A_{k}=A^{i_{k}}v_{0}\left|\psi_{0,0}^{\mathcal{K^{{}^{\prime}}},s}\left(\mathbf{r}_{k}\right)\right|^{2}$$

(9)

corresponds to the one-electron hyperfine interaction with the nuclear spin at site k with position ${r}_{k}$. Here, $v_{0}=\sqrt{3}\,a^{2}/4$ is the (two-dimensional) volume of a crystal unit cell containing one nucleus, a is the lattice parameter. Indeed, in QD with radius R includes $\mathcal{N}_{tot}=\pi R^{2}/v_{0}$ nuclei in total. For MoS₂ QD of size $R=26\,nm$ with lattice constant given by a = 3.19 $\AA$, there are $\mathcal{N}_{tot}\sim 10^{4}$ nuclei within the dot. However, only the isotopes that allowing non zero nuclear spin that’s will part of these nuclei, ⁹⁵Mo, ⁹⁷Mo and ³³S. Furthermore, the concentration of ³³S is negligible compared to that of Mo isotopes, and the decoherence of the electron spin mainly originates from the presence of ⁹⁵Mo and ⁹⁷Mo nuclear spins [25]. This leads to the number of nuclear spins within the QD, $\mathcal{N}=(\sum_{i}\nu_{i})\mathcal{N}_{tot}$, where $\nu_{i}$ is the natural abundance for different nuclear isotopic species $i$. $\psi_{0,0}^{\mathcal{K^{{}^{\prime}}},s}(\mathbf{r_{k}})$ is the envelope wave function of the localized electron, and $A^{i_{k}}$ is the total hyperfine coupling constant to a nuclear spin of species $i_{k}$ at site k, where is given by [21],

$$A^{i_{k}}=-\frac{\mu_{0}}{4\pi}\cdot\frac{8\pi}{3}\gamma_{S}\gamma_{i_{k}}\left|u_{i_{k}}\right|^{2}$$

(10)

where, $\mu_{0}$ is the vacuum permeability. $u_{i_{k}}$ is the amplitude of the periodic part of the Bloch function at the position of the nucleus of $i_{k}$ species. $\gamma_{S}$ is the gyromagnetic ratio of free electrons, and its value is always negative. However, $\gamma_{i_{k}}$ which is the nuclear gyromagnetic ratio can take either sign. This fact causes the hyperfine coupling constant $A^{i_{k}}$ to be either positive or negative. For convenience, in a material containing several different nuclear isotopic species $i_{k}$ we define an average hyperfine coupling constant. Here, we take the Root mean square(RMS) average [20, 26]

$$A=\sqrt{\sum_{i}\nu_{i}\left(A^{i}\right)^{2}}$$

(11)

For MoS₂, the hyperfine constant estimated for ⁹⁵Mo and ⁹⁷Mo is $A^{{}^{95}\text{Mo}}$=$A^{{}^{97}\text{Mo}}$ = - 0.57 $\mu$eV [27], by using these coupling constants with the abundances listed in Table (3) gives an RMS coupling strength A = 0.29 $\mu$eV. To be more efficient, in order to show the isotopes that have the main contrubition to the decoherence due to the hyperfine interaction in MoS₂, we plot the decoherence rate $1/T_{2}=\Gamma=\sum_{i}\Gamma_{i}$, See Fig. 4, may be defined as [20]

$$\Gamma_{i}=\frac{1}{T_{2}^{i}}=\nu_{i}^{2}\frac{\pi}{3}\left(\frac{I_{i}\left(I_{i}+1\right)A^{i}}{3\,\omega_{s}}\right)^{2}\frac{A^{i}}{\mathcal{N}}$$

(12)

where $\Gamma_{i}$ is the contribution from flip-flops between nuclei of the common species i. The quadratic dependence on isotopic abundance $\nu_{i}$, shown in Eq. (12), is particularly an important factor in the decoherence rate. Due to this dependence, electron spins in MoS₂, where Mo has two naturally occurring isotopic species, whereas S has only one, will show a decay mostly due to flip-flops between Mo spins, ⁹⁵Mo notably. Significantly, we note that the relatively large flip-flop rates for Mo isotopes, as a result of a large nuclear spin $5/2$ and isotopic natural abundance, respectively.

Now, we consider a localized electron in its orbital ground state. This can be seen, for instance, the wavefunction for the ground state of a single electron isolated in the conduction band (CB) about the valley $\mathcal{K}$ and $\mathcal{K^{{}^{\prime}}}$ under a magnetic field for a quantum dot (QD), with parabolic confinement, in monolayer MoS₂, which is given by,

$$\psi_{0,0}^{\mathcal{K^{{}^{\prime}}},s}(r_{k})=\frac{1}{\sqrt[4]{\pi}\ell_{0}}exp(-\frac{1}{2}(\frac{r_{k}}{\ell_{0}})^{2})\chi_{s}(\sigma_{s})$$

(13)

where, $\ell_{0}$ is the effective length scale, equal to the magnetic length $\ell_{B}=(\hbar/(eB_{z}))^{2}$ in the absence of the confining potential ($\omega_{0}\rightarrow 0$), giving by,

$$\ell_{0}=\sqrt{\frac{\hbar}{2m_{eff}^{\mathcal{K^{{}^{\prime}}},s}\Omega_{\mathcal{K^{{}^{\prime}}},s}}}$$

(14)

Due to our selection in Sec. (II) for the desired spin qubit, $\tau$ would take the value -1 (valley $\mathcal{K}^{\prime}$), therefore the ground state wavefunction that will be considerd for this work is $\psi^{\mathcal{K}^{\prime},s}_{0,0}(r_{k})$.

Adding and subtracting $\sum_{k}A_{k}I_{k}^{z}/2$ to the total Hamiltonian

$$\begin{split}\mathcal{H}_{\mathrm{tot}}&=\omega_{s}S_{z}+\sum_{k}\omega_{k}I_{k}^{z}+\sum_{k}\frac{A_{k}}{2}\left(S_{+}I_{k}^{-}+S_{-}I_{k}^{+}\right)+\sum_{k}A_{k}S_{z}I_{k}^{z}+\sum_{k}\frac{A_{k}}{2}I_{k}^{z}-\sum_{k}\frac{A_{k}}{2}I_{k}^{z}\\ &=\underbrace{\omega_{s}S_{z}+\sum_{k}(\omega_{k}-\frac{A_{k}}{2})I_{k}^{z}}_{\mathcal{H}_{0}}+\underbrace{\sum_{k}\frac{A_{k}}{2}\left(S_{+}I_{k}^{-}+S_{-}I_{k}^{+}\right)+\sum_{k}A_{k}(S_{z}+\frac{1}{2})I_{k}^{z}}_{\mathcal{H}_{1}}\\ &=\mathcal{H}_{0}+\mathcal{H}_{1}\end{split}$$

(15)

This Hamiltonian (15) can be split in two main part, an unperturbed part (longitudinal) noted $\mathcal{H}_{0}$ which consisting of all Zeeman terms and a perturbative part (transverse) $\mathcal{H}_{1}$ containing the virtual flip-flop processes of the hyperfine interaction (HI).
In the interaction picture with respect to the Hamiltonian $\mathcal{H}_{0}$

	$\displaystyle\mathcal{H}_{tot}^{I}$	$\displaystyle=e^{i\mathcal{H}_{0}t}\mathcal{H}_{1}e^{-i\mathcal{H}_{0}t}$		(16)
		$\displaystyle=\underbrace{S_{+}(t)h^{-}(t)+S_{-}(t)h^{+}(t)}_{\text{Transversal Hyperfine Term}}+\underbrace{|1\rangle\langle 1|h^{z}}_{\text{Longitudinal}}$

where

	$\displaystyle h^{\pm}(t)=\sum_{k}\frac{A_{k}}{2}I_{k}^{\pm}e^{\pm i\left(\omega_{k}-A_{k}/2\right)t},$		(17)
	$\displaystyle S_{\pm}(t)=S_{\pm}e^{\pm i\omega_{s}t},$
	$\displaystyle h^{z}=\sum_{k}A_{k}I_{k}^{z}$

Note that the last term in Eq. (16) is equivalent to $\sum_{k}A_{k}(S_{z}+\frac{1}{2}\mathds{1}_{2})I_{k}^{z}$, where $S_{z}=(|1\rangle\langle 1|-|0\rangle\langle 0|)/2$, with $|1\rangle=\left|\mathcal{K}^{\prime}\downarrow\right\rangle(|0\rangle=\left|\mathcal{K}^{\prime}\uparrow\right\rangle)$ is the down (upper) state of the central spin. Indeed, $S_{z}+\frac{1}{2}\mathds{1}_{2}=|1\rangle\langle 1|$. $\mathds{1}_{2}$ is the $2\times 2$ identity matrix. The transversal hyperfine term results in the off-resonant transitions between the system and environmental spins, while the longitudinal hyperfine term provides additional contributions to the energy splitting.

The time evolution of the combined system, consisting of the electron spin and $\mathcal{N}$ nuclear spins, which is given by the action of the total Hamiltonian $\mathcal{H}_{tot}^{I}$ in Eq. (16), is described as the following. To begin with, when $t<0$ we assume that the electron spin and the nuclear system are decoupled, and both of them prepared independently in the states described by the density operators $\rho_{S}(0)$ and $\rho_{E}(0)$, respectively. At $t$ $=0,$ the electron and nuclear spin system are brought into contact over a switching time scale $\tau_{sw}\ll 2\pi\hbar/\left|\omega_{s}-\omega_{k}+A\right|$ [10], which is sufficiently small, where $\left|\omega_{s}-\omega_{k}+A\right|$ is the largest energy scale in this problem. The state of the entire system, described by the total density operator $\rho(t)$, where it’s given at $t=0$,

$$\rho(0)=\rho_{S}(0)\otimes\rho_{E}(0)$$

(18)

with $\rho_{S}=\operatorname{Tr}_{E}(\rho)$ and $\rho_{E}=\operatorname{Tr}_{S}(\rho)$, this is called the reduced density matrix of the subsystem S and E, respectively. The evolution of the density operator $\rho(t)$ for $t\geq 0$ is governed by the Hamiltonian $\mathcal{H}^{I}_{tot}(t)$ for an electron spin coupled to an environment of nuclear spins.

This model that we consider in this work is similar to previous research studies of an electron spin confined to a Gallium arsenide(GaAs) QD [10] and Graphene QD [11], but there are MoS₂ specific properties, that lead to new physics. Given that the natural abundance $\nu_{i}$ of spin-carrying isotopes is small for molybdenum Mo and sulfur S, hence only $\mathcal{N}$ of all atoms $N_{tot}$ within the MoS₂ QD carry spin. However, in semiconducting materials such as GaAs, all isotopes possess a spin. To highlight these differences, we do a comparison of the most important characteristics of MoS₂, graphene, and GaAs which is given in Table (3). Indeed, the hyperfine interaction coupling constant A${}_{\text{MoS${}_{2}$}}$, is about approximately the same magnitude as A${}_{\text{Graphene}}$ and it’s about two orders of magnitude smaller than the constant A${}_{\text{GaAs}}$ in GaAs which even further reduce the nuclear magnetic field by the same amount. Moreover, The relatively small hyperfine energy, $\mathcal{E}_{HF}=\hbar/\tau_{HF}$ (in case of MoS₂ is on average of order $\mathcal{E}_{HF}\thickapprox 10^{-12}$eV), in MoS₂ increases the timescale $\tau_{HF}$, which is depending on both the hyperfine strength A and the isotopic abundance $\nu_{i}$, of this interaction significantly as compared to GaAs and graphene. The switching time scale $\tau_{sw}$ for several materials are listed in Table (3) as well.

Next, We drive into the method that used along side with this work. We First suppose at $t=0$, that the total system (electron and nuclear) describe by

$$\psi(0)=\psi_{S}(0)\otimes\psi_{E}(0)$$

(19)

where $\ket{\psi_{S}(0)}$, is the initial state for the central spin and $\ket{\psi_{E}(0)}$ is the initial state for the nuclear spin bath where we start with a perfectly polarized nuclear ensemble, as shown in Fig. (5). The Eq. (19) can be written in the subspace spanned by the bases $\{|0\rangle\otimes|0\rangle_{E},|1\rangle\otimes|0\rangle_{E},|0\rangle\otimes\left(I^{+}_{k}|0\rangle_{E}\right)\}$ as following:

	$\displaystyle\mid\psi(0)\rangle=$	$\displaystyle c_{0}|0\rangle\otimes|0\rangle_{E}+c_{1}(0)|1\rangle\otimes|0\rangle_{E}$		(20)
		$\displaystyle+\sum_{k}c_{k}(0)|0\rangle\otimes\left(I^{+}_{k}|0\rangle_{E}\right)$

where, $|0\rangle_{E}=\otimes_{k}^{{\mathcal{N}}}|0\rangle=|\underbrace{0,\cdots,0}_{{\mathcal{N}}}\rangle$ denote the vacuum state of the bath. Noting that, $c_{k}(0)=0\,\forall\,k$. This means that the environment is in the vacuum state initially. The time evolution of the total system, $t>0$, can be written by the following expression

	$\displaystyle\mid\psi(t)\rangle=$	$\displaystyle c_{0}|0\rangle\otimes|0\rangle_{E}+c_{1}(t)|1\rangle\otimes|0\rangle_{E}$		(21)
		$\displaystyle+\sum_{k}c_{k}(t)|0\rangle\otimes\left(I^{+}_{k}|0\rangle_{E}\right)$

where we have used the normalization condition $\mid c_{0}\mid^{2}+\mid c_{1}(t)\mid^{2}+\sum_{k}\mid c_{k}(t)\mid^{2}=1$. Let us introduce the states

$$\ket{\psi_{0}}=|0\rangle\otimes|0\rangle_{E},\>\ket{\psi_{1}}=|1\rangle\otimes|0\rangle_{E},\>\ket{\psi_{k}}=|0\rangle\otimes|k\rangle_{E}$$

(22)

where $|k\rangle_{E}=I^{+}_{k}|0\rangle_{E}=\ket{0_{1},\cdots,0_{k-1},1_{k},0_{k+1,\cdots}}$ denotes the state with only one nuclear spin in site k. We can now express the time evolved state (21) as:

$$\mid\psi(t)\rangle=c_{0}\ket{\psi_{0}}+c_{1}(t)\ket{\psi_{1}}+\sum_{k}c_{k}(t)\ket{\psi_{k}}$$

(23)

The amplitude $c_{0}$ is constant since $\mathcal{H}^{I}_{tot}(t)\ket{\psi_{0}}=0$, while the amplitudes $c_{1}(t)$ and $c_{k}(t)$ are time dependent. The time development of these amplitudes is governed by a system of differential equations which is obtained by the Hamiltonian in Eq. (16) and the Schrödinger equation,

$$\begin{split}i\partial_{t}|\psi(t)\rangle&=\dot{c}_{1}(t)\left|\psi_{1}\right\rangle+\sum_{k}\dot{c}_{k}(t)\left|\psi_{k}\right\rangle\\ &=\mathcal{H}_{tot}^{I}|\psi(t)\rangle=S_{+}(t)h^{-}(t)+S_{-}(t)h^{+}(t)+|1\rangle\langle 1|h^{z}\\ &\times\left(c_{0}\left|\psi_{0}\right\rangle+c_{1}(t)\left|\psi_{1}\right\rangle+\sum_{k}c_{k}(t)\left|\psi_{k}\right\rangle\right)\\ &=\sum_{k}\frac{A_{k}}{2}e^{-i\left(\omega_{s}-\omega_{k}+A_{k}/2\right)t}c_{1}(t)\ket{\psi_{k}}\\ &+\left(hc_{1}(t)+\sum_{k}c_{k}(t)\frac{A_{k}}{2}e^{-i\left(\omega_{k}-\omega_{s}-A_{k}/2\right)t}\right)\ket{\psi_{1}}\end{split}$$

(24)

where $h=\left\langle h^{z}(t)\right\rangle=\left\langle\psi(t)\left|h_{z}\right|\psi(t)\right\rangle$. Multiplying by $\bra{\psi_{1}}$ and $\bra{\psi_{k}}$ gives us two coupled differential equations for the amplitudes $c_{k}(t)$ and $c_{1}(t)$:

$$\begin{array}[]{c}\frac{d}{dt}c_{1}(t)=ihc_{1}(t)-i\sum_{k}\frac{A_{k}}{2}e^{i\left(\omega_{s}-\omega_{k}+\frac{A_{k}}{2}\right)t}c_{k}(t)\\ \frac{d}{dt}c_{k}(t)=-i\frac{A_{k}}{2}e^{-i\left(\omega_{s}-\omega_{k}+\frac{A_{k}}{2}\right)t}c_{1}(t)\end{array}$$

(25)

By integrating Eq. (25), the coefficient $c_{k}(t)$ can be formally written as

$$c_{k}(t)=-i\frac{A_{k}}{2}\int_{0}^{t}ds\,e^{-i\left(\omega_{s}-\omega_{k}+\frac{A_{k}}{2}\right)s}c_{1}(s)$$

(26)

Substituting it into Eq.(25), we obtain an exact time-convolution dynamical equation for the central spin

$$\frac{d}{dt}c_{1}(t)=ihc_{1}(t)-\int_{0}^{t}ds\left\langle h^{-}(t)h^{+}(s)\right\rangle_{E}e^{i\omega_{s}(t-s)}c_{1}(s)$$

(27)

where, $\left\langle h^{-}(t)h^{+}(s)\right\rangle_{E}e^{i\omega_{s}(t-s)}=\sum_{k}\left(A_{k}/2\right)^{2}e^{i\left(\omega_{s}-\omega_{k}+A_{k}/2\right)(t-s)}$ is the associated memory kernel function given by a two-point correlation function of the bath. $\langle\cdots\rangle_{E}$, is the expectation values over the state of the environment $\rho_{E}$, where $\rho_{E}(t)=\left(\ket{0}\bra{0}\right)_{E}$ is the vacuum state of the bath. The correlation function in the continuum limit assumes the form $\int d\omega\mathcal{J}(\omega)\,e^{i\left(\omega_{s}+A_{k}/2-\omega\right)(t-s)}$, with $\mathcal{J}(\omega)=\sum_{k}\left(A_{k}/2\right)^{2}\delta(\omega-\omega_{k})$ is the spectral density of the bath given by sum of (coupling strength)² × (density of modes), which is therefore simply the Fourier transform of the correlation function $\left\langle h^{-}(t)h^{+}(s)\right\rangle_{E}e^{i\omega_{s}(t-s)}$. The typical representative environments are described by the Lorentzian-type spectral functions, $\mathcal{J}(\omega)=\gamma_{0}\lambda^{2}/\left(2\pi\left(\omega_{s}-\omega\right)^{2}+\lambda^{2}\right)$. Here, the parameter $\lambda$ defines the width of the Lorentzian spectral density, indeed its the measure of the memory capacity or non-Markovianity of the environment [28] and is connected to the bath correlation time $\tau_{E}=1/\lambda$. On the other hand, $\gamma_{0}$ measures the strength of coupling between the qubit and its environment and hence the system characteristic time $\tau_{S}=1/\gamma_{0}$ denotes the relaxation time. By taking a Lorentzian spectral density in resonance with the transition frequency of the qubit we find an exponential two-point correlation function, expressed as

$$\left\langle h^{-}(t)h^{+}(s)\right\rangle_{E}e^{i\omega_{s}(t-s)}=\frac{1}{2}\gamma_{0}\lambda\mathrm{e}^{-\lambda\left|t-s\right|}$$

(28)

We define the function $\tilde{G}(t)=G(t)e^{-iht}$, where $c_{1}(t)=G(t)c_{1}(0)$, which can show that it defined as the solution of the integro–differential equation,

$$\partial_{t}\tilde{G}(t)=-\int_{0}^{t}ds\left\langle h^{-}(t)h^{+}(s)\right\rangle_{E}e^{i\omega_{s}(t-s)}e^{-ih(t-s)}\tilde{G}(s)$$

(29)

with initial condition $\tilde{G}(0)=G(0)=1$. Generally, $G(t)$ can be solved to give the exact solution by Laplace transform. Indeed, substituting Eq. (28) into Eq. (29) we obtain,

$$G(t)=e^{-\lambda t/2}\left(\frac{\sinh\left(\frac{\lambda\chi}{2}t\right)}{\chi}+\cosh\left(\frac{\lambda\chi}{2}t\right)\right)$$

(30)

Where $\chi=\sqrt{1-2(\gamma_{0}/\lambda)}$.
In order to get a master equation in differential form with a generator local in time, that is, the Time Convolutionless (TCL) master equation. We first give the exact time evolution mapping [16], which transforms the initial states into the states at time t

	$\displaystyle\Phi(t):\rho(0)\rightarrow\rho(t)$	$\displaystyle=\operatorname{Tr}_{E}\left\{\mid\psi(t)\rangle\langle\psi(t)\mid\right\}$		(31)
		$\displaystyle=\Phi(t)\rho(0),\>t\geq 0$

Then due to Eq. (21), the density matrix $\rho(t)$ it is expressed as

$$\begin{split}\rho(t)&=\operatorname{Tr}_{E}\{|\Psi(t)\rangle\langle\Psi(t)|\}\\ &=\left(\begin{array}[]{cc}{\rho}_{11}(t)&{\rho}_{10}(t)\\ {\rho}_{01}(t)&{\rho}_{00}(t)\end{array}\right)\\ &=\left(\begin{array}[]{cc}|G(t)|^{2}\rho_{11}(0)&G(t)\rho_{10}(0)\\ G^{\star}(t)\rho_{01}(0)&\rho_{00}(0)+\left(1-|G(t)|^{2}\right)\rho_{11}(0)\end{array}\right)\\ &=\left(\begin{array}[]{cc}|{c}_{1}(t{)|}^{2}&{c}_{0}^{\ast}{c}_{1}(t)\\ {c}_{0}{c}_{1}^{\ast}(t)&1-|{c}_{1}(t{)|}^{2}\end{array}\right)\end{split}$$

(32)

where $\rho_{ij}(t)=\langle i|\rho(t)|j\rangle$ for i, j=0,1. We can construct the exact TCL equation, $\partial_{t}\rho(t)=$ $\mathcal{K}_{\mathrm{TCL}}(t)\rho(t)$, exploiting the introduced relations [16, 29], we can introduce a time-local generator

$$\mathcal{K}_{\mathrm{TCL}}(t)=\dot{\Phi}(t)\Phi^{-1}(t)$$

(33)

We can obtain an exact TCL master equation to second order in the interaction picture [28, 30, 31]

	$\displaystyle\partial_{t}\rho(t)$	$\displaystyle=\mathcal{K}_{\mathrm{TCL}}(t)\rho(t)$		(34)
		$\displaystyle=-\frac{i}{2}\varepsilon(t)\left[S_{+}S_{-},\rho(t)\right]$
		$\displaystyle+\gamma(t)\left[S_{-}\rho(t)S_{+}-\frac{1}{2}\{S_{+}S_{-},\rho(t)\}\right]$

Consider the time derivative of $\rho(t)$

$${\partial}_{t}\rho(t)=\left(\begin{array}[]{cc}\Re(\frac{{\dot{c}}_{1}(t)}{{c}_{1}(t)})|{c}_{1}(t{)|}^{2}&(\frac{{\dot{c}}_{1}(t)}{{c}_{1}(t)}){c}_{0}^{\ast}{c}_{1}(t)\\ (\frac{{\dot{c}}_{1}^{\ast}(t)}{{c}_{1}^{\ast}(t)}){c}_{0}{c}_{1}^{\ast}(t)&-\Re(\frac{{\dot{c}}_{1}(t)}{{c}_{1}(t)})|{c}_{1}(t{)|}^{2}\end{array}\right)=\left(\begin{array}[]{cc}\gamma(t)\left|c_{1}(t)\right|^{2}&\left(\frac{i}{2}\varepsilon(t)-\frac{1}{2}\gamma(t)\right)c_{0}c_{1}^{*}(t)\\ \left(-\frac{i}{2}\varepsilon(t)-\frac{1}{2}\gamma(t)\right)c_{0}^{*}c_{1}(t)&-\gamma(t)\left|c_{1}(t)\right|^{2}\end{array}\right)$$

(35)

The dynamics of the exact TCL master equation, parameterized by $\epsilon(t)$ plays the role of a time-dependent Lamb shift, which is induced by coupling to the noisy environments, and $\gamma(t)$ plays the role of a time-dependent decay rate(decoherence rate). In addition, the decay rate $\gamma(t)$ can have negative values, which means that the dynamics of the system exhibit a strong non-Markovian behavior [30]. This parameters can be written as follows,

	$\displaystyle\gamma(t)+i\varepsilon(t)$	$\displaystyle=-2\left[\frac{\dot{G}(t)}{G(t)}\right],\>\varepsilon(t)=-2\Im\left[\frac{\dot{G}(t)}{G(t)}\right],$		(36)
		$\displaystyle\gamma(t)=-2\Re\left[\frac{\dot{G}(t)}{G(t)}\right]$

We use the fidelity [32, 28] $\mathcal{F}(t)=\sqrt{\langle\psi(0)|\rho(t)|\psi(0)\rangle}$ to measure the decoherence dynamics of the central spin. Indeed, the coherence of spin qubits is highly affected by the nuclear spins of the host material and their hyperfine coupling to the electron spin. When the system is prepared in the initial state $\ket{\psi(0)}=\ket{1}$, the fidelity can be written as

$$\mathcal{F}(t)=\sqrt{\left|c_{1}(t)c_{1}(0)\right|^{2}}=|G(t)|=|\tilde{G}(t)|$$

(37)

The fidelity illustrated in Fig. (6) and Fig. (6), can be considered in two main regimes. In the weak coupling case, which reflects that the decoherence dynamics of the quantum system is Markovian, $\lambda>2\gamma_{0}$, one has $\chi\,\epsilon\,\mathbb{R}$ ($\chi=\sqrt{1-2\gamma_{0}/\lambda}$) so that $G(t)$ is always positive. Indeed, when there is no more exchange between the qubit and his environment the coupling strength $\gamma_{0}\rightarrow 0$ the fidelity $\lim_{\gamma_{0}\to 0}\mathcal{F}(t)\sim 1$. However, in the strong coupling case, $\lambda<2\gamma_{0}$ one has $\chi\,\epsilon\,\mathrm{i}\mathbb{R}$, so that $G(t)$ oscillates between positive and negative values going through zero. This means the fidelity $\mathcal{F}(t)$ will then decay with an oscillating. Indeed in Fig. (6) by observing the contrast between the light blue and dark blue areas when $\lambda/\gamma_{0}<1$ , in another word the fidelity has a revival, which means the quantum information flow bounces from the spin bath back to the qubit system. This is the signature of non-Markovian behavior.