Coherence as a measure for non-Markovianity and its applications to color centers in diamond

TL;DR

This paper introduces an extended coherence-based measure for non-Markovianity applicable to structured phonon baths in color centers in diamond, providing a practical tool for analyzing quantum system-environment interactions at various temperatures.

Contribution

It proposes a generalized measure of non-Markovianity that overcomes limitations of previous methods and is adaptable to coherent dynamical maps in complex quantum systems.

Findings

The measure is well-behaved at arbitrary temperatures.

It reveals how non-Markovianity is influenced by bulk and localized phonon modes.

The method can eliminate false signs of non-Markovianity in coherent dynamical maps.

Abstract

The degree of non-Markovianity of a continuous bath can be quantified by means of the coherence. This simple measure is experimentally accessible through Ramsey spectroscopy, but it is limited to incoherent dynamical maps. We propose an extension of this measure and discuss its application to color centers in diamond, where the optical coherence between two orbital states is affected by interactions with a structured phonon bath. By taking realistic phonon spectral density functions into account, we show that this measure is well-behaved at arbitrary temperatures and that it provides additional insights about how non-Markoviantiy is affected by the presence of both bulk and quasi-localized phonon modes. Importantly, with only a little overhead the measure can be adapted to eliminate the false signs of non-Markovianity from coherent dynamical maps and is thus applicable for a large class of systems modeled by the spin-boson Hamiltonian.

Equations

$${\mathcal{N}}_C(T)=1+\frac{{\int }_0^{\infty }\dot{C}(\tau )d\tau }{{\int }_0^{\infty }|\dot{C}(\tau )|d\tau },0\le {\mathcal{N}}_C(T)\le 1.$$

$$H=-\frac{\Delta }{2}{\sigma }_z+\frac{\Omega }{2}{\sigma }_x+\sum _k{\omega }_k{a}_k^{†}{a}_k+{\sigma }_z\sum _k{g}_k\left(a_k+a_k^{†}\right),$$

$$\frac{d{\rho }_s}{dt}=-\frac{\gamma (t)}{2}({\rho }_s(t)-{\sigma }_z{\rho }_s(t){\sigma }_z),$$

$$\gamma (t)={\int }_0^{\infty }\frac{J(\omega )}{\omega }\coth \left(\frac{\omega }{2k_{B}T}\right)\sin (\omega t)d\omega ,$$

$$J_{bulk}(\omega )$$

$$J_{loc1}(\omega )$$

$$J_{loc2}(\omega )$$

$${\gamma }_{bulk}^{\downarrow }(t,d)$$

$${\gamma }_{loc1}^{\downarrow }(t)$$

$${\gamma }_{loc2}^{\downarrow }(t)$$

$$H=\frac{1}{2}{\omega }_0{\sigma }_z+\sum _k{\omega }_k{a}_k^{†}{a}_k+{\sigma }_z\sum _k{g}_k\left(a_k+a_k^{†}\right),$$

$$\frac{d⟨{\sigma }_{\pm }(t)⟩}{dt}=\left[\pm \frac{i{\omega }_0}{2}-\gamma (t)\right]⟨{\sigma }_{\pm }(t)⟩,$$

$$\frac{d{\rho }_s(t)}{dt}=-i[H_s,{\rho }_s(t)]-\frac{\gamma (t)}{2}\left(⟨{\sigma }_x(t)⟩{\sigma }_x+⟨{\sigma }_y(t)⟩{\sigma }_y\right),$$

$${\sigma }_z{\rho }_s(t){\sigma }_z=\frac{1}{2}\left(\text{mathds}1-⟨{\sigma }_x(t)⟩{\sigma }_x-⟨{\sigma }_y(t)⟩{\sigma }_y+⟨{\sigma }_z(t)⟩{\sigma }_z\right).$$

$$\frac{d{\rho }_s(t)}{dt}=-i[H_s,{\rho }_s(t)]-\frac{\gamma (t)}{2}\left({\rho }_s(t)-{\sigma }_z{\rho }_s(t){\sigma }_z\right).$$

$${\gamma }_{loc1}(t)=\frac{J_0\Gamma }{2}{\int }_0^{\infty }\frac{{\omega }^2\coth \left(\frac{\omega }{2k_{B}T}\right)\sin (\omega t)}{{\left(\frac{\omega }{{\omega }_{loc}}+1\right)}^2(\omega -{\omega }_{loc}{)}^2+(\Gamma /2{)}^2}d\omega .$$

$${\gamma }_{loc1}(t)$$

$$f(u)=\frac{{\left(u+{\omega }_{loc}\right)}^2\coth \left(\frac{u+{\omega }_{loc}}{2k_{B}T}\right)}{{\left(u/{\omega }_{loc}+2\right)}^2}\frac{\Gamma /2}{u^2+{\left(\Gamma /2\right)}^2}.$$

$$f(u)\approx \frac{1}{4}{\omega }_{loc}^2\coth \left(\frac{{\omega }_{loc}}{2k_{B}T}\right)\frac{\Gamma /2}{u^2+{\left(\Gamma /2\right)}^2},$$

$${\gamma }_{loc1}(t)\approx \frac{\pi J_0{\omega }_{loc}^2}{4}\coth \left(\frac{{\omega }_{loc}}{2k_{B}T}\right)\sin ({\omega }_{loc}t)e^{-\Gamma t/2},$$

$${\gamma }_{loc1}^{\uparrow }(t)\approx \left(\frac{2k_{B}T}{{\omega }_{loc}}\right){\gamma }_{loc1}^{\downarrow }(t).$$

$${\mathcal{N}}_{BLP}={\int }_{dD/dt>0}\frac{dD}{d\tau }d\tau ,D(t)=\frac{1}{2}\text{mbox}Tr|{\rho }_1(t)-{\rho }_2(t)|,$$

$$\rho_{1}(0)={1\over 2}\left(\begin{array}[]{ll}1&1\\ 1&1\end{array}\right),\quad\quad\rho_{2}(0)={1\over 2}\left(\begin{array}[]{rr}1&-1\\ -1&1\end{array}\right).$$

$$\dot{{\rho }_s}(t)=-i[H_s,{\rho }_s(t)]+D{\rho }_s(t),$$

$$|1⟩$$

$$|2⟩$$

$$D{\rho }_s(t)=\sum _{\xi =0,\pm {\omega }_0}\frac{{\gamma }_{\xi }(t)}{2}\left(L_{\xi }{\rho }_s(t)L_{\xi }^{†}-\frac{1}{2}\{L_{\xi }^{†}{L}_{\xi },{\rho }_s(t)\}\right).$$

$${\gamma }_{\xi }(t)$$

$${\dot{\rho }}_{ee}(t)$$

$${\dot{\rho }}_{eg}(t)$$

$${\rho }_{ee}(t)$$

$${\rho }_{eg}(t)$$

$$⟨S_z{⟩}_{({\rho }_{ee}=1)}$$

$$⟨S_z{⟩}_{({\rho }_{gg}=1)}$$

Full text

††preprint: APS/123-QED

I Introduction

By employing the principles of superposition and entanglement, quantum systems can outperform their classical counterparts in many applications such as computation, cryptography and high-precision measurements Nielsen ; Giovannetti ; Gross . However, to benefit from this quantum advantage, the systems must be protected from detrimental interactions with the environment using passive isolation as well as active techniques such as error correction or decoupling pulses. To implement efficient error-mitigation schemes it is crucial to have a precise understanding of the underlying system-environment interaction Breuerbook , in particular in realistic non-Markovian settings, where information can flow back from the bath to the system and is not immediately lost. These scenarios have been exploited, for example, for quantum metrology Chin , quantum channels Bylicka or quantum control Reich and led to a growing interest in the questions how different environments can be compared and how non-Markovianity (NM) can be quantified Breuer ; Luo ; Rivas . However, most of the proposed measures for NM so far are based on rather tedious mathematical constructs, such as the maximization of the trace distance between two different initial states Breuer , the addition of auxiliary systems Luo or the complete knowledge of the dynamical map Rivas . This makes proof-of-concept demonstrations and the broad use of such measures for modeling realistic applications very cumbersome.

In this work, we focus on a simple measure for NM, which can be directly derived from the coherence of the system of interest Chanda ; Passos . This quantity is most essential for quantum technology applications and in many cases it can be accessed through Ramsey measurements. To demonstrate its suitability for characterizing realistic system-bath interactions in a physically meaningful way, we discuss in more detail the application of this measure for color centers in diamond. Color centers such as the negatively charged silicon-vacancy (SiV-) and nitrogen-vacancy (NV-) center have attracted wide attention because they can be initialized, controlled and read out with high fidelities Gruber ; Oort ; Becker . These unique properties make them strong candidates for various quantum sensing and quantum information processing application Jero2008 ; Dutt2007 , but their optical properties are still limited by unavoidable interactions with phonons. The influence of a continuum of bulk modes and distinct quasi-localized resonances originates a rich and complex dynamics arising from different NM behavior beyond the extensively studied Ohmic environment Legget1987 . A deeper understanding of the NM dynamics induced by realistic phononic baths will be important for engineering and optimizing defect-phonon interactions in structured reservoirs such as cantilevers Kepesidis2016 , two-dimensional layers Abdi2019 , phonon waveguides Lemonde2018 or phononic crystals Mark2018 ; Nori2019 .

II Measure of non-Markovianity

One natural approach towards witnessing NM is to examine the back-flow of quantum information in terms of the Coherence LoFranco . For a given complete set of basis states $\{|i\rangle\}$, Coherence is commonly defined as $C(t)=\sum_{i\neq j}|\rho_{ij}(t)|$ Baumgratz , where the $\rho_{ij}(t)=\langle i|\rho(t)|j\rangle$ are the matrix elements of the system density operator. Although this definition is not unique, the choice of basis states usually follows naturally from the context and is very often taken as the eigenbasis of the bare system. For a two-level system the definition above reduces to $C(t)=\sqrt{\langle\sigma_{x}\rangle^{2}+\langle\sigma_{y}\rangle^{2}}$, where the $\sigma_{k}$ are the usual Pauli operators. In this case, $C(t)$ can be experimentally attained through standard Ramsey spectroscopy, i.e., by applying two $\pi/2$-pulses separated by a time $t$. A measurement of the final population in the excited state, $P_{e}=(1+\langle\sigma_{x}\rangle)/2$, then provides a measurement of $\langle\sigma_{x}\rangle$, or $\langle\sigma_{y}\rangle$, if an additional rotation between the two pulses is introduced. Similar strategies can also be applied to measure $C(t)$ for higher dimensional systems, where, however, the pulse sequences are slightly more involved.

The Coherence does not increase under incoherent completely positive and trace preserving maps ($\{\Lambda_{ICPTP}\}$) Baumgratz , in other words, it is a monotonically decreasing function in time $C_{\rho}(t)\geq C_{\Lambda_{ICPTP}(\rho)}(t)$. Hence, considering a free evolution that obeys $\Lambda_{t}\in\Lambda_{ICPTP}$, one can build up a measure of NM based on the oscillations in $C(t)$, that are triggered by the coupling to an external reservoir Chanda ; Passos . The measure of NM is then defined as

[TABLE]

Even though $\mathcal{N}_{C}(T)$ is only valid in the presence of incoherent dynamical maps, we show in Sec. V.1 that this constraint can be relaxed and that the measure can be extended to coherent maps, with practical applications on the generalized spin-boson Hamiltonian.

III Optical coherence in color centers

In this section, we investigate in more detail the application of the measure defined in Eq. (1) to the case of SiV*-* and NV*-* centers in diamond. By taking only a single ground state $|g\rangle$ and a single electronically excited state $|e\rangle$ into account, the optical properties of the center in the presence of a phonon reservoir are well-described by the spin-boson Hamiltonian. In the frame rotating with the frequency $\omega_{L}$ of the driving laser, the Hamiltonian reads ($\hbar=1$) Legget ; Weiss

[TABLE]

where $\Delta=\omega_{L}-\omega_{eg}$ is the detuning from the bare transition frequency $\omega_{eg}$ and $\Omega$ is the optical Rabi frequency. In Eq. (2), the $\omega_{k}$ are the phonon frequencies and $a_{k}^{\dagger}$ ($a_{k}$) the corresponding creation (annihilation) operators for the phonon modes. The $g_{k}=\lambda_{e,k}-\lambda_{g,k}$ Ariel2016 denote the effective electron-phonon coupling constants, which arise from the different deformation potentials, $\lambda_{e,k}$ and $\lambda_{g,k}$, in the ground and the excited state.

During the free evolution time $(\Omega=0)$ and assuming that the phonon reservoir is initially in a thermal state $\rho_{ph}$, the dynamics of the reduced color center state is described by the following time local master equation Luczka ; Haikka

[TABLE]

where we have omitted a Hamiltonian contribution $\sim[\sigma_{z},\rho_{s}(t)]$, which does not affect the coherence. Therefore, the system-environment interaction is fully determined by the time-dependent dephasing rate Luczka ($\hbar=1$)

[TABLE]

where $J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$ is the spectral density function (SDF), $k_{B}$ is the Boltzmann constant and $T$ is the reservoir temperature. In general, the SDF $J(\omega)$ satisfies two important properties: i) $J(0)=J(\omega>\omega_{max})=0$ and ii) $J(\omega)>0\;\forall\;\omega\in(0,\omega_{max})$, where $\omega_{max}$ is the largest phonon frequency of the reservoir. The formal solution of the off-diagonal elements of $\rho_{s}(t)$ is given by $\rho_{eg}(t)=\rho_{eg}(0)e^{-\Gamma(t)}$, where $\Gamma(t)=2\int_{0}^{t}\gamma(\tau)\,d\tau$ is a bounded function that satisfies $0\leq\Gamma(t)\leq 4\int_{0}^{\infty}J(\omega)/\omega^{2}\mbox{coth}(\omega/2k_{B}T)\,d\omega$.

IV Phonon spectral density and dephasing rate

From Eq. (4) we see that the dephasing rate $\gamma(t)$, and thus the degree of NM, depends only on the SDF, $J(\omega)$, and the temperature $T$. In the case of color centers or other solid-state emitters, information about the SDF can be obtained from the photoluminescence (PL) spectrum Dietrich2014 ; Ariel2016 , where the coupling to the phonons both reduces the bare resonance (zero-phonon line) and leads to additional phonon-sidebands. The experimental PL spectrum of the SiV*-* center exibits an isotopic shift feature in the prominent and narrow phonon sideband Dietrich2014 . This can be explained by the strong electron-phonon interaction with a quasi-localized phonon mode primarily composed of a large oscillation of the silicon atom Dietrich2014 ; Ariel2016 ; Londero2018 . In addition, lattice vibrations in the bulk lead to a smooth SDF, which typically scales as $\sim\omega^{3}$ for low frequencies in a three-dimensional lattice. The phononic SDF that reproduces the isotopic shift feature of the PL spectrum and the effect of acoustic phonons is phenomenologically given by $J(\omega)=J_{bulk}(\omega)+J_{loc1}(\omega)+J_{loc2}(\omega)$, with Ariel2016

[TABLE]

where $d$ is the dimension of the diamond lattice ($d=3$ in our case). Acoustic phonons are associated with low-energy vibrational excitations where the atoms of the color center are oscillating in phase and therefore experimenting a weak electron-phonon interaction. Thus, these phonons are reasonably well described by an intensity $\alpha$ and a cutoff frequency $\omega_{c}\simeq 1$ THz. However, quasi-localized phonon modes induce out-of-phase oscillations of the defect’s atoms with large amplitude leading to a strong electron-phonon interaction. This type of interactions is usually modeled by Lorentzian-like functions centered around the specific localized phonon frequency $\omega_{loc}$, with a characteristic width $\Gamma$ and an intensity $J_{0}$. For the region in between (from 1 THz to $\sim$14 THz) other vibrational modes participate Gergo2018 , which is captured by the Gaussian contribution $J_{loc2}(\omega)$. For the SiV*-* center, the following parameters of the SDF are chosen to accurately match the PL spectrum obtained from molecular dynamic simulations: $\alpha=0.0275$, $J_{1}=0.0025$ THz*-2*, $\sigma=2.4042$ THz, $\omega_{0}=9.35$ THz, $J_{0}=0.0235$ THz*-1*, $\Gamma=0.8414$ THz, and $\omega_{loc}=15.19$ THz. For the NV- center we do not make this decomposition and use the exact SDF obtained from a detailed first-principles calculation reported in Ref. Alkauskas . In Figure 1 we plotted the normalized spectral density functions of both NV*-* and SiV*-* centers.

By following the same partition as for the SDF, we write for the SiV- center the total dephasing rate defined in Eq. (4) as $\gamma(t)=\gamma_{bulk}(t)+\gamma_{loc1}(t)+\gamma_{loc2}(t)$. In order to better illustrate the behavior of the individual contributions we will first consider the limit of very low temperatures where $\coth(\omega/2k_{B}T)\approx 1$ and

[TABLE]

For a detailed derivation of the dephasing rate $\gamma_{loc1}^{\downarrow}(t)$ see Appendix B. We left $\gamma_{loc2}^{\downarrow}(t)$ in terms of the integral since it involves a rather complicated expression. Note that Eq. (8) is only valid for $d>-1$, which is satisfied in our case ($d=3$). In the opposite, high temperature regime, $\omega/2k_{B}T\ll 1$ and $\coth(\omega/2k_{B}T)\approx 2k_{B}T/\omega$. In this case we obtain $\gamma_{bulk}^{\uparrow}(t)=\left(2k_{B}T/\omega_{c}\right)\gamma_{bulk}^{\downarrow}(t,d-1)$, $\gamma_{loc1}^{\uparrow}(t)\approx\left(2k_{B}T/\omega_{loc}\right)\gamma_{loc1}^{\downarrow}(t)$ and $\gamma_{loc2}^{\uparrow}(t)=2k_{B}TJ_{1}\int_{0}^{\infty}\omega\mbox{exp}(-(\omega-\omega_{0})^{2}/(2\sigma^{2}))\sin(\omega t)\,d\omega$.

Figure 2 shows the expected dephasing rates for the SiV- center in diamond for low and high temperatures. Moreover, based on numerical simulations we set the low (high) temperature regime to be below (above) $T\approx 20$ K ($286$ K). These boundaries have been chosen such that the approximate low- and high-temperature limits for the dephasing rate given above match its exact numerical shape. From the statistical properties of the phonon reservoir, in particular, the mean number of phonons $n(w)=[\mbox{exp}\left(\hbar\omega/k_{B}T\right)-1]^{-1}$, we observe that $n(\omega_{loc})\approx 2$ for $T\approx 286$ K and $\omega_{loc}=15.19$ THz. In other words, the high-temperature regime is defined from the thermal activation of the strong electron-phonon coupling with the quasi-localized phonon mode.

If we individually look at the components of $\gamma(t)$, it is straightforward to notice that each one will get involved in the dynamics at a different temperature. For $\gamma_{bulk}(t)$ in Figure 2-(a) we observe a region where it takes negative values, and above a critical temperature ($T_{c}\approx 1.85$ K) it is always positive. However, $\gamma_{loc1}(t)$ and $\gamma_{loc2}(t)$ show negative values even at room temperature, with an amplitude that is non-negligible as compared to the low-temperature case, see Figures 2-(b) and 2-(c). This negative behavior has previously been connected to memory effects Luo , as we will detail in what follows. Finally, one can notice that $\gamma_{loc1}(t)$ has the leading contribution to the dephasing rate, see Figure 2-(d). For comparison, in Appendix C we also evaluate the dephasing rate for an NV*-* center, which also exhibits negative values at room temperature.

V Non-markovianity in color centers

Let us now discuss the degree of NM of the described phonon environment, as quantified through $N_{C}(T)$. To begin, we first introduce another related measure for NM, which we will use for comparison. Following the procedure in Ref. Hall , where the figure of merit is the canonical decay rate $\gamma^{c}(t)$ when the Lindbladian is written in an orthonormal basis for a $d$-level system $\mathcal{L}[\rho(t)]=-i[H(t),\rho(t)]+\sum_{j=1}^{d^{2}-1}\gamma_{j}^{c}(t)(L_{j}(t)\rho(t)L_{j}^{\dagger}(t)-\{L_{j}^{\dagger}(t)L_{j}(t),\rho(t)\}/2)$, where $\mbox{Tr}[L_{i}^{\dagger}(t)L_{j}(t)]=\delta_{ij}$. Based on this rate, the function $f(t)\equiv\max\{-\gamma^{c}(t),0\}=(|\gamma^{c}(t)|-\gamma^{c}(t))/2$ Hall is introduced. For our particular case where the effect of the environment represented in Eq. (3) only induces a pure dephasing dynamics and the operators involved form an orthogonal basis, we have $\gamma^{c}(t)=2\gamma(t)$ Hall ; Rivas . Therefore, a NM measure ($\mathcal{N}_{\gamma}$) as a function of the reservoir temperature is defined when integrating $f(\tau)$ over a bound time interval $\mathcal{N}_{\gamma}(T)=\int_{t}^{t^{\prime}}\left(|\gamma(\tau)|-\gamma(\tau)\right)\,d\tau$ Hall ; Rivas . At first glance, we can witness NM from the negative values of $\gamma(t)$ Luo ; Vacchini ; Haikka ; Rivas ; Breuer , which means that the previous discussion about the negative behavior of $\gamma(t)$ in Figure 2 stands as a proof of NM for the orbital states of the SiV- center,— similar conclusions are obtained for the NV- center. This result is the first evidence that the phononic contribution induces NM behavior in color centers in diamond, commonly modeled as purely Markovian Ralf2014 ; Kepesidis2016 ; Ariel2018 .

In Figure 3-(a) we show the temperature dependence of $\mathcal{N}_{C}(T)$ (dashed lines) and $\mathcal{N}_{\gamma}(T)$ (solid lines) for both color centers in diamond (NV*-* and SiV*-*). The system was prepared in the initial pure state $\rho_{0}=|\psi(0)\rangle\langle\psi(0)|$, where $\ket{\psi(0)}=\left(\ket{e}+\ket{g}\right)/\sqrt{2}$, and hereafter we will only focus on this condition. It is interesting that both measures are almost constant at low temperatures, but above $T\approx 100$ K, $\mathcal{N}_{\gamma}(T)$ starts to increase linearly with temperature, while in contrast, $\mathcal{N}_{C}$ goes to zero. The former can be explained by noticing that $\gamma_{loc1}$ increases with temperature, and would lead to the conclusion that the bath becomes more NM with increasing temperature. For the latter, we plotted the time evolution of the Coherence in Figure 3-(b). Here one immediately sees that at high temperature the NM disappears, as one would expect from the coupling to bulk phonons. To shed more light on this matter, the key is to look at the unusual and complex spectral density function of these systems. Onnone hand, at low temperatures the dynamics is ruled out by the SDF $J(\omega)$, since all phonons are frozen out. For this reason, the quasi-localized phonon has the leading contribution to NM, in agreement with the remarks given in Ref. Vasile for an engineered reservoir. On the other hand, at high temperatures, low-frequency (acoustic) phonons are dominant ($\gamma^{\uparrow}(t)\sim\int_{0}^{\infty}J(\omega)/(\hbar\omega)^{2}\sin\omega t\,d\omega$), and therefore the reservoir can be modeled by a super-Ohmic spectrum. In this scenario, NM is highly suppressed Haikka .

To further support the physical consistency of our coherence measure, we compare now $\mathcal{N}_{C}$ with the well-known measure proposed by Breuer, Laine, and Piilo ($\mathcal{N}_{BLP}$) Breuer , where the figure of merit is the trace distance between two quantum states $D(t)=\mbox{tr}|\rho_{1}-\rho_{2}|/2$, (see Appendix D for more details). Both measures have a similar behavior as a function of temperature, see Figure 5, although $\mathcal{N}_{C}$ neither requires a maximization over the set of initial states ($\rho_{1}(0)$ and $\rho_{2}(0)$) nor a full tomography of the density matrix. Moreover, it can be easily computed and experimentally measured through a system observable such as $\langle\sigma_{x}\rangle$. Also, it does not involve an auxiliary system (in contrast to the case of quantum mutual information measure Luo ) and reaches its maximum when $\mathcal{N}_{C}=1$ (supported by the results in Ref. Chruscinski for the decay rate).

V.1 Non-markovianity measure for coherent maps

In the current and many other settings of interest, the Rabi-frequency $\Omega$ in Hamiltonian (2) is controlled by an external drive and can be switched off. However, this might not be the case for other systems, where a non-vanishing coupling can induce oscillations between the states $|g\rangle$ and $|e\rangle$. As already pointed out in previous works Chanda ; Passos , this induces as well oscillations in the Coherence, leading to a false positive when detecting NM. Therefore, we propose here an experimental sequence to specifically filter the dynamics corresponding to the $(\Omega/2)\sigma_{x}$ component. This procedure is first demonstrated for the weak-coupling regime, where we can derive an exact master equation for $\rho_{s}(t)$, as detailed in Appendix E. The resulting time evolution for the off-diagonal elements is given by $\rho_{eg}(t)=\rho_{eg}(0)\mbox{exp}(-F_{0}(t)+F(t)/2)$, where $F_{0}(t)$ reduces to $\Gamma(t)$ for $\Omega=0$. The function $F(t)$ contains contributions, which are solely related to coherent population oscillations, and therefore must be cancelled. This can be done by using the fact that the quantity $S=(1/2)(\langle S_{z}\rangle_{(\rho_{ee}=1)}-\langle S_{z}\rangle_{(\rho_{gg}=1)})=\mbox{exp}(F(t))$ (as given in Eq. (35)) depends on $F(t)$ only. Here the subscripts indicate that the average dynamics of $\langle S_{z}\rangle$ is calculated starting from the initial excited ($e$) and ground ($g$) eigenstates of the Hamiltonian $H_{s}=-(\Delta/2)\sigma_{z}+(\Omega/2)\sigma_{x}$, respectively. The operators $S_{i}$ are the usual Pauli matrices written in the system eigenbasis. In summary, we perform the following three steps: $(1)$ perform a measurement of $\langle S_{z}\rangle$ for the initial condition $\rho_{gg}(0)=1$, $(2)$ repeat the measurement for $\rho_{ee}(0)=1$, $(3)$ calculate the Coherence by performing a measurement of $\langle S_{x}\rangle$ for the initial pure state $|\varphi(0)\rangle=(|1\rangle+|2\rangle)/\sqrt{2}$, and multiply the outcome by $S^{-1/2}$, where $|1\rangle$ and $|2\rangle$ are defined in Appendix E. From this procedure we define a renormalized coherence $\tilde{C}(t)=C(t)\times S^{-1/2}$ which can be used instead of $C(t)$ in Eq. (1). Note that this sequence is unbiased with the spectral density function and it takes advantage of the dynamics induced by the $\sigma_{x}$ component.

For further illustration, we plotted in Figure 3-(c) the coherence calculated with and without the filtering function $S^{-1/2}$ for the SDF $J(\omega)=\eta\omega\exp(-\omega/\omega_{c})$, which is well-know to lead to Markovian dynamics at arbitrary temperaurtures Haikka ; Fanchini2010 . First, because of the presence of $\sigma_{x}$-term in Hamiltonian (2), $C(t)$ shows an oscillation that the measure $\mathcal{N}_{C}=0.14$ detects as non-Markovianity. Secondly, our filter function eliminates the oscillation, and we obtained $\mathcal{N}_{\tilde{C}}=0$ for $\tilde{C}(t)$, as expected for this SDF. Our sequence takes advantage of the weak-coupling between the two-level system and phonons, and gets the information from $\sigma_{x}$ and uses it to cancel its contribution on the off-diagonal terms.

Beyond the weak coupling regime the generalized polaron transformation offers the possibility to apply a similar procedure for strong system-bath interactions Ahsan2011 ; Lukin2013 . This transformation is defined as $H^{\prime}=e^{S_{1}}He^{-S_{1}}$ with $S_{1}=\sigma_{z}\sum_{k}(f_{k}/\omega_{k})(b_{k}^{\dagger}-b_{k})$ and leads to $H^{\prime}=H_{0}+H_{ph}+V$, where $H_{0}=-(\Delta/2)\sigma_{z}+(\Omega_{R}/2)\sigma_{x}$ is the system Hamiltonian (neglecting the polaron shift), $H_{ph}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$ is the phonon Hamiltonian, $V=\sigma_{x}V_{x}+\sigma_{y}V_{y}+\sigma_{z}V_{z}$ is the interaction Hamiltonian (see Eqs. (38)-(40) in Appendix F) and $\Omega_{R}=B\Omega$, where $B$ is the renormalization factor Ahsan2011 ; Lukin2013 , see Appendix F. In Figure 3-(d) we plotted the renormalization factor $B$ for the full polaron transformation (FPT) ($f_{k}=g_{k}$) and the variational polaron transformation (VPT) ($f_{k}=F(\omega_{k})g_{k}$), where $F(\omega_{k})=[1+\Omega_{R}^{2}/(\omega_{k}\omega_{0})\mbox{coth}\left(\beta\omega_{k}/2\right)\mbox{tanh}\left(\beta\omega_{0}/2\right)]^{-1}$. We observed that for not too large values of $\Omega$, the renormalized term $(\Omega_{R}/2)\sigma_{x}$ has a negligible effect in the open dynamics above a temperature $T\approx 120$ K. This temperature is very close to the temperature associated with the quasi-localized phonon mode $T_{loc}=\hbar\omega_{loc}/k_{B}\approx 116$ K. Therefore, oscillations induced by $\sigma_{x}$ on the coherence will not be observed at high temperatures Lukin2013 .

VI Conclusions

In summary, we studied a simple measure to quantify the degree of Non-Markovianity (NM) that takes advantage of the pure dephasing dynamics and measures the back-flow of information from the environment to the system. We compare it with other established measures of NM Breuer ; Hall ; Rivas , and observed that it is well behaved both at low and high temperatures. Also, it is easy to compute and experimentally accessible through Ramsey spectroscopy. Moreover, we studied the dynamics of SiV- and NV- centers due to the vibrations of the diamond lattice and found that the competition of acoustic and quasi-localized phonon modes give rise to a NM dynamics with a rich thermal dependence. As a consequence, different measures of NM exhibit opposite dependence at high temperature, suggesting that the power of each measure will rely on how it is linked to a specific application. Furthermore, we extended this NM measure from incoherent (independent spin-boson model) to coherent dynamical maps (general spin-boson model), where for the weak coupling we found that following a specific experimental sequence one can filter out the undesired contribution of the $\sigma_{x}$-term, while for the strong coupling regime there is a temperature $T\approx 120$ K above which this term is negligible.

VII acknowledgments

We are grateful to Felipe Fanchini for fruitful discussions about non-Markovianity. AN acknowledges financial support from Universidad Mayor through the Postdoctoral fellowship. JRM acknowledges support from Fondecyt Regular No. 1180673 and AFOSR FA9550-18-1-0513. PR acknowledges support from the Austrian Science Fund (FWF) through Grant No. P 31701-N27. RC acknowledge financial support from Fondecyt Iniciación No. 11180143.

Appendix A Exact reduced dynamics for the independent spin-boson model

In this section we introduce the exact master equation associated with the reduced two-level system (TLS) of the independent spin-boson Hamiltonian given by

[TABLE]

where $\omega_{0}$ is the bare frequency of the TLS, $\omega_{k}$ are the phonons frequencies and $g_{k}$ are the electron-phonon coupling constants. First, we consider the Born approximation for the initial state $\rho(0)=\rho_{s}(0)\otimes\rho_{ph}$, where $\rho_{ph}=e^{-\beta H_{ph}}/\mbox{Tr}\{e^{-\beta H_{ph}}\}$ is a thermal bosonic state, with $\beta=(k_{B}T)^{-1}$, $H_{ph}=\sum_{k}\omega_{k}a_{k}^{\dagger}a_{k}$ and $\rho_{s}(t)=\mbox{Tr}_{ph}(\rho(t))$. In the Heisenberg picture, the exact equations of motion for the expectation values $\langle\sigma_{\pm}(t)\rangle$ are given by Luczka

[TABLE]

where $\sigma_{\pm}=(\sigma_{x}\pm i\sigma_{y})/2$, $\langle\sigma_{\pm}(t)\rangle=\mbox{Tr}_{s}\left(\sigma_{\pm}(t)\rho_{s}(0)\right)$, $\gamma(t)$ is the dephasing rate introduced in Eq. (4), $T$ is the reservoir temperature, $k_{B}$ is the Boltzmann constant, and $J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$ is the spectral density function of the reservoir. In general, the density matrix of a two-level system can be written as $\rho_{s}(t)=(1/2)\left(\mathds{1}+\langle\vec{\sigma}(t)\rangle\cdot\vec{\sigma}\right)$, where $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$, $\mathds{1}$ is the identity matrix, and $\sigma_{i}$ are the usual Pauli matrices. Using the last decomposition and Eq. (12), we derived the following equation

[TABLE]

where we have used $\sigma_{z}(t)=\sigma_{z}(0)$ (constant of motion) due to the fact that $[H_{s},\sigma_{z}]=0$ for $H_{s}=(\omega_{0}/2)\sigma_{z}$. Using the property $\sigma_{i}\sigma_{j}=\delta_{ij}\mathds{1}+i\epsilon_{ijk}\sigma_{k}$, we obtained the relation

[TABLE]

From the above relation and $\rho_{s}(t)=(1/2)\left(\mathds{1}+\langle\vec{\sigma}(t)\rangle\cdot\vec{\sigma}\right)$ we get $\langle\sigma_{x}(t)\rangle\sigma_{x}+\langle\sigma_{y}(t)\rangle\sigma_{y}=\rho_{s}(t)-\sigma_{z}\rho_{s}(t)\sigma_{z}$. Therefore, the master equation for $\rho_{s}(t)$ is given by

[TABLE]

By neglecting the first term we recover the dynamics in the interaction picture (see Eq. (2) of the main text).

Appendix B Dephasing rate induced by strong interactions with quasi-localized phonons

The dephasing rate associated with the strong interaction with a quasi-localized phonon mode is given by

[TABLE]

This integral can be solved analytically, however, we show next a method to obtain a good approximation that gives us a better understanding of the effect of the width $\Gamma$, the frequency of the quasi-localized phonon $\omega_{loc}$, the amplitude $J_{0}$, and temperature $T$. Using the change of variable $u=\omega-\omega_{loc}$ and extending the lower limit of the integration to $-\infty$ (assuming $\omega_{loc}\gg 1$), we obtain

[TABLE]

where

[TABLE]

The main contribution in both integrals given in Eq. (17) comes from the narrow Lorentzian function $L(u)=(\Gamma/2)/(u^{2}+(\Gamma/2)^{2})$ around the value $u=0$ (main peak around $\omega=\omega_{loc}$ for $J(\omega)$, see Figure 1). Using the approximation

[TABLE]

and the symmetry consideration $\int_{-\infty}^{\infty}f(u)\sin(ut)\,du=0$, we obtain

[TABLE]

which corresponds to a damped periodic oscillation, where $\omega_{loc}$ is the frequency, $\Gamma/2$ is the decay rate and both temperature $T$ and strength $J_{0}$ determine the maximum amplitude of the oscillations. At zero temperature, $\coth(\hbar\omega/2k_{B}T)=1$, and therefore, we recover the expression given in Eq. (9). At high temperatures, $k_{B}T\gg\hbar\omega_{loc}$, we have $\coth(\hbar\omega/2k_{B}T)\approx 2k_{B}T/\hbar\omega$, and then

[TABLE]

Appendix C Dephasing Rate for an NV- center

The dephasing rate $\gamma_{NV}(t)$ is calculated from Eq. (4) with the spectral density function given numerically in Ref. Alkauskas , and illustrated in Figure 1.

Appendix D Measures of the degree of non-Markovianity

For comparison, we introduce now the Breuer-Laine-Piilo (BLP) measure Breuer , which is based on the trace distance

[TABLE]

where $|A|=\sqrt{A^{\dagger}A}$. The BLP measure is plotted in Figure. 5 (dotted line) for the initial states

[TABLE]

In addition, we noted that the Rivas-Huelga-Plenio (RHP) measure is given by the relation $\mathcal{N}_{\gamma}=(d/2)\mathcal{N}_{RHP}$ Rivas , where $d=2$ is the dimension of the Hilbert space in our case (two-level system), and therefore $\mathcal{N}_{\gamma}$ and $\mathcal{N}_{RHP}$ are identical.

Notice that $\mathcal{N}_{\gamma}$ (solid line) shows a thermal dependence that is based on the dephasing rate $\gamma(t)$ given in Eq. (4), and it has a strictly increasing behavior leading to a large degree of NM when the temperature is above some critical value $T\sim 30$ K. However, we observed that $\mathcal{N}_{BLP}$ shows a strictly decreasing behavior leading to a small degree of NM at room temperature. Finally, $\mathcal{N}_{C}$ (dashed-line) shows a thermal dependence that is similar to $\mathcal{N}_{BLP}$, as one would expect since both NM measures are based on the back-flow of quantum information.

Appendix E Weak coupling model

For a better understanding of the dynamics of the system when it is weakly coupled to the environment, we traced over the reservoir in Eq. (2), assuming that the reservoir is in a thermal state. After the Born and the secular approximations, one gets the following master equation Laine ,

[TABLE]

where $H_{s}=(\omega_{0}/2)S_{z}$, $\omega_{0}=\sqrt{\Delta^{2}+\Omega^{2}}$, and $S_{z}=|1\rangle\langle 1|-|2\rangle\langle 2|$, being

[TABLE]

the system eigenstates. Notice that the Lamb shift Hamiltonian has been omitted since it is negligible according our numerical simulations. Furthermore, the Lindblandian is given by,

[TABLE]

Here, $L_{0}=(\Delta/\omega_{0})S_{z}$, $L_{+\omega_{0}}=(\Omega/\omega_{0})S_{-}$ and $L_{-\omega_{0}}=(\Omega/\omega_{0})S_{+}$, where $S_{+}=|1\rangle\langle 2|$ and $S_{-}=|2\rangle\langle 1|$ are the raising and lowering operators in the system eigenbasis. Moreover, for arbitrary temperatures, we have,

[TABLE]

where $n(\omega)=[\mbox{exp}(\hbar\omega/k_{B}T)-1]^{-1}$ is the mean number of phonons at thermal equilibrium and $J(\omega)$ is the spectral density function introduced in Eq. (4). One can immediately notice that for $\Omega=0$, we recover the dynamics in Eq. (3). Since the Hamiltonian evolution $\sim[S_{z},\rho_{s}(t)]$ does not affect the coherence, we will only focus on the dynamics induced by $D\rho_{s}(t)$. The time evolution of the density operator is given by,

[TABLE]

where $\rho_{ee}(t)=\langle 1|\rho(t)|1\rangle$, $\rho_{gg}(t)=\langle 2|\rho(t)|2\rangle$ and $\rho_{eg}(t)=\langle 1|\rho(t)|2\rangle$. The above equations can be formally integrated,

[TABLE]

where $F(t)=-(\Omega/\omega_{0})^{2}\int_{0}^{t}(\gamma_{-\omega_{0}}(t_{1})+\gamma_{\omega_{0}}(t_{1}))\,dt_{1}$, $g(t)=(\Omega/\omega_{0})^{2}\gamma_{-\omega_{0}}(t)$, and $F_{0}(t)=2(\Delta/\omega_{0})^{2}\int_{0}^{t}\gamma_{0}(t_{1})\,dt_{1}$. Notice that $F_{0}(t)$ corresponds to the dynamics induced by $S_{z}$, that does not generate coherence. Only the contribution of $F(t)$ in Eq. (32) will lead to a false non-Markovianity, and therefore we will cancel it. To do so, we first calculate the expectation values of $S_{z}$ for different initial conditions.

[TABLE]

From these results it is straight forward to obtain that

[TABLE]

Furthermore, we have that $\langle S_{x}\rangle=2\rho_{eg}$, given that $\rho_{eg}\in\mathbb{R}$ and $S_{x}=S_{+}+S_{-}$. Hence, for our generalize measure of NM, the Coherence will be defined as $\tilde{C}(t)=\langle S_{x}\rangle\times S^{-1/2}=2\rho_{eg}(0)e^{-F_{0}(t)}$, where for an incoherent dynamics $S=1$. Finally, one can observe that this expression for the Coherence has no contribution from $(\Omega/2)\sigma_{x}$ in Eq. (2) (as opposed to the case where $C(t)$ is calculated directly from Eq. (32) and even more, we have not done any assumption about the spectral density function $J(\omega)$, which makes this approach quite general.

Appendix F Strong coupling model

In order to study strong interactions between a two-level system and its phononic environment we introduce the following general polaron transformation Lukin2013 ; Ahsan2011 ($\hbar=1$)

[TABLE]

where $f_{k}=g_{k}$ corresponds to the the full polaron transformation (FPT). If we apply the above transformation on the general spin-boson Hamiltonian $H=-(\Delta/2)\sigma_{z}+(\Omega/2)\sigma_{x}+\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}+\sigma_{z}\sum_{k}g_{k}(b_{k}^{\dagger}+b_{k})$, we obtain

[TABLE]

The interaction Hamiltonian is given by $V=\sigma_{x}V_{x}+\sigma_{y}V_{y}+\sigma_{z}V_{z}$, where

[TABLE]

The bath operators are defined as

[TABLE]

with $\beta=(k_{B}T)^{-1}$ and the expectation value is calculated using the thermal phonon state $\rho_{ph}=\mbox{exp}(-\beta H_{ph})/\mbox{Tr}\{\mbox{exp}(-\beta H_{ph})\}$. To determine the optimal values of $f_{k}$ for the variational polaron transformation (VPT) it is necessary to minimize the free energy

[TABLE]

where $H_{0}=H_{A}+H_{B}$ is the non-interacting Hamiltonian with $H_{A}=-(\Delta/2)\sigma_{z}+(\Omega_{R}/2)\sigma_{x}+\sum_{k}(f_{k}/\omega_{k})\left(f_{k}-2g_{k}\right)$ and $H_{B}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$. Using $\langle V\rangle_{H_{0}}=0$ and the condition $dA_{B}/df_{k}=0$ we obtain that $f_{k}=g_{k}F(\omega_{k})$, where

[TABLE]

and $\omega_{0}=\sqrt{\Delta^{2}+\Omega_{R}^{2}}$. In the continuum limit, the renormalization factor $B$ is given by

[TABLE]

This renormalization factor depends on the shape of the SDF $J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$, the reservoir temperature $T$, and must be calculated from self-consistency between Eqs. (44)-(46).

Appendix G Role of different spectral density function on the Coherence

From Figure 3-(b), it is easy to see that Coherence $C(t)$ obeys two different regimes separated in temperature. At low temperature, the time evolution of $C(t)$ has contributions from of all the spectral density functions in Eqs. (5)-(7), however the strong oscillation in the main plot evidences that $\gamma_{loc1}$ is the leading contribution. Even when one would expect $\gamma_{loc1}$ to hold as the leading term in the dynamics at high temperatures, see Figure 2, this is not the case. In Figure 6 we showed $C(t)$ for two particular decay rates, namely $\gamma_{bulk}(t)$ (solid) and $\gamma_{loc1}(t)$ (dashed), at $300$ K. It is remarkable that $\gamma_{bulk}(t)$ reproduces the behavior (in terms of Non-Markovianity) depicted in the inset of Figure 3-(b), despite that it decays slower due to the absence of the other decay rates. This outcome supports the statement that only bulk phonons are relevant at high temperature Goldman2015 ; Goldman2017 . In contrast, $\gamma_{loc1}(t)$ shows a very interesting dynamics as well.

In particular, similar behavior have been observed in spin-echo spectroscopy for a single NV- Childress or an ensemble Stanwix interacting with a natural environment of ${}^{13}\mbox{C}$ nuclear spins, where the collapses and revivals originate from a coherent interaction with individual proximal nuclear spins, and the revivals decay comes from the interaction with the ${}^{13}\mbox{C}$ bath. In the same way, this behavior has been observed for a single NV interacting with a mechanical oscillator Bennet . In agreement with these results, we suggest that the interaction with this quasi-localized phonon mode can be though as an interaction with a single mode in a coherent state, with an added phenomenological decay that goes as $\exp(-\Gamma t)$, being $\Gamma=0.8414$ THz the width of the Lorentzian spectral density function. The composite state after this interaction is given by Montenegro ,

[TABLE]

where $\xi=1-e^{-it\omega_{loc}}$, $|\beta_{\uparrow}\rangle=|\beta e^{-it\omega_{loc}}+\lambda\xi\rangle$ ($|\beta_{\downarrow}\rangle=|\beta e^{-it\omega_{loc}}-\lambda\xi\rangle$) is the displaced coherent state, $\lambda=(\int_{0}^{\infty}J_{loc1}(\omega))^{1/2}\,d\omega$. Even more, we approximated $|\beta|^{2}$ to the thermal occupancy phonon number $n(\omega_{loc})=[\exp(\hbar\omega_{loc}/k_{b}T)-1]^{-1}$, given that at $T=300$ K the occupancy is small ($n(\omega_{loc})\approx 2$). Finally, after tracing out over the coherent state degrees of freedom we are able to calculate $C(t)$. This highly simplified model is capable to capture the collapses, revivals, and the overall decay, as illustrate by the dot-dashed curve.