Coherent backaction between spins and an electronic bath: Non-Markovian dynamics and low-temperature quantum thermodynamic electron cooling

TL;DR

This paper develops an analytical framework to model non-Markovian spin dynamics in contact with an electronic bath, revealing full decay behavior and proposing a quantum thermodynamic cooling scheme for electrons at low temperatures.

Contribution

It introduces a systematic method to analyze non-Markovian spin-bath interactions, extending beyond traditional Markovian approximations, and proposes a novel electron cooling technique.

Findings

Full analytic description of spin decay from non-Markovian to exponential regimes

Demonstration of temperature-insensitive decay useful for cooling

Proposal of a spin-based scheme to cool electrons at very low temperatures

Abstract

We provide a versatile analytical framework for calculating the dynamics of a spin system in contact with a fermionic bath beyond the Markov approximation. The approach is based on a second order expansion of the Nakajima-Zwanzig master equation but systematically includes all quantum coherent memory effects leading to non-Markovian dynamics. Our results describe, for the free induction decay, the full time range from the non-Markovian dynamics at short times, to the well-known exponential thermal decay at long times. We provide full analytic results for the entire time range using a bath of itinerant electrons as an archetype for universal quantum fluctuations. Furthermore, we propose a quantum thermodynamic scheme to employ the temperature insensitivity of the non-Markovian decay to transport heat out of the electron system and thus, by repeated reinitialization of a cluster of spins, to efficiently cool the electrons at very low temperatures.

Equations

$$H=\sum _{k\sigma }{ϵ}_k{c}_{k\sigma }^{†}{c}_{k\sigma }+b_z^{\text{el}}\sum _j{S}_{z,j}+b_z^I{I}_z+A\mathbf{I}\cdot {\mathbf{S}}_{j=0}.$$

$$I_{\pm }$$

$$I_z$$

$$I_{\pm }=\frac{1}{2}\left(I_x\pm i{I}_y\right),S_{\pm ,j}=\frac{1}{2}\left(S_{x,j}\pm i{S}_{y,j}\right).$$

$$I_{\uparrow }=\frac{1}{2}(1+I_z),I_{\downarrow }=\frac{1}{2}(1-I_z),$$

$$H$$

$${\rho }_I={\text{Tr}}_{\text{el}}\left[\rho \right],$$

$$\frac{d}{dt}{\rho }_I(t)=-i{\int }_0^{t}d{t}^{\prime }{\Sigma }_I\left(t-t^{\prime }\right){\rho }_I(t^{\prime }).$$

$${\Sigma }_I(t-t^{\prime })=-i{\text{Tr}}_{\text{el}}\left[L{e}^{-iQL(t-t^{\prime })}Q{L}_{\text{int}}{\rho }_{\text{el}}\right],$$

$${\tilde{\rho }}_I(s)$$

$${\tilde{\Sigma }}_I(s)=-i{\text{Tr}}_{\text{el}}\left[L{\left(s1+iQL\right)}^{-1}Q{L}_{\text{int}}{\rho }_{\text{el}}\right].$$

$${\tilde{\Sigma }}_I^{\text{Born}}(s)=-i{\text{Tr}}_{\text{el}}\left[L_{\text{int}}{\left(s1+Q{L}_0\right)}^{-1}Q{L}_{\text{int}}{\rho }_{\text{el}}\right].$$

$${\rho }_I={\rho }_{\uparrow }{I}_{\uparrow }+{\rho }_{\downarrow }{I}_{\downarrow }+{\rho }_{-}{I}_{-}+{\rho }_{+}{I}_{+},$$

$$\left[{\tilde{\Sigma }}_I^{\text{Born}}(s)\right]=\left(\begin{array}{cccc}{F}_1& -F_2& 0& 0\\ -F_1& F_2& 0& 0\\ 0& 0& F_{-}+F_z& 0\\ 0& 0& 0& F_{+}+F_z\end{array}\right).$$

$${\rho }_z(t)$$

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

$$F_1(s)$$

$$F_z(s)=\frac{A^2}{2i}{\int }_0^{\infty }dt{e}^{-ts}⟨\left\{S_{z,j=0}(t),S_{z,j=0}(0)\right\}⟩.$$

$$\alpha ={\left({\nu }_{0}A\right)}^2.$$

$$\displaystyle F_{z}(s)=i\alpha\bigl{[}F(s)+G(s)\bigr{]},$$

$$F(s)$$

$$+s\left[\ln \left(\frac{2\pi k_{B}T}{i\xi }\right)+\psi \left(1+\frac{s}{2\pi k_{B}T}\right)\right],$$

$$G(s)$$

$$+s\left[\ln \left(\frac{2\pi k_{B}T}{-i\xi }\right)+\psi \left(1+\frac{s}{2\pi k_{B}T}\right)\right],$$

$$F_z(s)\sim 2i\alpha s\ln \left(\frac{s}{\xi }\right).$$

$$F_z(s)\sim -2\alpha i\pi k_{B}T+2i\alpha s\ln \left(\frac{2\pi k_{B}T}{{\xi }^{\prime }}\right).$$

$$F_1(s)$$

$$F_2(s)$$

$$F_{-}(s)$$

$$F_{+}(s)$$

$$F_z(s)$$

$$T_1$$

$$T_2$$

$$T_1=T_2=\frac{\hslash }{4\alpha \pi k_{B}T},$$

$${\rho }_z^{\text{M}}(t)={\rho }_z^{\text{eq}}$$

Full text

I Introduction

In recent years, quantum effects have become controllable, to the extent that a proper quantum technology based on the core properties of quantum parallelism and entanglement is on the rise. The success of this technology will rely on our ability to understand and control more intricate quantum states and their evolution. One of the main challenges is to control the interaction of a quantum system with its environment Kloeffel et al. (2011); Yang and Sham (2013); Jun et al. (2016). Such an interaction is detrimental for a quantum application as it leads to the major issue of decoherence, but it also can be used positively, as it is on the basis of quantum thermodynamics Kosloff (2013); Goold et al. (2016); Vinjanampathy and Anders (2016), because it can be employed to manipulate the system via the environment Stace et al. (2005); Härtle et al. (2018) and to reveal information about the environment itself. The latter is, for instance, at the core of magnetic resonance techniques such as nuclear magnetic resonance and electron paramagnetic Resonance in which the measurable decay of individual quantum spins is the result of their interaction with their environment Ihara et al. (2010); Mamin et al. (2007); Augustyniak-Jabłokow et al. (2013); Hanson et al. (2006). Figure 1 shows a sketch of such a system.

It is the goal of this paper to provide an analytical framework to access the system-environment correlations, and thus to provide a direct access to the core physics that connects the broad range of physical disciplines from the mature field of magnetic resonance to the recent development of quantum thermodynamics. Although we restrict ourselves to the weak coupling limit, we provide a systematic approach to include all Markovian and non-Markovian dynamics. This allows us to bridge fully from the quantum coherent regime at short time scales to the thermodynamic regime at long time scales.

In common situations the number of degrees of freedom of the environment is macroscopic. Then, although the overall evolution remains unitary, the environment acts as a bath in which any transmitted information is effectively dispersed irreversibly. The environment is called Markovian or memoryless if during the contact with the system it remains unaffected by the interaction in its equilibrium state. Since a zero memory time cannot exist physically, the Markovian property is an approximation that is valid if the memory time is shorter than the characteristic time scales of the dynamics of the quantum system. While this approximation sufficiently describes quantum thermodynamics, decoherence or the behavior of magnetic resonance and similar techniques in many situations, it misses that during the memory time, even if it is short, the system and environment have a joint evolution in which in particular the excitation of the environment can act back on the system. This partially coherent backaction not only shapes the system’s short time dynamics but can also leave an imprint in the long time behavior in the form of a correction to the expected Markovian dynamics Hell et al. (2016); DiVincenzo and Loss (2005); Loss and DiVincenzo (1998) or a modification of the performance in thermal machines Watanabe et al. (2017). Such properties can thus be used passively as a diagnostic tool for the structure of the environment. But such non-Markovian behavior can also be used as a novel route to actively manipulate the environment through the quantum system Man et al. (2015) or to enhance the efficiency of quantum heat engines Abiuso and Giovannetti (2019).

The non-Markovian treatment becomes pertinent when memory times are enhanced, such as the modern developments of reaching sub-millikelvin temperatures even for electronic nanostructures Clark et al. (2010); Palma et al. (2017a, b) and by the design and discovery of strongly correlated materials. The latter exhibit collective responses to local excitations that naturally extend correlations in space and time, and thus the memory time. Within the Markov approximation, only the spatial signatures of strong correlations have an impact on the system’s dynamics, which is a property routinely investigated for signatures of the strong correlation physics Kölbl et al. (2012); Türeci et al. (2011); Latta et al. (2013). This neglects that the similarly strong temporal correlations delay the decay of memory and that the dynamics must be complemented by a concise modeling of the non-Markovian dynamics of, for instance, a spin’s free decay in magnetic resonance.

When evolving a system under a Markovian master equation the bath is treated as a large memoryless reservoir. The state of a system in contact with such a reservoir thermalizes and decays exponentially fast. In the non-Markovian regime, the backaction onto the system by fluctuations in the bath is important as the times are too short for the memory kernel to decay sufficiently.

As an exact solution can rarely be found, different approaches exist to incorporate the short time dynamics and non-Markovian effects based, for instance, on time-convolutionless master equations Timm (2011); Fischer and Breuer (2013); Jing and Wu (2018); Maniscalco and Petruccione (2006) or the Nakajima-Zwanzig equation Nakajima (1958); Zwanzig (1960). The Nakajima-Zwanzig equation is often used as starting point to set up hierarchical equation of motions Ishizaki and Tanimura (2005); Zhang and Yan (2016) or to modified projection methods DiVincenzo and Loss (2005); Degenfeld-Schonburg and Hartmann (2014) including nonperturbative expansions Coish and Loss (2004); Moustos and Anastopoulos (2017) and re-summation techniques Chen et al. (2016). Numerical algorithms include transfer tensors or path integral methods to reconstruct the memory kernel of the Nakajima-Zwanzig equation, which are even feasible beyond the weak coupling limit Pollock et al. (2018); Strathearn et al. (2018); Cerrillo and Cao (2014); Shi and Geva (2003). Another possibility is to calculate the ‘initial slip’ of the system and then evolving the state within a Markovian description with modified initial conditions Geigenmüller et al. (1983); Haake and Reibold (1985); Suárez et al. (1992); Yu et al. (2000). It should just be emphasised that any chosen approach has to maintain the positivity of the density matrix Loss and DiVincenzo (1998); Whitney (2008); Vacchini (2016); Shabani and Lidar (2005). Closest in spirit to our approach is perhaps a Liouville-space decomposition method Hell et al. (2016, 2014) with its emphasis on the importance of coherent backaction.

In this paper, we provide a theoretical framework for the entire cross-over from short to long times for a spin system coupled to a fermionic system, starting from the exact Nakajima-Zwanzig equation. We pursue three goals: (i) the development of the general formalism applicable to any spin and any environment system. (ii) the analysis of non-Markovian behavior when the environment is a Fermi gas, and (iii) the demonstration that non-Markovianity can be used to actively manipulate the environment through the proposition of a temperature independent quantum demagnetisation cooling protocol that goes beyond a standard thermodynamic cycle.

For goal (i), we set up a systematic analysis of the poles representing the various decay modes from the equation of motion of the density matrix in Laplace space. The important feature is that we identify the ensemble of decay modes but perform approximations to the exact pole positions. This means that although we will use a second order expansion of the memory kernel we pick up the entity of decay modes and can cover with them the full time dynamics through all times. This approach is versatile and independent of the type of the environment, but to give a concrete example we focus on fermionic environments for the further study of non-Markovian behavior. Although we expect the largest impact to arise from strongly correlated systems, we have chosen for goal (ii) an environment in the form of a noninteracting itinerant fermionic system, sketched in Fig. 1. This environment is as close to a memoryless fermionic bath as possible. Nevertheless, we demonstrate that it exhibits non-Markovian effects that have a considerable impact on both the short- and long time dynamics, and thus even for long times it has an effect that cannot be captured by a standard thermodynamic description. In addition, the simple fermionic gas shares the fact with strongly correlated (critical) systems that its properties depend only on a few global parameters, in this case the Fermi energy $E_{F}$ (or the band width $\xi$, assuming that $\xi\sim E_{F}$) and the temperature $T$. The obtained results therefore provide also a hint as to what could be expected as signatures of non-Markovian behavior of correlated systems. In a similar way a bosonic environment could be investigated. Goal (iii) involves a proof of principle that the non-Markovianity allows for an active manipulation of the environment. Since the non-Markovianity is driven by quantum fluctuations we propose to use it as a largely temperature independent way to transport heat out of the environment. Such a method would thus make it possible to overcome the bottleneck of diverging time scales in demagnetisation cooling methods which arises because the used relaxation times scale with $1/T$.

The paper is structured as follows. In Section II we introduce the generalised master equation within the projection operator approach and its formal solution for a particular system. This approach is a controlled expansion in the interaction strength, capturing the full coherent quantum dynamics for non-Markovian regime as well as the long time evolution. Taking into account temporal correlations in the bath we investigate how the dynamics of the spin system is affected if a backaction from the bath is present. In Section III we derive an analytical solution for the dynamics of the spin system, including the Markovian and non-Markovian decay. In Section IV we propose a cooling protocol to utilize the short time dynamics to overcome thermodynamic limitations in the context of cooling an electronic system.

II Tracking the full time evolution of a quantum system

Our aim is to determine the time evolution of a spin system coupled to a fermionic environment including the full quantum coherent dynamics. In particular, our work captures the influence of memory effects on the dynamics. The generic system considered here consists of an impurity spin-1/2, e.g. a nuclear spin or a localized paramagnetic electron spin, coupled to an electronic bath. It is described by the Kondo-type Hamiltonian

[TABLE]

The first term models the electronic environment with annihilation and creation operators $c_{k\sigma}$ and $c^{\dagger}_{k\sigma}$. The momenta $k$ label all possible states with energy dispersion $\epsilon_{k}$, and $\sigma$ is the spin index. The second and third term are the Zeeman terms for the electrons and the impurity spin. Here, $\mathbf{I}$ is the impurity spin operator, and $\mathbf{S}_{j}=\sum_{\sigma,\sigma^{\prime}}c_{j,\sigma}^{\dagger}\bm{\tau}_{\sigma,\sigma^{\prime}}c_{j,\sigma^{\prime}}$ are the electron spin operators, written in terms of the real space operators $c_{j,\sigma}$ and the Pauli matrix vector $\bm{\tau}$. The index $j$ labels the position of the electrons, where for the ease of notation we assume an underlying lattice, but this is not essential for the physics discussed here. For convenience we will normalise both spin operators to dimensionless $|\mathbf{S}_{j}|=|\mathbf{I}|=1$. We will also assume that $\mathbf{I}$ is a spin-$1/2$ operator since this will make the matrix structure of the formalism used below simpler, but this assumption is of no further importance, and the formalism and the results can easily be extended to larger spins. Both spin species interact through a contact interaction with strength $A$. A uniform magnetic field $B_{z}$ is applied along the $z$ axis and we define $b_{z}^{\text{el}}=g\mu_{B}B_{z}$ and $b_{z}^{I}=g_{I}\mu_{I}B_{z}$, with the g-factors $g$ for the electrons and $g_{I}$ for the impurity spins, $\mu_{B}$ for the Bohr magneton, and $\mu_{I}$ for the magnetic moment of the impurity spin.

Although the Hamiltonian in Eq. 1 is of the Kondo-type we emphasize that we do not consider it in the Kondo regime. We focus on small ratios of $A/E_{F}$, where $E_{F}$ is the Fermi energy, such that the Kondo temperature, $T_{K}\sim\exp(-E_{F}/A)$, is much below any realistic situation. The small ratio $A/E_{F}$ also allows us to expand the equation of motion to second order with which we thus neglect the higher order spin-flip terms driving the Kondo effect. Furthermore, we explicitly control definite initial values for the impurity spin and focus on the short time dynamics that would precede the development of Kondo correlations in any case. Nonetheless the elementary particle-hole fluctuations underlying the Kondo effect remain active and cause the logarithmic time dependences found below.

In this work we assume that the electron system remains nonmagnetic, $\langle S_{z}(t)\rangle=0$. This allows us to simplify the equations in the following while keeping all the relevant physical aspects. However, we then also explicitly exclude the Knight shift which describes the shift of the resonances of $\mathbf{I}$ by the effective magnetic field created by the paramagnetic polarization $\langle\mathbf{S}\rangle$ Abragam (1961); Slichter (1990). Although the latter is notable for embedded spins in many electron conductors our focus is on the irreversible behavior from quantum and thermal fluctuations, and the Knight shift can effectively be absorbed by shifting the Zeeman field acting on $\mathbf{I}$. Therefore, to maintain clarity in the discussion we prefer to neglect this shift. Nonetheless we shall in the discussion occasionally come back to its influence.

Changing into a rotating frame of reference eliminates the impurity Zeeman term, $H\rightarrow H-b_{z}^{I}J_{z}$, with the total angular momentum $J_{z}=I_{z}+\sum_{j}S_{z,j}$. The spin operators transform as

[TABLE]

where as usual

[TABLE]

Later we shall also use the notation

[TABLE]

with $\mathbb{1}$ being the identity operator. The Hamiltonian Eq. 1 in the new reference frame can then be written as

[TABLE]

with $b_{z}=b_{z}^{\text{el}}-b_{z}^{I}$.

The full dynamical behavior of the system is encoded in the time evolution of the reduced density matrix of the system

[TABLE]

where $\rho$ is the full density matrix of the system and the environment. We furthermore define the projection operator $P$ such that $PO=\rho_{\text{el}}\otimes\text{Tr}_{\text{el}}[O]$ for any operator $O$, and $Q=\mathbb{1}-P$ as the complement of $P$. The equilibrium density matrix of the electronic system is denoted as $\rho_{\text{el}}$ in its initial state and $\text{Tr}_{\text{el}}$ is the trace over the electronic degrees of freedom. The Nakajima-Zwanzig equation provides a framework to formulate an equation of motion for the reduced density matrix $\rho_{I}(t)$ Breuer and Petruccione (2007); Fick and Sauermann (1990)

[TABLE]

In the latter equation the reduced memory kernel

[TABLE]

carries the information of the system’s history and its time dependence needs to be carefully taken into account to capture the non-Markovian behavior. We split the Hamiltonian $H=H_{0}+H_{\text{int}}$ into two parts. The noninteracting part $H_{0}=~{}\sum_{k,\sigma}\epsilon_{k}c^{\dagger}_{k\sigma}c_{k\sigma}+b_{z}\sum_{j}S_{z,j}$ describes the fermionic bath, and $H_{\text{int}}=A\mathbf{I}\cdot\mathbf{S}_{j=0}$ is the interaction with the spin system, according to Eq. 5. Based on these Hamiltonians we define in the Nakajima-Zwanzig equation the Liouvillian superoperators $L_{0}O=[H_{0},O]$ and $L_{\text{int}}O=[H_{\text{int}},O]$, where $O$ is any operator, and $L=L_{0}+L_{\text{int}}$.

To solve the integro-differential equation of motion, we analyze Eq. 7 in Laplace space. If $f(t)$ is any function of time, the Laplace transform is given by $\tilde{f}(s)=\int_{0}^{\infty}dt\,\exp{(-ts)}f(t)$, for $\text{Re}(s)>0$, and its application on Eq. 7 leads to

[TABLE]

with the initial condition $\rho_{I}(t=0)$ and

[TABLE]

Up to this point, the Nakajima-Zwanzig approach is exact for any Hamiltonian that can be split into interacting and noninteracting parts. A detailed derivation can also be found in Appendix A. To make progress we shall assume that $A<E_{F}$ such that we can expand the memory kernel $\Sigma_{I}$ in Laplace space. We stress that this expansion is only perturbative in the system-bath coupling, but nonperturbative in the time domain. As we will see below we use the expansion of the memory kernel only to facilitate the identification of the decay modes of the equation of motion in Laplace space, but we then will consider the contribution of the entire ensemble of modes to the time evolution. With this we go far beyond the Markov approximation or any expansion in the memory time. Instead the approach shares some aspects of introducing an approximate self-energy term in the Dyson equation for many-body Green’s functions, and indeed the memory kernel $\Sigma$ can be seen as equivalent to the self-energy of nonequilibrium many-body theory. Therefore, it preserves the possibility of the superposition of the infinite number of quantum fluctuations which govern the non-Markovian behavior. The latter is encoded in the pole structure of the correlators in Laplace space and analyzing those presents a systematic way of fully including the non-Markovian contribution which will be discussed in detail later on. It is worth noting that the formalism is independent of the nature of the interaction or the actual structure of the bath and can be applied equally to fermionic or bosonic environments.

The approach provides a controlled expansion of the memory kernel in the interaction while keeping all the information about the system’s past. Within the Born approximation the correction of the memory kernel is quadratic in the interaction term $L_{\text{int}}$, and is thus quadratic in the coupling strength $A$

[TABLE]

Choosing the impurity spin as a basis, i.e. the space spanned by the spin operators $\{I_{\uparrow},I_{\downarrow},I_{-},I_{+}\}$, allows us to decompose any operator into a part that only acts on the bath and a part that only acts on the system. Within this basis the reduced density matrix decomposes as

[TABLE]

such that $\langle I_{\beta}\rangle=\text{Tr}_{I}\left[I_{\beta}\rho_{I}\right]=I\rho_{\beta}$, where $\beta=\uparrow,\downarrow,-,+$ and $\text{Tr}_{I}$ is the trace over the system’s degrees of freedom.

The superoperators can then be represented by $4~{}\times~{}4$ matrices acting on the reduced density matrix vector $(\rho_{\uparrow},\rho_{\downarrow},\rho_{-},\rho_{+})^{T}$ Coish and Loss (2004). In the following, we use square brackets $[O]$ to denote this matrix representation of a superoperator $O$. The memory kernel $\left[\tilde{\Sigma}_{I}^{\text{Born}}(s)\right]$ in its matrix representation takes the form

[TABLE]

The entries of the memory kernel are the Laplace transforms of spin-spin correlation functions $F_{j}(s)$. The derivation of Eq. 13 including the spin-spin correlators $F_{j}$ can be found in Appendix B.

Since the Hamiltonian $H_{0}$ describing the electronic environment is spin conserving the matrix entries in Eq. 13 that do no conserve the electron spin are zero, and the remaining terms $F_{1,2,\pm,z}$ describe variants of $\langle S_{\mp}S_{\pm}\rangle$ or $\langle S_{z}S_{z}\rangle$ electron spin correlators (see Appendix B). The memory kernel $\tilde{\Sigma}_{I}$ takes into account electronic quantum fluctuations induced by an excitation of the impurity spin. Because of the Pauli principle such fluctuations are dominated by particle-hole excitations at the Fermi level. These propagate through the bath and act back onto the impurity spin, creating an effective time-retarded coupling of the impurity spin with itself. Thus, before reaching the thermodynamic equilibrium regime the impurity spin is correlated with its initial state.

With the memory kernel in the matrix representation of Eq. 13 the solution of Eq. 10 is a simple matrix inversion. The application of the inverse Laplace transform on this solution then provides the full time evolution of the density matrix $\tilde{\rho}_{I}(s)\rightarrow\rho_{I}(t)$. For the longitudinal and transverse components, $\rho_{z}(t)=\rho_{\uparrow}(t)-\rho_{\downarrow}(t)$ and $\rho_{\pm}(t)$, this yields

[TABLE]

where $\lambda$ is a real number such that all singularities of the integrand lie to the left of the integration contour. Notice that at this order of approximation the equations for $\rho_{z}$ and $\rho_{\pm}$ decouple. Any cross dependence would require a further expansion of the memory kernel which is beyond our current interest.

The quantities $\rho_{z,\pm}(t=0)$ are the initial states described by the corresponding density matrix. The full time evolution of the reduced density matrix is determined by the complex contour integration along the Bromwich contour shown in Fig. 2. The dynamical behavior is dominated by the location of the integrand’s poles or branch cuts in Laplace space. Such singular structures in Laplace space correspond directly to a collective response of the physical system, and in much of our following discussion we shall read off directly the physical consequences from the position of the poles.

In general we must distinguish between two types of singular structures. The Markov approximation consists of neglecting any $s$ dependence of the $F_{1,2,\pm,z}(s)$ functions and picking up the simple pole at $s=-i[F_{1}(0)+F_{2}(0)]$ in Eq. 14 or $s=-i[F_{\pm}(0)+F_{z}(0)]$ in Eq. 15. If we keep the full $s$ dependence there remains a pole near these values which we shall continue calling the Markov pole. From the residuum theorem this pole leads to a purely exponential decay of the initial state, characteristic for the memoryless Markov behavior. The time scale for this decay is set by the residue of the Markovian pole which corresponds to the spin-lattice relaxation time $T_{1}$ for the longitudinal decay and the spin-spin relaxation or decoherence time $T_{2}$ for the transverse decay.

In addition to the Markov pole, the full $s$ dependence of the denominators can lead to further singularities. At zero temperature, these can take the form of branch cuts but at any finite temperature these cuts split into infinite sequences of poles with a spacing proportional to the temperature. We shall call these poles the non-Markovian poles. Although each pole leads to an exponential decay, the superposition of an infinity of them results in an overall decay that is slower than exponential and this represent the non-Markovian memory effect from the backaction of the bath. This memory effect is eventually suppressed by thermal fluctuations so that the non-Markovian decay is only substantial on times shorter than the memory time $\tau_{T}\sim\hbar/k_{B}T$ set by the temperature. The crossover behavior with respect to the time $\tau_{T}$ can indeed be understood as follows: The separation between the non-Markovian poles is proportional to the temperature $T$. At short times, a large number of these poles contribute to the integral and the time evolution is nonexponential. In particular, in the limit $T\to 0$ the poles merge to a branch cut, causing the typical algebraic decay of a fermionic response. At larger times, at $T>0$, the number of poles contributing with a non-negligible amplitude shrinks until at $t>\tau_{T}$ the decay is essentially described by the exponential arising from the pole closest to the Markov pole. Since the non-Markovian decay is then faster than the Markovian decay the non-Markovian behavior becomes invisible when $t$ passes the time set by the thermal fluctuations.

This description of how the singular structure of the denominator of Eqs. (14) and (15) shapes the time range over which non-Markovian effects are notable is general, and the same qualitative behavior will occur for any bath. To become quantitative we need to focus on a specific type of bath.

III Memory effects in a simple metal

The existence of the rich pole structure discussed above indicates already that the non-Markovian behavior is a universal feature and does not depend on specific system and bath types. To illustrate this generality we show that the effect already appears in the simplest possible fermionic bath, a Fermi gas. Although a simple system the latter fully encodes a many-body response since the spin fluctuations are strongly constrained by the Pauli principle. As such a free Fermi gas is already prototypical to more involved correlated electron systems and has the advantage that all calculations can be performed explicitly. The full system is described by the Hamiltonian of Eq. 5. For simplicity we assume that the system is not magnetized, $\langle S_{x,y,z}\rangle=0$ even if a small magnetic field is applied. This is a reasonable assumption for metals in which $E_{F}$ largely exceeds the Zeeman energy but this also excludes explicitly the Knight shift. As mentioned earlier this makes the discussion of the irreversible behavior of the spin dynamics more transparent, and a later inclusion of the Knight shift is straightforward.

To start, we shall assume that the magnetic field is zero. Due to the full $SU(2)$ symmetry of the combined system, we then have

[TABLE]

which allows us to express all correlators in terms of $F_{z}$ (see Appendix B),

[TABLE]

Here $\langle O\rangle$ denotes the average over the electronic equilibrium state, $\text{Tr}_{\text{el}}[O\rho_{\text{el}}]$, for any operator $O$ and $\{\cdot,\cdot\}$ is the anti-commutator. A detailed derivation of the analytic expression for the spin-spin correlators is presented in Appendix C. Since the electronic Hamiltonian is quadratic, the spin-spin correlators decouple into simple fermionic expectation values and Eq. 17 reduces essentially to integrations over Fermi functions and exponentials. For a conventional metal with an approximately constant density of states about the Fermi energy, the only relevant parameter for the universal low-energy behavior is the electron temperature $T$. This results in two distinct regimes for the spin-spin correlator $F_{z}(s)$. In the region where $\text{Re}(s)<k_{B}T/\hbar$, thermal fluctuations determine the system’s dynamics. On the other hand, for $\text{Re}(s)>k_{B}T/\hbar$, quantum fluctuations are the dominant contributions and memory effects arise from particle-hole fluctuations about the Fermi level. These fluctuations are universal in that they depend only on the existence of the Fermi surface, and we evaluate the correlators of Eq. 17 under the aspect of this universality. However, the general treatment is not limited to this type of environment and as long as a computation of the correlators $F_{j}(s)$ is possible the method can be applied to any other bath as well.

The correlators $F_{j}(s)$ are proportional to the interaction $A^{2}$ due to the Born approximation of the memory kernel $\tilde{\Sigma}_{I}(s)$. In the following we use a small expansion parameter $\alpha$ set by the interaction strength $A$ and the density of states at the Fermi level $\nu_{0}$

[TABLE]

Generally, $\nu_{0}\sim E_{F}^{-1}$ with $E_{F}$ being the Fermi energy, such that $\alpha\sim\left(A/E_{F}\right)^{2}$. We assume that $\alpha$ is small; in fact, we already assumed the interaction strength $A$ between system and environment is small, such that the Born approximation for the equation of motion is justified. Compared to the Kondo problem we are thus in a regime, where the Born approximation does not break down and the effective coupling does not diverge. The spin-spin correlation functions for a bath modelled as a Fermi gas up to $O(\alpha)$ (see Appendix C for details) can be written as

[TABLE]

with

[TABLE]

where $k_{B}$ is the Boltzmann constant, $\psi(z)=~{}[\frac{d}{dz}\Gamma(z)]/\Gamma(z)$ is the digamma function with $\Gamma(z)$ being Euler’s Gamma function, and $\xi\sim E_{F}$ is a high energy cutoff on the order of the Fermi energy (or the bandwidth). The details on the derivation can be found in Appendix C. Although $F(s)$ and $G(s)$ are almost identical, we shall need both of them to unambiguously express the correlators at nonzero magnetic field (see below).

Equations (20) and (21) are universal in that they capture exactly the low-energy fluctuations of the electronic excitations. All high-energy fluctuations depending on the nonuniversal details of the band structure are absorbed in the cutoff $\xi$ within the logarithms. The latter contribute only weakly to the dynamics of the spin $I$ at very short times $t<\hbar/\xi$. This allows us to focus on the universal behavior at longer times in the following, fully encoded in $F(s)$ and $G(s)$.

The existence of quantum critical correlations in $F_{z}(s)$ is best visible in the limit $T\to 0$. Then the digamma function $\psi$ has the asymptotes of a logarithm and the $T$ dependence in the logarithm cancels, which leads to

[TABLE]

This result also conveniently bypasses a direct $T=0$ calculation, which would be much more involved. A logarithmic behavior of spin-spin correlators as in Eq. 22 signals also the existence of an underlying Kondo effect. However, as mentioned earlier, by maintaining a small $\alpha$ and preparing a well-defined state $\rho_{I}(t=0)$, the onset of Kondo physics even at very low temperatures is either absent or would occur at rather long times. Therefore, even in then regime where the correlation functions reflect the effective $T\to 0$ limit through their logarithmic form, the temperatures are still high enough to neglect any Kondo physics.

The coherence disappears with increasing temperature, and indeed at large $k_{B}T\gg\text{Re}(s)$ the digamma function $\psi$ tends to a constant that can be absorbed in the cutoff $\xi\to\xi^{\prime}$, and we find that

[TABLE]

The $s$ independent term corresponds to the Markov approximation and will thus cause the standard exponential decay. The second term provides the further renormalization from the electron fluctuations. In this high-temperature limit, this correction is proportional to $s$ and thus causes only a reduction of the Markovian decay amplitude. But for $\text{Re}(s)>k_{B}T$, which means for times $t<\hbar/k_{B}T$, the effects caused by the quantum correlations cause a significant deviation from the Markovian exponential behavior.

Before evaluating this dynamics through Eqs. 14 and 15, we should recall that this $F_{z}(s)$ incorporates the dynamics at zero field. An investigation for $b_{z}^{I}\neq 0$ requires us to reintroduce field-dependent phase factors in the correlation functions, leading to shifts of the arguments of $F(s)$ and $G(s)$. From Eqs. 83, 84, 85, 86 and 87, we have

[TABLE]

with $F_{z}(s)$ as given in Eq. 23. The full time evolution of the reduced density matrix is obtained through Eq. 14 and Eq. 15. In the latter equations, the contour integration is evaluated through the zeros of the denominators $s+iF_{1}(s)+iF_{2}(s)$ and $s+iF_{\pm}(s)+iF_{z}(s)$, respectively. Since $F_{z,1,2,\pm}$ are proportional to $\alpha$, there are two types of zeros. First, there is a zero at $s\sim\alpha$, leading to what we have called the Markov pole in the previous section. Second, there is a sequence of zeros near the poles of the digamma function $\psi$. Indeed, in the vicinity of a pole, $\psi$ can raise to $\psi\sim 1/\alpha$ such that $\alpha\psi\sim 1$, which in turn can compensate for the remaining terms in the denominator. Since the digamma function $\psi(z)$ has poles at $z=0,-1,-2,\dots$, this leads to a dynamics governed by the quantum fluctuations of the Fermi gas expressed by the $\psi$ term, and above we called these poles the non-Markovian poles. The resulting behavior will be calculated explicitly in the next subsections.

For $\rho_{z}$, it should also be noted that, in Eq. 14, there is a further pole at $s=0$. This corresponds to the equilibrium value $\langle I_{z}\rangle=I\rho_{z}(t\to\infty)$ and we will include it in the discussion of the Markov behavior.

III.1 Markovian decay

The Markov approximation neglects the $s$ dependence of the spin-spin correlation functions Eqs. 24, 25, 26, 27 and 28 by assuming that the system does not have a memory time, i.e., the equation of motion Eq. 7 is time local. The corresponding Markov pole is on the order of $s\sim\alpha$ and has a negative real part. This leads to an exponential decay of the density matrix components governed by the decay times $T_{1,2}=-1/\mathrm{Re}(s)$, with the relaxation time $T_{1}$ and the dephasing time $T_{2}$ describing the evolution of $\rho_{z}$ and $\rho_{\pm}$ respectively.

We have thus to solve for $s+iF_{1}(s)+iF_{2}(s)=0$ and $s+iF_{+}(s)+iF_{-}(s)=0$ for general $b_{z}^{I}$ under the condition $s\sim\alpha$. As a general strategy, we shall profit from the smallness of $\alpha$ to expand the position $s$ of the poles and their residues to order $\alpha$. However, to capture the full time dynamics nonperturbatively, we must not expand the time-dependent exponentials $e^{st}$. The details of this calculation for the Markov part are given in Appendix D. From the real part of the corresponding Markov poles, we obtain the decay times

[TABLE]

In the zero field limit $b_{z}^{I}=0$, these expressions become

[TABLE]

which recovers the Korringa relation $TT_{1}=\kappa$ of Fermi liquids Korringa (1950); Moriya (1963); Abragam (1961); Slichter (1990), with $\kappa=\hbar/4\alpha\pi k_{B}$ being the Korringa constant. Expressed in $\kappa$, Eq. 29 coincides with the magnetic field dependent decay times found in the literature Pobell (2007); Slichter (1990); Abragam (1961). Furthermore, the equality $T_{1}=T_{2}$ is the consequence of the $SU(2)$ symmetry of the system. This is in contrast to systems with a broken $SU(2)$ symmetry, in which $T_{2}=2T_{1}$ can be achieved Golovach et al. (2004); Borhani et al. (2006); Slichter (1990); Abragam (1961). The latter result can also obtained in the present case by introducing symmetry breaking through the magnetic field. Indeed, in the limit $b_{z}^{I}\gg k_{B}T$, Eqs. 29 and 30 become $T_{1}=\hbar/2\alpha\pi b_{z}^{I}$ and $T_{2}=\hbar/\alpha\pi b_{z}^{I}$, and we recover $T_{2}=2T_{1}$. In Fig. 3, we show the universal behavior of $TT_{1,2}/\kappa$ as a function of $b_{z}^{I}/2k_{B}T$. It is notable that $T_{2}$ fulfils for any nonzero field $T_{1}<T_{2}<2T_{1}$ and that the upper limit $T_{2}=2T_{1}$ is reached only rather slowly for large $b_{z}^{I}$ values, which is best seen through the ratio $T_{2}/T_{1}$ shown in the inset of Fig. 3.

From Appendix D, we obtain then the Markovian part of the time evolution

[TABLE]

where $\rho_{z}^{\text{eq}}$ is the paramagnetic equilibrium magnetization corresponding to complete thermalization with the electronic bath,

[TABLE]

and where the transverse component contains a spin precession described by the frequency

[TABLE]

In the latter expressions, $\xi^{\prime}$ is an inessential renormalization of $\xi$ obtained by absorbing a constant $\psi(1)/2$ in the cutoff. Through a similar further shift of the cutoff to $\xi^{\prime\prime}$, we can rewrite the zero field expressions as

[TABLE]

with $\xi^{\prime\prime}=\xi\exp(-\psi(1))$.

Notice that in contrast to the standard Markovian decay we have kept an $O(\alpha)$ correction in the amplitudes of Eqs. 32, 33 and 36. This results from keeping the $s$ dependence in $F_{z,1,2,\pm}(s)$ instead of setting in the latter functions $s=0$ and solving, for instance, $s+iF_{1}(0)+iF_{2}(0)=~{}0$. As a consequence, the weight of the amplitudes is reduced from 1 to $[1-O(\alpha)]$. This effect alone already indicates the presence of the further non-Markovian decay terms and that the latter have an amplitude of order $\alpha$.

The Markovian decay according to Eqs. 32 and 33 is shown in Fig. 4 for initial conditions $\rho_{z}(t=\hbar\xi^{-1})=1$ and $\rho_{\pm}(t=\hbar\xi^{-1})=1$, respectively. In the zero-field limit, the evolution of $\rho_{z}(t)$ and $\rho_{\pm}(t)$ are identical (orange curve). For $b_{z}^{I}>0$, both $T_{1}$ and $T_{2}$ times decrease. In the figure, we show the decay for a rather large field $b_{z}^{I}=2\pi k_{B}T$ in which $\rho_{z}^{\text{eq}}=-\tanh(1)\approx-0.76$ is significant (blue curve), and the effect of $\exp(i\omega_{\pm}t)$ in $\rho_{\pm}^{\text{M}}$ becomes visible since $\omega_{\pm}\sim 1/T_{2}$ (purple dashed curves). Such oscillations renormalize the normal precession from the spin in the magnetic field that we have removed by going to the rotating frame. As such, their role is similar to the Knight shift, which we have neglected in the present treatment but which also adds to the eigenfrequency. But in contrast to the Knight shift, the $\omega_{\pm}$ have a strong nonlinear magnetic field dependence. They are only proportional to $b_{z}^{I}$ at low fields but for $b_{z}^{I}>2\pi k_{B}T$ are strongly sub-linear and eventually, at very strong fields, change sign.

III.2 Non-Markovian contributions

Apart from the isolated Markov pole, we find an infinite sequence of poles close to the singularities of the spin-spin correlators Eqs. 24, 25, 26, 27 and 28. While for the Markov poles $s\sim\alpha$, the $s$ values of these new poles are to leading order independent of $\alpha$. Expressions like $s+iF_{1}(s)+iF_{2}(s)$ can then only be zero if a divergence of order $1/\alpha$ compensates their small $\alpha$ amplitude. This means that the $s$ have to lie near the singularities of $F$ and $G$. For $F(s)$ and $G(s)$ as given by Eqs. 20 and 21, the divergences are governed by the singularities of the digamma functions $\psi\left(1+s/2\pi k_{B}T\right)$ and $\psi\left(1+(s\pm ib_{z}^{I})/2\pi k_{B}T\right)$. The digamma function $\psi(1+z)$ has simple poles at $z=-1,-2,\dots$ and in their vicinity we find the necessary $\psi(1+z)\sim 1/\alpha$ behavior. This allows us to systematically expand the denominators of Eqs. 14 and 15 and pick up the residues to order $\alpha$. The detailed calculation is done in Appendix E. Remarkably, we can do the summation over the residues exactly and the resulting non-Markovian contributions to $\rho_{z}$ and $\rho_{\pm}$ are given by

[TABLE]

with $\tau_{T}=\hbar/2\pi k_{B}T$ being the thermal time and $h(t,b_{z}^{I})$ being the function

[TABLE]

Notice that since $\ln(1-e^{-t/\tau_{T}})\sim e^{-t/\tau_{T}}$ for $t>\tau_{T}$ the characteristic decay time for these non-Markovian parts is $\tau_{T}$ and not $T_{1,2}$ which would be by $1/\alpha$ longer. This is not surprising as temperature fluctuations erase any memory effect. However, it is notable is that we have obtained an explicit prescription of the non-Markovian quantum behavior within the memory time. At low magnetic fields, the scale-free nature of the Fermi gas becomes manifest in that temperature is the only scaling parameter for the dynamic behavior. Larger magnetic fields introduce a further scale to the system, but we see that their effect is limited, entering in a nontrivial way only through $h(t,b_{z}^{I})$. The latter causes either only a small renormalization of the overall decay for $b_{z}^{I}<4\alpha\pi k_{B}T$, or the addition of further, rather slow oscillations for $b_{z}^{I}>4\alpha\pi k_{B}T$. These oscillations are similar to the Knight-shift-type effect of $\omega_{\pm}$ for the Markov contribution but they appear here for both $\rho_{\pm}^{\text{nM}}$ and $\rho_{z}^{\text{nM}}$ (see Fig. 5), and with different onset fields and frequencies due to the dependence through $h(t,b_{z}^{I}/2)$ and $h(t,b_{z}^{I})$, respectively. At zero field, the evolution of $\rho_{z}^{\text{nM}}$ and $\rho_{\pm}^{\text{nM}}$ coincide (Fig. 5, orange curve). It should finally be noted that Eq. 39 has been derived under the assumption that the maximally considered $b_{z}^{I}$ does not excessively exceed $2\pi k_{B}T$. Nevertheless, we expect that (39) remains a good approximation even in the limit $b_{z}^{I}\gg 2\pi k_{B}T$ because corrections to the derivation done in Appendix E can only weakly renormalise the position of the non-Markovian poles.

III.3 Total decay and limits

The full time evolution of the reduced density matrix is the sum of the Markovian and the non-Markovian contributions, given for $\rho_{z}$ by Eqs. 32 and 37 and for $\rho_{\pm}$ by Eqs. 33 and 38. These results cover the universal dynamics for all times $t>\hbar/\xi$, with a domination of the non-Markovian decay for $t<\tau_{T}=\hbar/2\pi k_{B}T$ and a crossover to the standard Markovian exponential decay at $t>\tau_{T}$. For times $t<\hbar/\xi$, the evolution is non-universal but from an expansion of Eq. 7 around $t=0$ for small times we see that the onset of the decay is quadratic in time, $\rho\approx 1-i\Sigma(0)t^{2}/2$, and hence causes only a very small lowering of the amplitudes before the onset of the universal behavior. The latter starts logarithmically $\sim\ln(\xi t/\hbar)$ for both $\rho_{z}$ and $\rho_{\pm}$. This is indeed the signature of a Fermi edge singularity many-body reaction Mahan (1967a, b); Anderson (1967); Nozières and De Dominicis (1969) that is triggered by the local spin-spin interaction, and this behavior would turn into Kondo correlations if the interaction could cause an arbitrary number of spin flips Anderson and Yuval (1969); Anderson et al. (1970).

In Fig. 6, we show the full decay of the different components in the cross-over region from non-Markovian to Markovian behavior around $t\sim\tau_{T}$. Figure 6(a) shows the comparison of $\rho_{z}=\rho_{\pm}$ at $b_{z}^{I}=0$ with a pure Markovian exponential decay $\rho(t)=e^{-t/T_{1,2}}$. Since the Markovian decay times $T_{1,2}$ are on the order of $\tau_{T}/\alpha$ (far outside the plotted range) any focus on times around $T_{1,2}$ misses the non-Markovian behavior. However, since the non-Markovian amplitude is on the order of $\alpha$, the amplitude of the Markovian decay is reduced from $1$ to $[1-O(\alpha)]$. The fast non-Markovian decay thus results in a fast initial slip Geigenmüller et al. (1983); Haake and Reibold (1985); Suárez et al. (1992); Yu et al. (2000); Hell et al. (2016), a systematic offset of the amplitude of $O(\alpha)$ at all times $t>\tau_{T}$, which we show at zero field in Fig. 6(a) and at nonzero field in Fig. 6(b). The offset is detectable by extrapolating the Markovian behavior back to time $t=0$ [light dashed line in panel (a)]. In panel (c) we use a logarithmic time scale to enhance the visibility of the non-Markovian features at a field $b_{z}\sim k_{B}T$ (solid lines), and at a field $b_{z}\gg k_{B}T$ (dashed lines) at which the oscillations in $h(t,b_{z}^{I})$ become visible.

IV Manipulating the environment:

Cooling protocol

The notable feature of the non-Markovian decay is that it is fast and initially always temperature independent. It may thus be possible to put this property to use in an application that depends on the decay of the spin system but is in practice limited by the diverging $T_{1,2}$ times at very low temperatures. Such a situation indeed arises in cooling techniques based on a cold spin system. Adiabatic demagnetization is a standard technique that allows reaching very low temperatures of a spin system such that the latter can be used as a refrigerant for cooling another system Pobell (2007). The technique relies on the fact that the entropy of the spin system depends on the ratio of magnetic field and spin temperature, $B/T_{I}$, such that the adiabatic reduction of an initially strong magnetic field leads to lower spin temperatures $T_{I}$. But the efficiency as a refrigerant then depends further on the thermal equilibration between the spin system and the rest of the system, including the lattice as well as the electrons. This process relies on the relaxation of the spin system, which is governed by the $T_{1}$ time. The Korringa relation $T_{1}\propto 1/T$ provides the crucial bottleneck in that the increasingly long relaxation times at low temperatures make the cooling ineffective against heat leaks. For instance, nuclear spins in rhodium can be cooled to temperatures below 100 pK but the electron and phonon temperatures are limited to about 0.1 mK, likely due to this effect Knuuttila et al. (2001); Hakonen et al. (1993). In other bulk metals such as platinum, electron temperatures of $1.5$ $\mu$K can be reached Wendler et al. (1997), but also here a notable discrepancy with the nuclear spin temperature of 0.3 $\mu$K persist. More challenging is the cooling of semiconductors or nanostructures because of their much lower thermal conductivity, and the lowest reached electron temperatures are here in the mK range Clark et al. (2010); Scheller et al. (2014); Palma et al. (2017a, b); Zyuzin et al. (2014).

In the following, we thus examine the possibility of using the temperature-independent non-Markovian decay to transport heat out of the electron bath to speed up the slow thermal relaxation process described by the Markovian decay. Quantum thermodynamic cooling protocols explicitly based on nonequilibrium physics have been suggested before (see, e.g., Refs. Härtle et al. (2018); Chen et al. (2010); del Campo et al. (2014); Funo et al. (2019)). We complement these proposals by using our explicit knowledge of the non-Markovian decay to design a protocol with a strong focus on speeding up the cooling process. We should also note that such a cooling protocol goes beyond a thermodynamic cycle because intensive thermodynamic quantities such as the spin temperature cannot be defined during the quantum coherent evolution. Since the non-Markovian decay has an amplitude of order $\alpha$, this approach needs repeated reinitialization of the spin system. This could be done, for instance, through femtosecond resolved optical pumping of electron spins Kikkawa and Awschalom (1998) or nuclear spins Salis et al. (2001a, b), by optical pumping of hole spins Gerardot et al. (2008), or by partial measurements Blok et al. (2014).

In Fig. 7, this idea is sketched. We show the heat $\Delta Q$ transported out of the electronic system as a function of time. We will see in Eq. 42 below that $\Delta Q$ is proportional to the decay of $\rho_{z}$. Therefore, $\Delta Q$ has a fast initial non-Markovian decay followed by a slow relaxation set by the relaxation time $T_{1}$ (dark purple curve). Through a repeated resetting of the spin polarization at the end of the non-Markovian regime, as shown in the inset, we can repeat the fast initial drop of $\Delta Q$ and transport heat out of the system more quickly (bright orange curve).

The curve shows, of course, an idealized situation and the cooling effect must be complemented by the possible repeat time, the reheating by the pumping, and the influence of external heat leaks. Our goal in the following analysis is to demonstrate under which conditions such a cooling protocol can become effective, and to show that this can be achievable with state-of-the-art techniques.

Instead of a single impurity spin $I$, we require now a macroscopic ensemble of such spins. But if the direct interaction between these spins is weak, we can treat each spin individually and use the results found in the previous sections. Without spin-spin interaction, the impurity spin Hamiltonian $H_{I}$ becomes

[TABLE]

which is the same as the Zeeman term in Eq. 1 but now with the index $i=1,\dots,N_{I}$ labeling the impurity spins $\mathbf{I}_{i}$. Under the assumption of independent spins and the assumption that we can neglect the coupling between the impurity spins mediated through the electron system, we can write the reduced density matrix for the spin system as the product $\rho_{I}=\rho_{I,1}\otimes\dots\otimes\rho_{I,N_{I}}$. Assuming identical initial conditions, the time evolution of $\rho_{I}$ consists of $N_{I}$ copies of the same evolution of a single spin.

The energy current $J_{I}$ of the spin system can then be defined by

[TABLE]

such that $J_{I}>0$ corresponds to an energy flow from the environment into the spin system. The heat $\Delta Q(t)$ transferred up to a time $t$ is given by

[TABLE]

where $\rho_{z}(t)$ is the decay of a single representative spin and the factor $N_{I}$ takes the ensemble into account. To focus on the non-Markovian decay in $\Delta Q$, we consider time scales $t\ll\tau_{T}$ at which, from Eq. 37, the dynamics of $\rho_{z}(t)$ becomes temperature independent and leads to

[TABLE]

with $Q_{z}=N_{I}\mu_{n}g_{I}B_{z}I=N_{I}b_{z}^{I}$ and $\rho_{0}=[\rho_{z}(0)-\rho_{z}^{\text{eq}}]$.

The spins absorb energy if $\Delta Q<0$ and therefore the spins should initially minimize their Zeeman energy such that $B_{z}[\rho_{z}(0)-\rho_{z}^{\text{eq}}]<0$. For $B_{z}>0$, this means we should choose $\rho_{z}(0)\approx-1$, corresponding indeed to the ground state of ultracold spins, instead of the highly excited $\rho_{z}(0)=+1$ chosen for illustration in the previous sections. We then obtain $\rho_{0}=-[1-\tanh(b_{z}^{I}/2\pi k_{B}T)]=-|\rho_{0}|$. For an optimal cooling, the amplitude $B_{z}|\rho_{0}|$ should be as large as possible, but large $B_{z}$ lead to a small $|\rho_{0}|$. Since $|\rho_{0}|$ depends on the ratio $b_{z}^{I}/2\pi k_{B}T$, the optimum is obtained by tuning $b_{z}^{I}$ to the order of $2\pi k_{B}T$ and consequently $B_{z}|\rho_{0}|\sim T$.

We now want to consider a cooling cycle where we utilize the temperature-independent non-Markovian decay to transfer heat from the electronic bath into the impurity spin system. The relevant time scale for the duration of one cycle is $\Delta t+\tau_{r}$, where $\Delta t\ll T_{1}$ is the time interval of the spin decay in the fast non-Markovian regime and $\tau_{r}$ is the time needed to reinitialize the impurity spin. In the inset of Fig. 7, a sketch of the cooling cycle and the corresponding state of the impurity spin on the Bloch sphere is shown. In the short time $\Delta t$, the spin decays by a small amplitude proportional to $\alpha$, according to Eq. 37. After the initial slip, the impurity spin needs to be projected back onto its initial state, during which the system is necessarily reheated by an amount $Q_{r}$.

The total heat $Q_{\text{el}}$ transferred out of the electron system over $N_{P}$ cooling cycles can then be written as

[TABLE]

where $-\Delta Q(\Delta t)$ is the non-Markovian cooling per cycle, $Q_{r}$ is the heat deposited during each reinitialization process, and $tJ^{Q}_{\text{ext}}$ takes into account external heat leaks generating a continuous inflowing heat current $J^{Q}_{\text{ext}}$. The total time $t$ of the process is $N_{P}(\Delta t+\tau_{r})$.

Cooling is possible if $Q_{\text{el}}<0$. By dividing Eq. 44 by $N_{P}|4\alpha Q_{z}\rho_{0}|$ and absorbing $Q_{r}/|4\alpha Q_{z}\rho_{0}|$ in the logarithm, we can write this condition as

[TABLE]

where

[TABLE]

is a dimensionless quantity measuring the efficiency of the cooling protocol,

[TABLE]

measures the influence of the external heat leaks, and

[TABLE]

sets the characteristic time for the protocol. The time $\tau_{0}$ should fulfill $\tau_{0}<\tau_{T}$ as only then the non-Markovian decay is effective. This provides a condition on $Q_{r}$,

[TABLE]

Since $Q_{z}\rho_{0}$ should be chosen proportional to $T$, this condition puts a not surprising constraint on the lowest reachable temperatures set by the heating by reinitialization. Since the dependence of the protocol on the reinitialization time $\tau_{r}$ is less important, an optimization through the interaction time $\tau_{r}$ should be possible. Notice also that a tuneability exists through the electron density as $\xi\sim E_{F}$.

The influence of the external heat leaks is assessed through Eq. 45. Figure 8 shows $q_{\text{cp}}$ as a function of $\Delta t/\tau_{0}$ in comparison with a choice of $q_{\text{ext}}$. The shaded area marks the parameter region where Eq. 45 is fulfilled and cooling the system with the pumping protocol is possible. There is naturally a maximum heat leak $q_{\text{ext}}=q_{\text{ext}}^{\text{max}}$ beyond which cooling is no longer possible, indicated by the maximum of the curve. Maximizing the amount of heat that can be carried out of the electronic system throughout the whole cooling process leads to an optimal time $\Delta t$ between pulses

[TABLE]

where $e$ is Euler’s number and $W(z)$ is the Lambert $W$ function, defined as the inverse of the function $z(W)=W\exp(W)$. If $z=\tau_{r}/e\tau_{0}\ll 1$, large heating by the pumping is implied and $W(z)\approx z/e$, i.e., $\Delta t\sim\tau_{0}$. In the opposite limit $z=\tau_{r}/e\tau_{0}\gg 1$, we have $W(z)\approx\ln(z)$, which leads to $\Delta t\sim\tau_{r}$. Upon the constraint of Eq. 49, we see that $\Delta t_{\text{opt}}$ is thus within the range set by $\tau_{r}$ and $\tau_{T}$.

Taking Eq. 44, we can estimate how many pump cycles are possible by looking at the temperature of the electronic system $T_{\text{el}}(t)$, assuming an equilibrium state of the bath itself,

[TABLE]

where $C_{\text{el}}$ is the specific heat for the electronic system. Before we start pumping the system, we can assign an initial temperature $T_{I}(0)$ to the impurity spins. Although through the pumping and reinitializations the concept of a spin temperature $T_{I}$ becomes no longer meaningful, we shall use $T_{I}(0)$ as a lower bound, but due to the temperature dependence of $q_{\text{cp}}$ and $\tau_{0}$ it may be an optimistic bound for the final electronic temperature. Setting thus the minimum temperature the electrons can reach to $T_{\text{el}}(t)=T_{I}(0)$ and assuming that the temperature dependence of the parameters is otherwise weak, we are led to

[TABLE]

During each cycle, a small amount of heat proportional to $\alpha$ is transferred from the electrons into the impurity spin system via spin flips. Thus, the efficiency of the cooling process relies on a fast repetition such that $N_{P}\sim 1/\alpha$ repetitions can be made faster than the $T_{1}$ time.

The small parameter $\alpha=(\nu_{0}A)^{2}$ is set by the ratio of the coupling constant $A$ and the Fermi energy, but $\nu_{0}$ contains also the ratio $n_{e}/n_{I}$ of electron density $n_{e}$ over impurity spin density $n_{I}$ Simon and Loss (2007); Simon et al. (2008); Braunecker et al. (2009a, b). The parameter is thus highly dependent on the considered system and can vary from very small $\alpha\sim 10^{-10}$, as is characteristic for nuclear spins in bulk metals Pobell (2007), to $\alpha\sim 10^{-4}$ for paramagnetic spins in correlated metals Bonville et al. (1991), and to $\alpha>1$ in magnetic semiconductors Okabayashi et al. (1998); Jungwirth et al. (2006), which is, however, beyond the validity of our approach.

As an example, let us consider a semiconductor with nuclear spins, e.g., GaAs, with $A=90$ $\mu$eV and $E_{F}$ on the order of meV, such that $\alpha\sim 10^{-8}$. For such small $\alpha$, a direct observation of the non-Markovian behavior is elusive but it provides a good illustration that the cooling protocol could even then become useful. The corresponding $T_{1}=\hbar/4\alpha\pi k_{B}T$ is at $T\sim 0.1$ mK on the order of seconds. The suggested heat transfer protocol becomes effective if we can repeat it $N_{P}\sim 1/\alpha$ times within the $T_{1}$ time. Thus, the maximum time $\Delta t+\tau_{r}$ between two pulses should be limited by $\Delta t+\tau_{r}<\alpha T_{1}\sim 10$ ns, which means that both $\Delta t$ and $\tau_{r}$ can be on the order of 10 ns.

Recall that $\tau_{0}$ in Eq. 48 is given by $\hbar/\xi\sim\hbar/E_{F}$, which is in the picosecond range, times an exponential. The latter is then allowed to grow to $10^{3}$-$10^{4}$ such that the requirement on the exponent $Q_{r}/|4\alpha Q_{z}\rho_{0}|$ is not too stringent. Furthermore, a reinitialization time $\tau_{r}\sim$ 10 ns is short but not beyond the reach of modern experimental techniques. We must also recall that it corresponds to a reinitialization of a tiny decay of amplitude $\alpha$ such that the limiting factor will be given by the electronics for the repetition rate rather than the physical manipulation of the spin.

Moreover, tuning $\alpha$ to larger values works in favour of this cooling protocol. Instead of influencing $\alpha$ through $E_{F}$, we can also substantially enhance it by considering different materials as outlined above. Values of $\alpha\sim 1$ can become rather common, but we must emphasize that our theory depends on the smallness of $\alpha$ and that at the resulting small $T_{1}$ times the standard cooling techniques may remain more effective.

V Conclusions

In this work, we discussed a unifying analytical approach to nonequilibrium quantum thermodynamics, coherent non-Markovian evolution, and backaction from the system on the bath for quantum spins embedded in an electron conductor. We addressed three goals in particular: (i) the presentation of a general framework allowing us to obtain the full dynamics of a spin system in a bath with non-Markovian memory effects; (ii) the application when the bath is an itinerant fermionic system; and (iii) the proposition of a cooling protocol based on the temperature insensitive non-Markovian quantum fluctuations. The development of the framework was based on a careful expansion of the generalized master equation for the spin’s density matrix and on an extraction of the physical consequences from the pole structure in the complex Laplace space. The memory effect arising from the backaction of the fermionic system on the spin dynamics is indeed governed by the coherent excitations of the fermionic bath, which each are represented by one of the poles. The collective of the poles then defines the non-Markovian dynamics, and we have presented a systematic method for accessing the dynamics. This method is limited to a small coupling strength between the spin and the bath degrees of freedom, encoded in the parameter $\alpha$ defined in Eq. 18. But the method is not limited to the considered fermionic system, and an extension of this work to strongly correlated systems will be of special interest. In the latter systems, the internal interactions lead to strong time correlations and therefore to an enhanced non-Markovian contribution to the spin dynamics. An analysis of the memory effects will thus provide a time-resolved access to the strong correlation physics. Furthermore, the method is not restricted to fermionic environments. As long as the bath correlators can be computed, a similar analysis of the pole structure can follow for a bosonic environment as well.

As an archetype for a correlated response, we have considered the example of a free itinerant electron gas. Such an electron gas resembles a critical correlated system in the sense that its quantum and thermal fluctuations are governed by a small number of parameters, here the electron temperature $T$ and the Fermi energy $E_{F}$. For such an electronic bath, we have provided an explicit solution to the full free induction decay of the spin system, crossing over from the initial fast non-Markovian decay to the conventional Markovian exponential decay. The latter is characterized by the times $T_{1}$ and $T_{2}$, whereas the non-Markovian decay occurs on the much shorter electronic thermal time and remains detectable at longer times through the initial slip, the systematic offset of the Markovian decay from the expected initial value at time zero.

At short times, the non-Markovian decay is temperature independent and determined entirely by quantum fluctuations. Since this decay can be tuned to transport heat from the electron to the spin system, we have finally proposed a cooling protocol based on a repeated triggering of the non-Markovian decay by reinitializations. Such a method could overcome the diverging time scale $T_{1}\propto 1/\alpha T$ that limits the efficiency of adiabatic demagnetization cooling. Although the method is probably not suitable to concurrence with the standard cooling methods for bulk metals, we have provided an estimate that it could become effective in cooling semiconductor structures. For semiconductors, the $T_{1}$ times are much longer and are responsible for the bottleneck that has prevented experiments from reaching electron temperatures below a millikelvin so far.

Acknowledgements.

We thank S. Parameswaran, R. Whitney, F. Pollock and T. Stace for useful discussions. S.M. acknowledges the support by the EPSRC under Grant No. EP/L015110/1. D.L. and J.K. acknowledge the support by the Swiss National Science Foundation and NCCR QSIT. The work presented in this paper is theoretical. No data has been produced and supporting research data is not required.

Appendix A Master equation for reduced density matrix

To gain access to the temporal correlations or memory effects, we use the projection operator method Breuer and Petruccione (2007); Fick and Sauermann (1990). In the following, we provide a short derivation of Eq. 7 in the main text. The starting point is the Nakajima-Zwanzig equation, a generalized exact master equation for the density matrix $\rho(t)$

[TABLE]

with the memory kernel $\Sigma(t-t^{\prime})$ . The superoperator $L$ is the Liouvillian defined by the Hamiltonian $H$, i.e., $LO=~{}[H,O]$ for any operator $O$. The memory kernel,

[TABLE]

captures the whole history of the system up to time $t$. The projection operators $P$ and $Q$ obey $P+Q=~{}\mathbb{1}$, $P^{2}=P$, $Q^{2}=Q$, and $PQ=0$. We choose the projection operator $P$ as $PO=~{}\rho_{\text{el}}\otimes~{}\text{Tr}_{\text{el}}\left[O\right]$, with $\rho_{\text{el}}$ being the equilibrium density matrix of the bath. The trace over the bath degrees of freedom $\text{Tr}_{\text{el}}[\cdot]$ defines the reduced density matrix describing the spin system $\rho_{I}=\text{Tr}_{\text{el}}\left[\rho\right]$. Splitting the Hamiltonian into $H=H_{0}+H_{\text{int}}$, we can write the Liouvillian $L=L_{0}+L_{\text{int}}$ accordingly. Then, the first term in Eq. 53 including $L_{0}$, which describes the eigendynamics of the spin system, drops out due to the choice of the rotating frame of reference. The part containing $PL_{\text{int}}P$ vanishes since the system is spin conserving and does not acquire a net magnetic moment. Using the definition of the projector $P$ leads to the master equation, Eq. 7, solved in the main text for the reduced density matrix $\rho_{I}$

[TABLE]

with the reduced memory kernel

[TABLE]

To solve the master equation in Eq. 55 for the full time evolution we analyze it in Laplace space. The Laplace transform of the generalized master equation, $f(t)\to~{}\tilde{f}(s)=~{}\int_{0}^{\infty}dt\,e^{-ts}f(t)$, with $\text{Re}(s)>0$, is given by

[TABLE]

with the identity operator $\mathbb{1}$, and the Laplace transform of the reduced memory kernel $\tilde{\Sigma}_{I}(s)$ is

[TABLE]

The operator $\left(s\mathbb{1}+iLQ\right)^{-1}$ obeys a Schwinger-Dyson equation and iteration of this identity provides an expansion in powers of the interacting part $L_{\text{int}}$ of the full Liouvillian Fick and Sauermann (1990); Coish and Loss (2004)

[TABLE]

For odd $n$, the corresponding term in the memory kernel $\tilde{\Sigma}_{I}(s)$ is zero since the Hamiltonian $H_{0}$ is spin conserving. In the main text, we used the Born approximation up to second order in the interaction. This corresponds to the first term in the series expansion in Eq. 59.

Appendix B Expression for the memory kernel $\tilde{\Sigma}(s)$

Assuming spin-1/2 for the impurity spin, every operator $O$ can be decomposed in the basis of the spin system spanned by $\{I_{\uparrow},I_{\downarrow},I_{-},I_{+}\}$. The operator $O$ can then be written as

[TABLE]

Using the impurity spin basis in the rotating frame of reference, the superoperators, $[L_{0}],\ [L_{\text{int}}],\ [\tilde{\Sigma}(s)]$, can be expressed as a $4\times 4$ matrix. Here, we denote the matrix representation of the superoperators with $[\cdot]$. To find expressions for the superoperators, we use lower case operators which act only on the bath. The Hamiltonian for the bath can then be written as $H_{0}=h_{0}\otimes\mathbb{1}_{I}$. For the interacting part, we decompose the spin-spin interaction

[TABLE]

In the case of the spin-spin interaction discussed here, $H_{\text{int}}=~{}A\mathbf{S}\cdot\mathbf{I}$ the bath operators are $h_{\uparrow}=-h_{\downarrow}=\pm AS_{z}/2=h_{z}$ and $h_{\pm}=2AS_{\pm}$. Within the basis of the spin-1/2, $\{I_{\uparrow},I_{\downarrow},I_{-},I_{+}\}$, the interaction part of the Liouvillian $[L_{\text{int}}]$ is given by

[TABLE]

Here, superscript $-$ refers to the commutator, $L_{h_{z}}^{-}o=~{}[h_{z},o]$, and the superscript $+$ refers to the anticommutator, $L_{h_{z}}^{+}o=~{}\{h_{z},o\}$. $L,\ R$ denote if the operator acts from the left- or right-hand side. Finally, the memory kernel within the spin system’s basis can be expressed as

[TABLE]

The entries of the memory kernel $F_{1,2,z,\pm}(s)$ with $\text{Re}(s)>0$ are the Laplace transforms of spin-spin correlation functions

[TABLE]

for $t\geq 0$. For a system with $SU(2)$, symmetry all these expressions are equal, and an example for such a situation is the simple Fermi gas considered in the main text. But it should be emphasized that the $SU(2)$ symmetry of the Fermi gas holds only in the laboratory frame since the transformation to the rotating frame Eq. 2 explicitly breaks it. For the evaluation of the fermionic correlators, it is therefore important to go back to the laboratory frame, in which the latter equations become

[TABLE]

Here, $b_{z}^{I}$ is the prefactor for the Zeeman term of the impurity spin in Eq. 1. Even with the underlying $SU(2)$ symmetry, these functions are generally distinct, and only at $b_{z}^{I}=0$ will the $SU(2)$ symmetry of the Fermi gas cause $F_{1}(s)=~{}F_{2}(s)=~{}F_{\pm}(s)=~{}F_{z}(s)$.

Appendix C Spin-spin correlator for a Fermi gas

To calculate the Laplace transform of the spin-spin correlation functions Eqs. 64, 65, 66, 67 and 68 we rewrite each spin $\mathbf{S}$ as

[TABLE]

where $c_{k\sigma}$ are the electron operators and $\bm{\tau}=({\tau^{x},\tau^{y},\tau^{z}})$ is the vector of Pauli matrices. Let us focus first on the zero-field case, $b_{z}^{I}=0$, in which from the $SU(2)$ symmetry of the Fermi gas it follows that $F_{1}(s)=~{}F_{2}(s)=~{}F_{\pm}(s)=~{}F_{z}(s)$. This allows us to restrict the calculation to $F_{z}(s)$, which for the Fermi gas is given by

[TABLE]

In this expression, the sum runs over all $k,k^{\prime}$ within the first Brioullin zone with the $d$-dimensional unit cell volume $a^{2d}$, and $\sigma$ is the spin index. To evaluate the $k$ summations, we introduce the density of states $\nu(\epsilon)=~{}a^{d}\lvert d^{d}k/d\epsilon\rvert/(2\pi)^{d}$ and integrate over $\epsilon$ instead. To capture the important physics about the Fermi energy $E_{F}$ without being troubled by the nonuniversal high-energy contributions near the band edges, we write the density of states as $\nu(\epsilon)=\nu_{0}\exp{\left(-\lvert\epsilon\rvert/\xi_{0}\right)}$, where we choose to set the zero of $\epsilon$ to the Fermi energy and where $\xi_{0}\sim E_{F}$ is a high-energy cutoff. This approximation maintains a constant $\nu(\epsilon)=\nu_{0}$ near the Fermi surface but the exponential provides a finite bandwidth $\xi_{0}$ without introducing high-energy artifacts. As a consequence, all nonuniversal (very short time) behavior that depends on the structure of the entire band will be absorbed in the $\xi_{0}$, and the $\xi_{0}$-independent part of the response function represents the universal low-energy physics. Further corrections of order $\xi_{0}^{-1}$ will be neglected. It is then convenient to introduce the function

[TABLE]

where $f(\epsilon)$ is the Fermi-Dirac distribution. For the last identity, we have used that when $\epsilon$ is measured off the Fermi surface that $\nu(\epsilon)=\nu(-\epsilon)$ and $f(-\epsilon)=1-f(\epsilon)$. With $f(\epsilon)=[1+\exp(\epsilon/k_{B}T)]^{-1}$, the energy integration leads to

[TABLE]

The cutoff $\xi_{0}$ regularizes the short time divergence of the remaining time integral associated with the nonuniversal high-energy sector, and the $\xi_{0}$-independent contribution encodes the universal low energy physics. As mentioned, we neglect in the remaining time integral all orders $\xi_{0}^{-1}$ or smaller. We then obtain

[TABLE]

where $\psi(z)$ is the digamma function, and $\xi$ is a slight renormalization of $\xi_{0}$ absorbing the constant $\psi(1)/2$ in the logarithm. For the response function $F_{z}$, this $F$ contributes the first term of the commutator of Eq. 68. The second term has a reverse time evolution in the expectation values, leading to $e^{-i\epsilon t}$ instead of $e^{i\epsilon t}$ in Eq. 76. As a result, the integral evaluates to the same $F(s)$ as in Eq. 78 but with a conjugation in the cutoff $i\xi_{0}\to-i\xi_{0}$. To avoid any ambiguity, we define the function

[TABLE]

In terms of $F$ and $G$, we then have

[TABLE]

where $\alpha=(\nu_{0}A)^{2}$ is the small dimensionless coupling parameter.

Notice that all correlation functions are of the form $F+G$ with various arguments (see below), and therefore the nonuniversal terms $\pm i\xi_{0}$ always cancel and the only dependence on the cutoff remains in the logarithms. This allows us to drop the $\pm i\xi_{0}$ and use henceforth

[TABLE]

Since for the free Fermi gas at zero magnetic field $F_{z}=~{}F_{1}=~{}F_{2}=~{}F_{\pm}$, Eq. 80 provides all spin-spin correlators. For a finite field $b_{z}^{I}\neq 0$, we have seen through Eqs. 69, 70, 71, 72 and 73 that we can still express the correlators in terms of the $SU(2)$ symmetric laboratory frame, which means in terms of $F(s)$ and $G(s)$ but with arguments $s$ that are shifted by $\pm ib_{z}^{I}$ as follows:

[TABLE]

Notice that only $F_{z}$ remains unchanged from the zero-field expression. Notice furthermore that because of the different cutoffs $\pm i\xi$ in $F$ and $G$ the two correlators $F_{1}(s)$ and $F_{2}(s)$ are different. Because of this, we obtain the correct equilibrium magnetization $\rho^{\text{eq}}_{z}$ derived in Eq. 89.

Appendix D Markov poles

The full time evolution for the reduced density matrix is given by Eqs. 14 and 15, which are evaluated through the sum over the residues at the poles of the integrands. One of the poles of $\rho_{z}$ is always trivial, $s=0$, and leads to the long time equilibrium value. The remaining poles can be split in one Markov pole and an infinite number of non-Markovian poles. In this appendix we evaluate the Markov pole for the different components of the reduced density matrix.

For the longitudinal component $\rho_{z}(s)$, we split off the $s=0$ pole, allowing us to write

[TABLE]

The first term produces in Eq. 14 a time-independent contribution and hence the equilibrium value $\rho_{z}^{\text{eq}}$ obtained in the long time limit $t\to\infty$. Using the results of Appendix C, we obtain

[TABLE]

corresponding indeed to the conventional magnetization of a magnetic moment in a magnetic field.

The second pole of $\rho_{z}$ responsible for the dynamics is located at

[TABLE]

where for the second equality we have used the identity $2\text{Im}\psi(1+ix)=-1/x+\pi/\tanh(\pi x)$ (DLMF, , Sec. 5.4.18). For the transverse components $\rho_{\pm}(t)$, the Markov poles sit at

[TABLE]

where $\xi^{\prime}$ is a slightly shifted cutoff that absorbs a constant $\psi(1)$ contribution. The residues for these poles $s_{z,\pm}^{\text{M}}$ are

[TABLE]

The real part of the pole $s_{z,\pm}^{\text{M}}$ defines the relaxation time $T_{1}$ for the longitudinal component $\rho_{z}$ of Eq. 14 and the decoherence time $T_{2}$ for the transverse component $\rho_{\pm}$ of Eq. 15 for the Markovian part of the decay (reintroducing $\hbar$ in the following expressions):

[TABLE]

These results are reported in Eqs. 29 and 30 in the main text. Since $s_{\pm}^{\text{M}}$ is complex, it causes oscillations in addition to the exponential decay, the spin’s precession in the magnetic field. Writing $s_{\pm}=-\frac{1}{T_{2}}+i\omega_{\pm}$, we capture this precession in the phase factors

[TABLE]

Combining all the results leads to the Markovian decay of the reduced density matrix stated in Eqs. 32 and 33.

Appendix E Non-Markovian poles

Apart from the isolated Markov pole, we expect to find an array of poles close to the poles of the digamma function $\psi$. Using the series expansion for $\psi(1+z)$,

[TABLE]

we see that the poles lie at negative integers $z=-n$. Again, we solve $s+iF_{1}(s)+iF_{2}(s)=0$ up to $O(\alpha)$. For $b_{z}^{I}=0$, the digamma function to investigate is $\psi(1+~{}s/2k_{B}T)$ and therefore we seek poles of the form $s_{n}=-2\pi k_{B}Tn+p$, where $n\geq 1$ is an integer and $p$ is a small correction. Using the expansion of Eq. 97, we find to order $\alpha$ that $p$ is independent of $n$ and given by

[TABLE]

The residue corresponding to $s_{n}$ is, to order $\alpha$,

[TABLE]

Here and for the remainder of this appendix we set the constant amplitudes $[\rho_{z}(t=0)-\rho_{z}^{\text{eq}}]=\rho_{\pm}(t=0)=1$ and reintroduce them only in the final results.

The sum over all residues provides the non-Markovian time evolution. Since at $b_{z}^{I}=0$ the evolution of $\rho_{z}(t)$ and $\rho_{\pm}(t)$ through Eqs. 14 and 15 is identical we find that (reintroducing $\hbar$)

[TABLE]

where $\tau_{T}=\hbar/2k_{B}T$ is the thermal time.

For $b_{z}^{I}>0$ we have to investigate in addition combinations with the shifted digamma functions $\psi(1+(s\pm ib_{z}^{I})/2k_{B}T)$. Each pole splits now into two and we must distinguish between the different $b_{z}^{I}$ dependences for $\rho_{z}$ and $\rho_{\pm}$. If we start with $\rho_{z}$, we write the pair of poles near $-2\pi k_{B}Tn$ as $s_{n,r}^{z}=-2\pi k_{B}Tn+p_{r}^{z}$, where $r=\pm$ and

[TABLE]

The two corresponding residues become

[TABLE]

The limit $b_{z}^{I}=0$ is correctly recovered since one of the residues then vanishes. The summation over all $n$ provides then the non-Markovian evolution of $\rho_{z}$,

[TABLE]

where we have defined

[TABLE]

Reintroducing now the amplitude $[\rho_{z}(t=0)-\rho_{z}^{\text{eq}}]$, we obtain the result

[TABLE]

For $\rho_{\pm}$, the only difference is that the arguments $s+ib_{z}^{I}$ and $s-ib_{z}^{I}$ of the digamma function terms are replaced by $s$ and $s\mp ib_{z}^{I}$. This is exactly the same situation under the change of $b_{z}^{I}\to b_{z}^{I}/2$ and a shift of $s$ by $\mp ib_{z}^{I}/2$. Consequently, by exactly the same analysis,

[TABLE]

The behavior of $h(t,b_{z}^{I})$ on time scales relevant to the non-Markovian decay is shown in Fig. 9. For $b_{z}^{I}<~{}4\alpha\pi k_{B}T$, it contributes to the exponential decay but for larger magnetic fields it becomes oscillating with the frequency $\frac{1}{\hbar}{\sqrt{(b_{z}^{I})^{2}-(4\alpha\pi k_{B}T)^{2}}}$, which at large fields or low $T$ saturates at the nuclear cyclotron frequency $b_{z}^{I}/\hbar$.