**5. Discussion and Conclusions**

In this paper, we reviewed studies which go beyond the conventional adiabatic approximation for open quantum systems to transfer energy quanta and electron spins with using the full counting statistics, which could provide conditions to show quicker transport. We considered a setup consisting of a two-level system representing an anharmonic junction or a quantum dot and its environment(s) representing a canonical or grand canonical ensemble of the energy quanta and the electron to be transferred. We needed to take into account relaxation phenomena in discussing the transfer. In this case, the adiabatic approximation corresponded to the situation where the relevant system such as the two-level system approaches its stationary state faster than the period of modulation, that is, *τ*<sup>−</sup><sup>1</sup> *r* Ω with *τr* the relaxation time of the two-level system and Ω the modulation frequency. Because the relaxation time is finite, the condition for which the adiabatic approximation is valid corresponds to the much longer period of the modulation than *τr*. This means that we can analyze systematic features including adiabatic as well as nonadiabatic features by changing the ratio of the modulation period and *τr*. To clarify the relationship between modulation period and *τr*, we discretized the external modulation thereby permitting a systematic analysis of the ratio by changing each interval while retaining the validity of the Born–Markov approximation. For energy quanta pumping, we showed that the nonadiabatic effect contributes a new term to the formula for the pumped quantity under the adiabatic approximation. For spin pumping, we showed that adiabaticity made no contribution but nonadiabaticity is essential. Comparing these features, we showed that the adiabatic contribution can vanish when the stationary state does not depend on the external modulation as for spin pumping. This means that we need to pay attention to the feature of the stationary state in using the adiabatic approximation in describing relaxation phenomena. (We would draw the reader's attention to the differences in the meaning of nonadiabaticity which has been used in the electron charge pumping by modulation of single gate voltage [61,62].)

With the same setup, the role of nonadiabaticity in pumping phenomena involving energy quanta was discussed more extensively under continuous modulation [63] where the relaxation of the two-level system is treated within the Born–Markovian approximation. In recent work of the present authors on the role of the non-Markovian effect on spin pumping phenomena [60], we found that a nonzero impetus of the dynamics of the pumped quantity under the Born–Markovian approximation shows an unphysical effect, especially for higher modulation frequencies or for the short time regimes. Because the instantaneous impetus contributes strongly under continuous modulation, including the non-Markovian effect would also be necessary in pumping phenomena of energy quanta, especially in evaluating the feature under continuous modulation. This situation remains an open problem. In addition, we described in this work the relaxation process with ordered cumulants of up to second order in the system–environment interaction. An extension to higher orders of cumulants is necessary if we are to discuss relaxation phenomena under strong system–environment interactions. The inclusion of cumulants up to infinite order within the Markovian approximation has been discussed for spin pumping phenomena within the linear response regime using the Green functions [19]. To discuss the

non-Markovian effect, it would be necessary to include higher-order cumulants, a topic that remains for a future study.

We can find recent extensions of the treatments with full counting statistics into the strong system-environment coupling for heat transfer [64] and electron pumping [65]. The essential idea to go beyond the weak coupling is to use the similarity (unitary) transformations: the polaron transformation (the reaction coordinate mapping) is used in the former (latter) studies, respectively. As mentioned in the general formalism of FCS, it should be noted that we need careful treatments on the joint probability, Equation (2), when we use the similarity (unitary) transformations on the time evolution operator. The transformation of the projective measurement is also necessary to recover the original joint probability (See Reference [45]). It might be necessary to compare the dynamics of transported quantity with and without the transformation of the projective measurement.

Since the treatment of FCS to discuss the nonadiabatic effects on quantum pumping is general, we can apply it to many other cases: One of the most interesting issues is to study the non-adiabatic treatment on the combined effect caused by multiple external parameters such as in References [66–68] where adiabatic transport of charge and/or heat is discussed under time-dependent potential and two reservoirs with biased potentials. We can find other issues to remove the adiabatic approximation in spin pumping via a quantum dot between reservoirs with biased chemical potentials [69] and in the quantum transport and/or quantum pumping under dynamical motion of quantum dot [70] based on the recent developments of experimental techniques on microelectromechanical systems [71]. Further, it would be interesting to discuss the non-adiabatic effect on ac-driven electron systems coupled to multiple reservoirs at finite temperature whose adiabatic treatment is discussed in Reference [72]. We expect that these treatments could provide insights to find new applications, such as the design of nanomachines and understanding of the quantum thermodynamics, as well as quicker transport.

**Author Contributions:** The authors contributed equally to this work.

**Funding:** This work was supported by a Grant-in-Aid for Challenging Exploratory Research (Grant No. 16K13853), partially supported by a Grant-in-Aid for Scientific Research on Innovative Areas, Science of Hybrid Quantum Systems (Grant No. 18H04290), and the Open Collaborative Research Program at National Institute of Informatics Japan (FY2018).

**Acknowledgments:** The authors thank Gen Tatara for valuable discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A. Derivation of Quantum Master Equation for FCS**

When we consider the FCS [44], the density operator *W*(*λ*)(*t*) evolves in time in accordance with the modified Liouville–von Neumann equation

$$
\dot{\mathcal{W}}^{(\lambda)}(t) = -i\mathcal{L}^{(\lambda)}\mathcal{W}^{(\lambda)}(t),\tag{A1}
$$

where L(*λ*) is the Liouville operator defined as L(*λ*)*A* = 1 *h*¯ [*<sup>H</sup>λ<sup>A</sup>* − *AH*−*<sup>λ</sup>*] ≡ 1 *h*¯ [*<sup>H</sup>*, *<sup>A</sup>*]*λ* for arbitrary operator *A* with *Hλ* = *e*(*i*/2)*<sup>λ</sup>HE He*−(*i*/2)*<sup>λ</sup>HE* = *H*0 + *HSE*,*λ*. With these relations, L is divided into

$$
\mathcal{L}^{(\lambda)} = \mathcal{L}\_0^{(\lambda)} + \mathcal{L}'^{(\lambda)} \tag{A2}
$$

where

$$
\mathcal{L}\_0^{(\lambda)} A = \frac{1}{\hbar} [H\_0 A - A H\_0] \quad , \ \mathcal{L}^{\prime(\lambda)} A = \frac{1}{\hbar} [H\_{SE,\lambda} A - A H\_{SE,-\lambda}].\tag{A3}
$$

We eliminate the variables of the environment using a projection operator P, which satisfies the idempotent relation, P<sup>2</sup> = P. We also introduce a complementary operator Q ≡ 1 − P. Denoting the relevant and irrelevant parts of the time evolution operator as [52]

$$\mathbf{x}(t) \equiv \mathcal{P}e^{-i\mathcal{L}^{(\lambda)}t}, \quad \mathbf{y}(t) \equiv \mathcal{Q}e^{-i\mathcal{L}^{(\lambda)}t},\tag{A4}$$

with an initial time *t*0 = 0, we obtain

$$\frac{d}{dt}\mathbf{x}(t) = \mathcal{P}(-i\mathcal{L}^{(\lambda)})\mathbf{x}(t) + \mathcal{P}(-i\mathcal{L}^{(\lambda)})y(t) \,, \tag{A5}$$

and

$$\frac{d}{dt}y(t) = \mathcal{Q}(-i\mathcal{L}^{(\lambda)})x(t) + \mathcal{Q}(-i\mathcal{L}^{(\lambda)})y(t) \,. \tag{A6}$$

The formal solution of Equation (A6) is given by

$$y(t) \quad = \int\_0^t e^{-\underline{Q}i\mathcal{L}^{(\lambda)}(t-\tau)} \mathcal{Q}(-i\mathcal{L}^{(\lambda)}) x(\tau) d\tau + e^{-\underline{Q}i\mathcal{L}^{(\lambda)}t} \mathcal{Q}.\tag{A7}$$

Using *x*(*τ*) = <sup>P</sup>*ei*L(*λ*)(*<sup>t</sup>*−*<sup>τ</sup>*)*e*<sup>−</sup>*i*L(*λ*)*<sup>t</sup>* = <sup>P</sup>*ei*L(*λ*)(*<sup>t</sup>*−*<sup>τ</sup>*)(*x*(*t*) + *y*(*t*)) and

$$\theta(t) = 1 - \int\_0^t e^{-\mathcal{Q}i\mathcal{L}^{(\lambda)}\tau} \mathcal{Q}(-i\mathcal{L}^{(\lambda)}) \mathcal{P}e^{i\mathcal{L}^{(\lambda)}\tau} d\tau \equiv 1 - \sigma(t),\tag{A8}$$

we rewrite the formal solution of *y*(*t*) in the form

$$\mathbf{y}(t) = \theta(t)^{-1}((1-\theta(t))\mathbf{x}(t) + \mathbf{e}^{-\underline{\Omega}i\mathcal{L}^{(\lambda)}t}\mathbf{Q}).\tag{A9}$$

By substituting Equation (A9) into Equation (A5), we obtain

$$\frac{d}{dt}\mathbf{x}(t) = \mathcal{P}(-i\mathcal{L}^{(\lambda)})\theta(t)^{-1}\mathbf{x}(t) + \mathcal{P}(-i\mathcal{L}^{(\lambda)})\theta(t)^{-1}e^{\mathcal{Q}(-i\mathcal{L}^{(\lambda)})t}\mathcal{Q},\tag{A10}$$

which holds for arbitrary projection operator and initial condition. Using the relation *<sup>θ</sup>*(*t*)−<sup>1</sup> = ∑∞*<sup>n</sup>*=<sup>0</sup> *<sup>σ</sup>*(*t*)*<sup>n</sup>*, the first term on the right-hand side of Equation (A10) is rewritten as

$$\mathcal{P}(-i\mathcal{L}^{(\lambda)})\theta(t)^{-1}\mathbf{x}(t) = \mathcal{P}(-i\mathcal{L}^{(\lambda)})\mathbf{x}(t) + \mathcal{P}(-i\mathcal{L}^{(\lambda)})\sigma(t)\mathbf{x}(t) + \cdots \tag{A11}$$

To pick out the lower order of L(*λ*), we use the relation

$$e^{-\mathcal{Q}i\mathcal{L}^{(\lambda)}t}\mathcal{Q} = e^{-i\mathcal{L}\_0^{(\lambda)}t}\mathcal{Q}T\_+\exp[\int\_0^t dt'e^{i\mathcal{L}\_0^{(\lambda)}t'}\mathcal{Q}(-i\mathcal{L}'^{(\lambda)})\mathcal{Q}e^{-i\mathcal{L}\_0^{(\lambda)}t'}],\tag{A12}$$

and PL(*λ*) 0 = L(*λ*) 0 P, which gives

$$\begin{array}{rcl} \mathcal{P}(-i\mathcal{L}^{(\lambda)})\theta(t)^{-1}\mathbf{x}(t) &=& \mathcal{P}(-i\mathcal{L}^{(\lambda)})\mathbf{x}(t) + \mathcal{P}(-i\mathcal{L}^{(\lambda)})\int\_{0}^{t}e^{-i\mathcal{L}^{(\lambda)}\_{0}\tau}\mathcal{Q}(-i\mathcal{L}^{(\lambda)})\mathcal{P}e^{i\mathcal{L}^{(\lambda)}\_{0}\tau}d\tau\mathbf{x}(t) + \cdots \\ &=& \mathcal{P}(-i\mathcal{L}^{(\lambda)})\mathbf{x}(t) + \mathcal{P}(-i\mathcal{L}^{(\lambda)})\int\_{0}^{t}\mathcal{Q}(-i\mathcal{L}^{(\lambda)}(-\tau))\mathcal{P}d\tau\mathbf{x}(t) + \cdots \end{array} \tag{A13}$$

where we have used the definition

$$
\hat{\mathcal{L}}\_1^{(\lambda)}(t) = e^{i\mathcal{L}\_0^{(\lambda)}t} \mathcal{L}'^{(\lambda)} e^{-i\mathcal{L}\_0^{(\lambda)}t}.\tag{A14}
$$

When we multiply the initial condition *W*(*λ*)(0) by the right-hand side of Equation (A8), we obtain the TCL equation for reduced density operator under FCS.

Let us consider a projection operator P = *ρ*ETrE where TrE refers to a trace operation over the environment. When we multiply the initial condition of the density operator of the total system, *<sup>W</sup>*(*<sup>λ</sup>*, 0) from the right by *xS*(*t*) in Equation (A2), we obtain

$$\mathbf{x}(t)\mathcal{W}^{(\lambda)}(0) = \rho\_\mathcal{\mathbb{E}}\mathrm{Tr}\_\mathcal{\mathbb{E}}\mathcal{W}^{(\lambda)}(t) \tag{A15}$$

Defining TrE*W*(*λ*)(*t*) ≡ *<sup>ρ</sup>*(*λ*)(*t*), we obtain the TCL equation for the reduced density operator under FCS,

$$\frac{d}{dt}\rho^{(\lambda)}(t) = \text{Tr}\_{\mathbb{E}}[(-i\mathcal{L})\rho^{(\lambda)}(t)] + \zeta^{(\lambda)}(t)\rho^{(\lambda)}(t) + \psi^{(\lambda)}(t),\tag{A16}$$

with

$$\mathcal{Z}^{(\lambda)}(t) \equiv \sum\_{n=2}^{\infty} \mathcal{Z}\_n^{(\lambda)}(t),\tag{A17}$$

$$\psi^{(\lambda)}(t) \equiv \text{Tr}\_{\mathbb{E}}[(-i\mathcal{L})\theta(t)^{-1}e^{\mathbb{Q}(-i\mathcal{L})t}\mathcal{Q}\mathcal{W}^{(\lambda)}(0)].\tag{A18}$$

Equation (9) with (10) is obtained by taking the lower term of *ξ*(*λ*)(*t*) up to second order in L and replacing *ξ*(*λ*) 2(*t*) with *ξ*(*λ*)(*t*),

$$\xi^{(\lambda)}(t)\rho(\lambda,t) = \int\_0^t \text{Tr}\_{\mathbb{E}}[(-i\mathcal{L}')\mathcal{Q}(-i\mathcal{L}'(-\tau))\rho\_{\mathbb{E}}\rho(\lambda,t)]d\tau. \tag{A19}$$

with the factorized initial condition of the total system as *W*(*λ*)(0) = *<sup>ρ</sup>*E*ρ<sup>λ</sup>*(0), which makes the third term on the right-hand side of Equation (A16) vanish. With the assumption TrE*HSE*,*<sup>λ</sup>* = 0, we have

$$\begin{array}{lcl} \mathfrak{F}^{(\lambda)}(t)\rho^{(\lambda)}(t) & \\ = (\frac{i}{\hbar})^2 \text{Tr}\_{\text{E}}[H\_{SE}[H\_{SE}(-\tau), \rho\_{\text{E}}\rho^{(\lambda)}(t)]\_{\lambda}]\_{\lambda}]\_{\lambda} \\ = (\frac{i}{\hbar})^2 \text{Tr}\_{\text{E}}[H\_{SE,\lambda}H\_{SE,\lambda}(-\tau)\rho\_{\text{E}}\rho^{(\lambda)}(t) - H\_{SE,\lambda}\rho\_{\text{E}}\rho^{(\lambda)}(t)H\_{SE,-\lambda}(-\tau)}{-H\_{SE,\lambda}(-\tau)\rho\_{\text{E}}\rho^{(\lambda)}(t)H\_{SE,-\lambda} + \rho\_{\text{E}}\rho^{(\lambda)}(t)H\_{SE,-\lambda}(-\tau)H\_{SE,-\lambda}]\_{\lambda}} \\ \end{array} \tag{A20}$$

which coincides with Equation (110) in [44].

### **Appendix B. Connection between Two Formalisms of the FCS**

We next explain the relationship between the two formalisms of the FCS provided by Esposito et al. in Reference [44] and by Sinitsyn and Nemenman in Reference [27]. To establish the connection between the two formalisms, we consider the number of quanta *N* as the observable to be measured to identify the number of quantum transfers from S to E, as stipulated in Reference [27]. For a large class of open quantum systems, the system–environment interaction is described by an interaction Hamiltonian of the form *H*SE = *V*+ + *V*− ≡ *A* ⊗ *B*† + *A*† ⊗ *B*, where *B*† and *B* are creation and annihilation operators of a quanta in the environment, respectively, and *A* is either a Hermitian or non-Hermitian operator acting on the relevant system. *V*+ ≡ *A* ⊗ *B*† or *V*− = *V*†+ describes transfer of a quanta from S to E or from E to S. In this instance, the number operator of the quanta is given by *N* = *B*†*B*. The generalized Liouvillian is expressed as L(*λ*) = L0 + <sup>L</sup>+*eiλ*/2 + L−*e*<sup>−</sup>*iλ*/2, where <sup>L</sup>0*W*(*λ*) ≡ *h*¯ <sup>−</sup><sup>1</sup>(*<sup>H</sup>*0*W*(*λ*) − *<sup>W</sup>*(*λ*)*H*0) is the unperturbed Liouvillian, <sup>L</sup>+*W*(*λ*) ≡ *h*¯ <sup>−</sup><sup>1</sup>(*<sup>V</sup>*+*W*(*λ*) − *W*(*λ*)*V*†+) and L−*W*(*λ*) ≡ *h*¯ <sup>−</sup><sup>1</sup>(*<sup>V</sup>*−*W*(*λ*) − *W*(*λ*)*V*† −) are Liouvillians describing the transfer of a quanta from S to E and from E to S, respectively. Using the formal solution of Equation (8) with the given Liouvillian in Equation (6), we obtain a formal expression for the moment generating function

$$Z(\lambda) = \sum\_{n=-\infty}^{\infty} P\_n e^{in\lambda},\tag{A21}$$

with *Pn* the probability of having *n* net transitions from S to E, e.g.,

$$\begin{split} P\_{1} &= \text{Tr}\Big[\Big(\int\_{t\_{i}}^{t} dt\_{1} \mathbb{L}l\_{0}(t\_{i},t\_{1})(-i\mathcal{L}\_{+})\mathbb{L}l\_{0}(t\_{1},t\_{i}) + \int\_{t\_{i}}^{t} dt\_{3} \int\_{t\_{i}}^{t\_{3}} dt\_{2} \int\_{t\_{i}}^{t\_{2}} dt\_{1} \mathbb{L}l\_{0}(t\_{i},t\_{3})(-i\mathcal{L}\_{+}) \\ & \times \mathbb{L}l\_{0}(t\_{3},t\_{2})(-i\mathcal{L}\_{-})\mathbb{L}l\_{0}(t\_{2},t\_{1})(-i\mathcal{L}\_{+})\mathbb{L}l\_{0}(t\_{1},t\_{i})\cdots \Big)\bar{\mathbb{W}}(t\_{i})\Big], \end{split} \tag{A22}$$

where *<sup>U</sup>*0(*<sup>t</sup>*1, *<sup>t</sup>*2) = exp[−*i*L0(*<sup>t</sup>*<sup>1</sup> − *<sup>t</sup>*2)] is the unperturbed time evolution operator. The expression for the moment generating function corresponds to Equation (9) in Reference [27].

### **Appendix C. An Analytical Proof of the Absence of Adiabatic Contribution to the Spin Pumping**

Following the procedure by Sinitsyn and Nemenman in Reference [27], we divide the cycle of precession into intervals *δt*, which correspond to the step-like changes of *M*. For the step-like precession, the density matrix of the dot during *ti* ≤ *t* ≤ *ti*+<sup>1</sup> is given by

$$|\rho\_i^{(0)}(t)\rangle = e^{\Xi\_i^{(0)}(t-t\_i)} \prod\_{j=1}^{i-1} e^{\Xi\_j^{(0)}\delta t} |\rho\_1^{(0)}(0)\rangle,\tag{A23}$$

where we denote the density matrix and the generator in the *i*th interval as |*ρ*(0) *i* (*t*) and Ξ(0) *i* , respectively, and |*ρ*(0) 1 (0) is the initial condition in the first interval. Taking the spectral decomposition of Ξ(0) *i* and assuming a quick approach of the system to its steady state in each interval in Equation (A23), as in Reference [27], the terms remaining in the decomposition are those that contain the steady state |*u*00(*ti*) satisfying Ξ(0) *i* |*u*00(*ti*) = 0. Evaluating the density matrix up to first order in *δt*, we find the first-order term vanishes, implying the invariance of the steady state populations of the dot under the step-like rotation of *φ* in our model (see Appendix D in Reference [58]). Thus, we find that the density matrix at time *t* under the adiabatic limit can be approximated by the steady state as |*ρ*(0) 0 (*t*) ≈ |*u*(0) 0 (*ti*). For this steady state |*u*(0) 0 (*ti*) given by Equation (A25), we also find that there is no electron transfer between dot and lead; specifically, we find that

$$<\langle \Delta n\_{\sigma,i} \rangle \approx \int\_{t\_i}^{t\_{i+1}} dt' \langle \langle 1 | \left[ \frac{\partial \Xi\_i^{(\lambda\_\sigma)}}{\partial (i \lambda\_\sigma)} \right]\_{\lambda\_\sigma = 0} | u\_0^{(0)}(t\_i) \rangle \rangle = 0,\tag{A24}$$

indicating that there is no net electron transfer in the interval and hence the generated spin current represented by Equation (36) is totally absent in the adiabatic limit. We therefore need to include the nonadiabatic effect to generate a finite spin current. The result is in marked contrast to the previous example of the energy pumping, where the adiabatic contribution to the energy pumping G*<sup>ν</sup>ad*is finite.

### **Appendix D. Steady State of the Minimum Model**

The steady state |*u*(*λσ*=<sup>0</sup>) 0 (*ti*), satisfying Ξ(*λσ*=<sup>0</sup>) *i* |*u*(*λσ*=<sup>0</sup>) 0 (*ti*) = 0 is analytically obtained using a graphical method discussed in Reference [73]. In Appendix C of Reference [58], we provide a detailed derivation of the steady state. For use in the present paper, here we simply present the result.

The dynamics of the populations described by the TCL master equation is closed for the six components of the reduced density matrix, *ρ*00(*t*) ≡ 0, <sup>0</sup>|*ρ*(*λσ*=<sup>0</sup>)(*t*)|0, <sup>0</sup>, *ρ*01(*t*) ≡ 0, <sup>1</sup>|*ρ*(*λσ*=<sup>0</sup>)(*t*)|0, <sup>1</sup>, *ρ*0110(*t*) ≡ 0, <sup>1</sup>|*ρ*(*λσ*=<sup>0</sup>)(*t*)|1, <sup>0</sup>, *ρ*1001(*t*) ≡ 1, <sup>0</sup>|*ρ*(*λσ*=<sup>0</sup>)(*t*)|0, <sup>1</sup>, *ρ*10(*t*) ≡ 1, <sup>0</sup>|*ρ*(*λσ*=<sup>0</sup>)(*t*)|1, <sup>0</sup>, and *ρ*11(*t*) ≡ 1, <sup>1</sup>|*ρ*(*λσ*=<sup>0</sup>)(*t*)|1, <sup>1</sup>. By arranging these connected components as |*ρ*(*λσ*=<sup>0</sup>)(*t*) = [*ρ*00(*t*), *ρ*01(*t*), *ρ*0110(*t*), *ρ*1001(*t*), *ρ*10(*t*), *<sup>ρ</sup>*11(*t*)]t, where [··· ]t denotes transposition, an analytic expression of its steady state obtains,

$$
\left|u\_{0}^{(\lambda\_{\ell}=0)}(t\_{i})\right\rangle = \begin{pmatrix} f^{-}(\boldsymbol{\varepsilon}\_{\uparrow})f^{-}(\boldsymbol{\varepsilon}\_{\downarrow}) \\ \cos^{2}\frac{\theta}{2}f^{+}(\boldsymbol{\varepsilon}\_{\uparrow})f^{-}(\boldsymbol{\varepsilon}\_{\downarrow}) + \sin^{2}\frac{\theta}{2}f^{-}(\boldsymbol{\varepsilon}\_{\uparrow})f^{+}(\boldsymbol{\varepsilon}\_{\downarrow}) \\ e^{+i\phi\_{\uparrow}}\cos\frac{\theta}{2}\sin\frac{\theta}{2}[f^{+}(\boldsymbol{\varepsilon}\_{\uparrow})f^{-}(\boldsymbol{\varepsilon}\_{\downarrow}) - f^{-}(\boldsymbol{\varepsilon}\_{\uparrow})f^{+}(\boldsymbol{\varepsilon}\_{\downarrow})] \\ e^{-i\phi\_{\uparrow}}\cos\frac{\theta}{2}\sin\frac{\theta}{2}[f^{+}(\boldsymbol{\varepsilon}\_{\uparrow})f^{-}(\boldsymbol{\varepsilon}\_{\downarrow}) - f^{-}(\boldsymbol{\varepsilon}\_{\uparrow})f^{+}(\boldsymbol{\varepsilon}\_{\downarrow})] \\ \sin^{2}\frac{\theta}{2}f^{+}(\boldsymbol{\varepsilon}\_{\uparrow})f^{-}(\boldsymbol{\varepsilon}\_{\downarrow}) + \cos^{2}\frac{\theta}{2}f^{-}(\boldsymbol{\varepsilon}\_{\uparrow})f^{+}(\boldsymbol{\varepsilon}\_{\downarrow}) \\ f^{+}(\boldsymbol{\varepsilon}\_{\uparrow})f^{+}(\boldsymbol{\varepsilon}\_{\downarrow}) \end{pmatrix},\tag{A25}$$

where *f* +( *k*) ≡ Trl[*c*†*σ*,*kcσ*,*k<sup>ρ</sup>*eql ], *f* <sup>−</sup>( *k*) ≡ Trl[*<sup>c</sup>σ*,*kc*†*σ*,*k<sup>ρ</sup>*eql ], ↑ ≡ d − *M* and ↓ ≡ d + *M*. From the expression, we find that the steady-state values of the populations *ρ*00, *ρ*01, *ρ*10 and *ρ*11 are independent of angle *φ*. Thus, the steady state populations remain unchanged by changing *φ*.
