Accessing the bath information in open quantum systems with the stochastic c-number Langevin equation method

TL;DR

This paper introduces a stochastic c-number Langevin equation method that enables direct access to bath information in open quantum systems, providing a more comprehensive understanding of system-bath interactions.

Contribution

The work develops a novel stochastic approach that allows simultaneous analysis of system and bath dynamics, extending traditional methods that only focus on the system.

Findings

Successfully applied to models like pure dephasing and spin-boson.

Able to extract bath information alongside system dynamics.

Demonstrated effectiveness in quantum dot-phonon interactions.

Abstract

In traditional open quantum systems, the baths are usually traced out so that only the system information is left in the equations of motion. However, recent studies reveal that using only the system degrees of freedom can be insufficient. In this work, we develop a stochastic c-number Langevin equation method which can conveniently access the bath information. In our method, the studied quantities are the expectation values of operators which can contain both system operators and bath operators. The dynamics of the operators of interest is formally divided into separate system and bath parts, with auxiliary stochastic fields. After solving the independent stochastic dynamics of the system part and the bath part, we can recombine them by taking the average over these stochastic fields to obtain the desired quantities. Several applications of the theory are highlighted, including the pure dephasing model, the spin-boson model, and an optically excited quantum dot coupled to a bath of phonons.

Tables

Table 1. Table 1: Numbers of trajectories for different results

Pure dephasing model (low temperature)	$4.5\times {10}^7$
Pure dephasing model (high temperature)	$5\times {10}^6$
Pure dephasing model with control	$4\times {10}^7$
Spin-boson model	$4\times {10}^7$
Spin-boson model with pump	${10}^7$
Quantum dot with optical pump	${10}^6$

Equations

$$H_0$$

$$H_{\mathrm{b}}$$

$$H_{\mathrm{I}}$$

$${\dot{A}}_s(t)$$

$${\dot{a}}_k(t)$$

$$a_k=e^{-i{\omega }_{k}t}{a}_k(0)-i{g}_k{\int }_0^{t}ds{e}^{-i{\omega }_k(t-s)}S(s).$$

$${\dot{A}}_s(t)$$

$$\alpha (t-s)\equiv \sum _k{g}_k{g}_k^{\ast }\exp (-i{\omega }_k(t-s)),$$

$$\hat{\xi }(t)\equiv \sum _k\left[g_k{a}_k(0)\exp (-i{\omega }_{k}t)+g_k^{\ast }{a}_k^{†}(0)\exp (i{\omega }_{k}t)\right],$$

$$B(t)=U^{†}(t)B(0)U(t),$$

$$U(t)$$

$$U_{\mathrm{sys}}(t;z_{1,\tau },z_{2,\tau })$$

$$U_{\mathrm{b}}(t;z_{1,\tau },z_{2,\tau })$$

$$B(t)$$

$$B_{\mathrm{sys}}(t;z)$$

$$B_{\mathrm{b}}(t;z)$$

$$⟨B(t)⟩=\mathrm{Tr}\left\{\mathrm{B}(\mathrm{t}){\rho }_{\mathrm{tot}}(0)\right\}.$$

$$⟨B(t)⟩$$

$$⟨B_{\mathrm{sys}}(t;z)⟩\equiv \mathrm{Tr}\left\{{\mathrm{B}}_{\mathrm{sys}}(\mathrm{t};\mathrm{z}){\rho }_{\mathrm{sys}}(0)\right\},$$

$$⟨B_{\mathrm{b}}(t;z)⟩\equiv \mathrm{Tr}\left\{{\mathrm{B}}_{\mathrm{b}}(\mathrm{t};\mathrm{z}){\rho }_{\mathrm{b}}(0)\right\}.$$

$$B_{\mathrm{sys}}(0;z)$$

$$B_{\mathrm{sys}}(t;z)$$

$$\mathrm{with},$$

$$Y_l(t;z)$$

$$⟨B_{\mathrm{sys}}(t;z)⟩=\sum _l{b}_l⟨Y_l(t;z)⟩.$$

$$\frac{\partial }{\partial t}{Y}_l(t;z)$$

$$D(0;z)$$

$$x_{1,t}$$

$$x_{2,t}$$

$$D(0;x)$$

$$[H_{\mathrm{sys}},Y_l]$$

$$[S_{\mathrm{sys}},Y_l]$$

$$\{S_{\mathrm{sys}},Y_l\}$$

$$\frac{\partial }{\partial t}{Y}_l(t;x)=i\sum _m\left({\mathcal{H}}_{lm}+\frac{x_{1,t}}{\sqrt{2}}{\mathcal{S}}_{lm}^c+\frac{x_{2,t}}{\sqrt{2}}{\mathcal{S}}_{lm}^a\right)Y_m(t;x).$$

Full text

I Introduction

Open quantum systems are usually divided into the most relevant parts (the system) and the secondary parts (the bath). Therefore one can use a partial trace to eliminate the bath. However, recent studies make the division between the system and the bath somewhat unclear. For example, the closed many-body systems do not have bath parts, but can still thermalize mb2 ; mb3 ; mb4 ; mb5 ; mb6 ; mb7 ; mb8 ; mb9 ; mb10 . Large coupling strengths can make the correlation between the system and the bath important in thermodynamics qtd1 ; qtd3 ; qtd4 ; qtd5 ; qtd6 ; qtd7 ; qtd8 ; qtd9 . In photosynthesis, the mixing of the system modes and the bath modes can suppress decoherence qbio1 ; qbio2 ; qbio3 ; qbio4 ; photosyn1 ; photosyn2 ; photosyn3 ; photosyn4 ; photosyn5 ; photosyn5 . Moreover, when the bath is influenced by the pump exerted on the system, tracing out the bath is difficult and usually contains approximations drib1 ; drib2 ; drib3 . Thus a method without performing a partial trace is very desirable.

To retain the potential information lost in tracing out the baths, one approach is to introduce generating functions mebi1 ; mebi2 ; mebi3 . However, it is difficult to find the generating functions in general cases. The Heisenberg-Langevin approach qle1 ; qle2 ; qle3 ; qle4 ; qle5 ; qle6 can be another good choice for such problems without apparent divisions between the systems and the baths. Although the Langevin equations consist of the system parts and the bath parts, no partial trace is included. Therefore, the accessible quantities are not restricted to the system parts qleb1 ; qleb2 ; qleb3 . Nevertheless, the bath parts become the quantum noise terms in the equations, which makes the numerical simulations of the Langevin equations very difficult. Also, in some models, e.g., multi-level systems, the Langevin equations may contain the nonlinear time-non-local terms. Thus the use of the Heisenberg-Langevin approach in open quantum systems has limitations.

In the Schr$\rm\ddot{o}$dinger picture, nonlinear time-non-local terms can be avoided. Therefore, researchers have developed many successful methods, like the quantum state diffusion method Diosi ; qsd1 ; qsd2 ; qsd3 ; qsd4 , polaron transform methods polatran1 ; polatran2 ; polatran3 ; polatran4 ; polatran5 ; polatran6 ; polatran7 ; polatran8 ; polatran9 , the hierarchy equation methods hie1 ; hie2 ; hie3 ; hie4 ; hie5 ; hie5.5 ; hie6 ; hie7 , path integral methods pathi , and the stochastic Liouville equation methods sfd0 ; sfd1 ; sfd2 ; sfd2.5 ; sfd3 ; sfd4 ; sfd5 ; oqs10 . These methods are suitable for different situations, but with various limitations. For example, the quantum state diffusion method can efficiently calculate large systems, but its application is difficult for many models. The polaron transform methods, which can provide compact master equations, usually include perturbations in either the transformations or the derivations of the master equations (which treats only certain parts of the bath coupling non-perturbatively). One can calculate the long-time results with hierarchy equation methods if the correlation functions of the baths are not too complicated. The path integral methods are powerful but complicated, and restricted when adding in other effects to the model. The stochastic methods are not restricted to any bath spectrum or temperature, but at the cost of poor long-time performance. However, most of the recent methods contain a partial trace.

In this article, we combine the Heisenberg-Langevin method and the stochastic Liouville equation to overcome some of the mentioned disadvantages of these recent methods. To avoid the nonlinear time-non-local terms in the equations, we separate the dynamics of the system and the bath with the Hubbard-Stratonovich transformation hstrans (by introducing auxiliary stochastic fields). Then the Langevin equations become linear and time local at the cost of containing classical noise terms. Moreover, we change the stochastic equations to c-number equations by taking expectation values. The final result can be obtained by taking the average over the noise terms. Our method is not based on the perturbation or the Markovian approximation, so it is, in principle, numerically exact. Similar to the stochastic Liouville equation, the stochastic c-number Langevin equations can deal with very complicated bath spectra at any temperature, but may have poor numerical convergence for long time simulations. In addition, our method also has the advantage of the Heisenberg-Langevin method , in that it can conveniently obtain the quantities of the bath.

Our paper is organized as follows. In Sec. II, we derive the stochastic c-number Langevin equations. We first separate the dynamics of the system and the bath by introducing the stochastic noise terms. Then, we solve the stochastic dynamics of the system part and the bath part separately. Finally, we combine the system part and the bath part by taking the average over the noise terms to obtain the desired results. In Sec. III, several numerical examples of two-level systems are calculated. The first example is the pure dephasing model which can be used to check the method. The second example, the spin-boson model, has been calculated with many numerical methods, but the bath quantities are seldom studied. Then a realistic and practical system, a quantum dot system coupled to phonons, is chosen as the third example. Section IV presents our conclusions. We also include four appendices. Appendix A provides the derivation of the stochastic bath evolution operator. In Appendix B, we describe how to generate the bath operators from the stochastic bath evolution operator. In Appendix C, we derive the expectation value of the stochastic identity operator of the bath. In Appendix D, our method is extended to the time-dependent Hamiltonian case.

II The stochastic c-number Langevin equations

For the total system, we consider the Calderia-Leggett type model c-l model , with $H_{\rm tot}=H_{0}+H_{\rm I}$ (setting $\hbar=1$),

[TABLE]

Here, $H_{0}$ is composed of the free Hamiltonian of the system, $H_{\rm sys}$, and the multi-mode bosonic bath $H_{\rm b}$. The coupling Hamiltonian $H_{\rm I}$ describes the coupling between the system and the bosonic bath. In the coupling Hamiltonian, $S$ is a Hermitian operator of the system, and $a_{k}^{\dagger}$ ($a_{k}$) is the bosonic creation (annihilation) operator of the $k$th-mode in the bath.

Solving such an open system with the Heisenberg-Langevin approach is usually difficult. Here, we take the standard generalized quantum Lanvegin equation as an example osbook1 . Consider the Heisenberg equations of an arbitrary system operator $A_{s}(t)$, and the bath operator $a_{k}$, then

[TABLE]

The bath equations can be formally solved,

[TABLE]

Substituting this formal solution for the bath operators in the equation of the system operator $A_{s}(t)$, the Lanvegin equation of $A_{s}(t)$ can be derived:

[TABLE]

where,

[TABLE]

is the zero-temperature correlation function, and,

[TABLE]

is the quantum noise. In this article, noise terms with hat correspond to quantum noise terms of the bath, and noise terms without hat are classical noise terms terms from the auxiliary field. In general, the commutator $[S,A_{s}(t)]$ is some operator $B_{s}(t)$. Thus the last term in Eq. (LABEL:qlefms) has the form $\int_{0}^{t}ds[\alpha(t-s)-\alpha^{*}(t-s)]B_{s}(t)S(s)$, which is both nonlinear and time non-local. Such a term makes it difficult to apply the Heisenberg-Lanvegin approach to general cases. In addition, the quantum noise further increases the difficulty of solving these equations.

II.1 Separating the system and the bath with stochastic noise terms

To avoid dealing with the nonlinear terms directly, we formally separate the dynamics of the system and the bath sfd1 ; sfd2 ; sfd2.5 ; sfd3 ; sfd4 ; sfd5 . The basic idea is similar to generating indirect interactions with intermediate fields idcoupling1 ; idcoupling2 ; idcoupling3 ; idcoupling4 . Instead of the original system-bath coupling, we consider the equivalent coupling induced by auxiliary fields with trivial dynamics. This produces stochastic, but uncoupled system dynamics and bath dynamics. By averaging over the noise terms, which resembles tracing out the intermediate field, we can recover the original system-bath coupling.

An operator at time $t$ in the Heisenberg picture $B(t)$ is related to its initial value $B(0)$ by

[TABLE]

with the evolution operator $U(t)\equiv\exp(-iH_{\rm tot}t)$. It is difficult to obtain $U(t)$ for general cases. However, by applying the Hubbard-Stratonovich transformation, the evolution operator $U(t)$ can be divided into two parts,

[TABLE]

where,

[TABLE]

and,

[TABLE]

Here, $U_{\rm sys}(t;z_{1,\tau},z_{2,\tau})$ is the stochastic evolution operator of the system, and $U_{\rm b}(t;z_{1,\tau},z_{2,\tau})$ is the stochastic evolution operator of the bath; $T_{+}$ is the time ordered operator, and the average over the noise terms is denoted by $\mathcal{M}_{z}\{\}$. We first assume $B(0)=B_{\rm sys}(0)\otimes B_{\rm b}(0)$ so that the dynamics of the total system can be divided into two parts. Thus:

[TABLE]

In more general cases, with $B(0)=\sum_{i}B_{i,\rm sys}(0)\otimes B_{i,\rm b}(0)$, we can deal with these terms one by one.

In the master equation method, the bath part is usually traced out to obtain the reduced density matrix $\rho_{\rm sys}(t)$. However, some information of the bath is lost in this procedure. To avoid such a disadvantage, we take the expectation value of the operator, instead of tracing out the bath:

[TABLE]

Now we assume that the system and the bath are factorized at the initial time $\rho_{\rm tot}(0)=\rho_{\rm sys}(0)\otimes\rho_{\rm b}(0)$. Then, the expectation value can also be divided into the system part and the bath part,

[TABLE]

with,

[TABLE]

and,

[TABLE]

According to Eq. (11), we can separately calculate the stochastic expectation value of the system part $\langle B_{\rm sys}(t;z)\rangle$ and the bath part $\langle B_{\rm b}(t;z)\rangle$. Then, the expectation value of the concerned operator $\langle B(t)\rangle$ can be obtained by taking the noise average of the stochastic ones. Note that these noise terms come from the auxiliary fields instead of the bath. The quantum noise induced by the bath is included in $\langle B_{\rm b}(t;z)\rangle$.

II.2 The stochastic dynamics of the system part

We first consider a system of finite dimension. In such a case, a set of basis operators $Y_{l}$ can be found to represent the $B_{\rm sys}(0;z)$,

[TABLE]

The stochastic system operator at time $t$ can also be expressed in the same way

[TABLE]

The expectation value of the system part can then be expressed as

[TABLE]

According to Eq. (16), once all the stochastic expectation values of the basis $Y_{l}(t;z)$ are known, the stochastic expectation value of the system part can be easily calculated. The equation for $Y_{l}(t;z)$ can be obtained by directly taking the time derivative of it, so that

[TABLE]

According to the form of the stochastic evolution operator in Eq. (LABEL:sfd), Eq. (LABEL:efsb) can be written as

[TABLE]

where,

[TABLE]

If we define two complex noise terms:

[TABLE]

and,

[TABLE]

We can express $D(0;z)$ in Eq.(LABEL:ssle1) in a compact form,

[TABLE]

The commutators and anti-commutators can also be expressed by the basis operators $Y_{m}$,

[TABLE]

When the dimension of the system Hilbert space is infinite, our method can also be applied if a set of basis operators which satisfies Eq. (21) can be found. With Eq. (LABEL:ssle1) and Eq. (21), we can obtain the following equation for $Y_{l}(t;x)$

[TABLE]

The equation for the expectation values is straightforward to obtain,

[TABLE]

Equation (LABEL:ssle2) can also be written in a vector form if we define the stochastic expectation value vector as

[TABLE]

where $n$ is the number of the basis operators of the system and $T$ means transpose. The vector form of the equation reads

[TABLE]

The elements of the matrices $\mathcal{H}$, $\mathcal{S}^{c}$ and $\mathcal{S}^{a}$ are, respectively, the terms $\mathcal{H}_{lm}$, $\mathcal{S}^{c}_{lm}$, and $\mathcal{S}^{a}_{lm}$ in Eq. (LABEL:ssle2). With Eq. (16) and Eq. (24), the stochastic expectation value of the system part can be calculated. Note that the choice of the basis operators is arbitrary. For example, we can take $B_{\rm sys}(t;z)$ as the first one of them, and find enough operators to satisfy Eq. (21).

II.3 The Stochastic Dynamics of the Bath Part

Next, we come to the stochastic evolution operator for the bath $U_{\rm b}(t;z_{3,\tau},z_{4,\tau})$. In the following, for simplicity, the stochastic evolution operators for the bath will be written as $U_{\rm Ib}(t;z)$ and $U_{\rm Ib}^{{\dagger}}(t;z)$. Note that the noise terms are different in $U_{\rm Ib}(t;z)$ and $U_{\rm Ib}^{{\dagger}}(t;z)$. We first change to the interaction picture,

[TABLE]

Then the bath part of the stochastic operator becomes $B_{\rm b}(t;z)=U_{\rm Ib}^{{\dagger}}(t;z)B_{\rm Ib}(t)U_{\rm Ib}(t;z)$, where $B_{\rm Ib}(t)\equiv\exp({iH_{\rm b}t})B_{\rm b}(0)\exp({-iH_{\rm b}t})$ is the initial value of the bath part of the stochastic operator in interaction picture. The equation for ${U}_{\rm Ib}(t;z)$ is

[TABLE]

Here the $\hat{\xi}(t)$ is the quantum noise term which appears in the quantum Langevin equation. Equation (26) can be solved by the Magnus expansion mge ; mge2 (also see Appendix A),

[TABLE]

Now we consider the stochastic expectation value of the bath part $\langle B_{\rm b}(t;z)\rangle$. We express the bath part of the operator with a Taylor expansion,

[TABLE]

The stochastic bath part of operator at time $t$ is,

[TABLE]

Conceptually Eq. (29) can be directly calculated from Eq. (27). However, the bath operators can also be generated approximately from the stochastic evolution operator when the time is not too short (see Appendix B). In this way,

[TABLE]

where $\mathcal{A}_{k}(t;z_{3,\tau},z_{4,\tau})$ has the following form,

[TABLE]

The term $\alpha(s_{1}-s_{2})\equiv\sum_{k}g_{k}g_{k}^{*}\exp(-i\omega_{k}(s_{1}-s_{2}))$ in Eq. (31) is the zero-temperature correlation function mentioned in Eq. (LABEL:qlefms). With Eq. (31), Eq. (29) can be converted to

[TABLE]

Here, $I_{\rm b}(t;z)\equiv U_{\rm b}^{{\dagger}}(t;z)U_{\rm b}(t;z)$ is the bath part of the stochastic identity operator. Note that the stochastic identity operator is no longer an identity operator. Also, we abbreviate the $\mathcal{A}_{k}(t;z_{3,\tau},z_{4,\tau})$ [$\mathcal{A}^{*}_{k}(t;z_{1,\tau},z_{2,\tau})$] to $\mathcal{A}_{k}(t;z)$ [$\mathcal{A}^{*}_{k}(t;z)$]. Since the stochastic bath operators can be generated from the bath part of the stochastic identity operator, the next step is to obtain the stochastic expectation value of $I_{\rm b}(t;z)$. If the initial state of the bath is the thermal state, then $\langle I_{\rm b}(t;z)\rangle$ can be exactly calculated sfd5 (also see Appendix C):

[TABLE]

with,

[TABLE]

The $\alpha_{T}(t-s)$ term in Eq. (LABEL:seobu) is the finite temperature correlation function of the bath. The effects of the bath temperature $T$ is described by the parameter $\beta\equiv 1/({k_{\rm B}T})$. With the value of $\langle I_{b}(t;z)\rangle$, the stochastic expectation value of the bath part $\langle B_{\rm b}(t;z)\rangle$ can be calculated from:

[TABLE]

Using the exponential property of $\langle I_{\rm b}(t;z)\rangle$, we can assume that $\langle B_{\rm b}(t;z)\rangle$ has the following form,

[TABLE]

The $f(t,x_{3})$ term is some stochastic function, where $x_{3}$ represents one or more noise terms. The property of $x_{3}$ depends on the form of $B_{\rm b}(t;z)$. If we are only interested in the system operators, then the bath part is just $B_{\rm b}(t;z)=I_{\rm b}(t;z)$.

II.4 Noise average of the stochastic operators

To obtain the expectation value of the operator $B(t)$, we need to eliminate the auxiliary stochastic fields by taking a noise average over the product of the system part and the bath part,

[TABLE]

Direct calculation of Eq. (LABEL:saoo1) is numerically inefficient. We can introduce a transformation of the noise terms to absorb $\langle I(t;z)\rangle$ into the relevant distribution function sfd2 ,

[TABLE]

Notice that this transformation only changes the noise terms in the system part. With the transformation in Eq. (LABEL:tron), Eq. (24) becomes,

[TABLE]

By defining two colored noise terms $\xi_{t}=x_{1,t}^{\prime}+g(t,z)$ and $\eta_{t}=ix_{2,t}$, we obtain the following stochastic equation:

[TABLE]

Here, $\theta(t-s)$ is the step function, which is 1 for $t>s$ and 0 for $t<s$; the value for $t=s$ is not important as Im$(\alpha_{T}(0))=0$. Each line of Eq. (LABEL:ssle4) is a stochastic c-number Langevin equation of a stochastic basis operator $Y_{m}(t;\xi,\eta)$. Note that the noise terms in Eq. (LABEL:ssle4) represent classical noise terms, therefore different from the quantum noise terms. By taking the noise average of these expectation values of the stochastic basis operators, we obtain the expectation values of the basis operators $\langle Y_{m}(t)\rangle=\mathcal{M}_{z}\{\langle Y_{m}(t;\xi,\eta)\rangle\}$. With these $\langle{Y_{m}(t)}\rangle$, the expectation value of any system operator can be conveniently calculated with the following relation:

[TABLE]

When the bath part of the operator $B(t)$ is not an identity operator, we just need to add an additional term in the noise average. Consequently, the expectation value can be calculated in a similar way,

[TABLE]

As we have mentioned, the effects of the bath are contained in the stochastic bath part of the operators instead of the artificial noise terms $z_{i,s_{i}}$, $i=1,2,3,$ and $4$. However, after the transformation in Eq. (LABEL:tron), the information of the bath is absorbed into the noise terms $\xi_{t}$, $\eta_{t}$, and $f(t,\zeta)$. The influence of the bath on the system dynamics is described by $\xi_{t}$ and $\eta_{t}$; concerned bath operators are provided by $f(t,\zeta)$; other details of the bath are eliminated as in the master equation.

III The dynamics of two-level systems

In this section, we will calculate the dynamics of different two-level systems as applications of our method. Two-level systems are the simplest kind of multi-level systems, and have a wide range of applications (see, e.g., mtt ; uosct5 ; uosct2 ). Meanwhile, the Langevin equation of a two-level system can also have the problem of nonlinear time-non-local terms. We will consider the coupling energy $\langle H_{\rm I}(t)\rangle$ and the bath displacement $\langle x(t)\rangle\equiv\langle\sum_{k}(g_{k}^{*}a_{k}^{{\dagger}}(t)+g_{k}a_{k}(t))\rangle$ as two simple cases of quantities which contain bath operators. Note that the bath displacement has units of frequency instead of length, because it is multiplied by the coupling coefficient. Such a displacement can better reveal the influence of the bath on the system. The stochastic expectation values of the system part can be obtained from Eq. (LABEL:ssle4), so we will directly come to the bath part. According to Eq. (1) and Eq. (LABEL:pebot2), we have

[TABLE]

Equation (LABEL:tlsbp) can be evaluated with Eq. (31) and Eq. (LABEL:seobu),

[TABLE]

From Eqs. (LABEL:tlsbpe), we can find that the contribution of the bath part can be represented by a new noise term $\zeta_{t}$. The correlation among $\zeta_{t}$, $\xi_{t}$, and $\eta_{t}$ is described through

[TABLE]

Then, the coupling energy can be calculated with Eq. (LABEL:saoo2) by introducing an additional noise to the noise average,

[TABLE]

where, $S_{\rm s}(t;\xi,\eta)$ is the system part stochastic expectation value of the coupling operator, $S(t)$, in Eq. (1). The bath part of the bath displacement is just the same as the one of the coupling energy. Therefore, we have the following relation:

[TABLE]

We will mainly consider two kinds of spectral densities for the baths. The first spectral density is the Ohmic form with a Debye regulation (see, e.g., ods ):

[TABLE]

where, $\Gamma$ is the coupling strength, and $\omega_{c}$ is the cut-off frequency. We will calculate several well-known simple models with this spectral density. In these models, we will express all the parameters in normalized units, as the ratio to the frequency of the two-level system $\omega_{0}$.

The second kind of spectral density we consider is the super-Ohmic spectrum with an exponential cut-off:

[TABLE]

where the coupling strength is described by $\alpha$. This spectral function is frequently used in solid-state quantum dot systems, e.g., when coupled to acoustic phonons. We will first reproduce some known results, with our alternative techniques, and then obtain some new results. In addition, the bath displacement, which cannot be obtained with former methods, will also be calculated. The parameters in this part will be expressed with units of picosecond, like other works modeling real quantum dot systems papi1 ; papi2 . In our numerical calculations, the noise functions are generated with the FFTW3 pack fft .

III.1 Pure dephasing model

The first model we consider is the so-called pure dephasing model:

[TABLE]

We use such an analytically solvable model to check our method. According to Eq. (LABEL:ssle4), the stochastic c-number Langevin equations of the pure dephasing model are

[TABLE]

The noise averages of $\langle\sigma_{xs}(t;\xi,\eta)\rangle$, $\langle\sigma_{ys}(t;\xi,\eta)\rangle$, $\langle\sigma_{zs}(t;\xi,\eta)\rangle$, and $\langle I_{\rm sys}(t;\xi,\eta)\rangle$ correspond to the expectation values of three Pauli matrices and the identity operator. We assume the initial state of the system to be the eigenstate of $\sigma_{x}$ $|\psi_{0}\rangle=\frac{1}{\sqrt{2}}(|e\rangle+|g\rangle)$, where $|g\rangle$ ($|e\rangle$) is the ground state (excited state). The initial state of the bath is assumed to be a thermal state. The coupling energy between the system and the bath is

[TABLE]

From the analytical solution of this model, we obtain the expectation values of the following quantities,

[TABLE]

With Eq. (52), the stochastic results can be compared with the analytical results. We first show the results at low temperature $\beta=1000$ ($\beta$ has units of $1/\omega_{0}$).

From Fig. 1, we can see that the converged part of the stochastic results (solid curves) agree well with the analytical results (dotted curves). Figure 1(a) shows the expectation values of the system operators; there, the stochastic results are nearly converged except for some tiny deviations at long times. However, the result in Fig. 1(b), which shows the coupling energy, is not that good. For comparatively long times, $\omega_{0}t>15$, the stochastic result oscillates severely around the analytical result. Also note that the interaction energy agrees well with the analytical even for short times, so the approximation in getting the bath operators is reliable. Then, we introduce the accumulated errors of $\langle\sigma_{x}\rangle$ and $\langle H_{\rm I}\rangle$ to characterize the convergence of our method. The accumulated error of some operator $B(t)$ is calculated from

[TABLE]

The $\langle B(s)\rangle_{\rm a}$ and $\langle B(s)\rangle_{\rm s}$ correspond to the analytical result and the stochastic result, respectively. As shown in Fig. 1(c), the error decreases significantly with the number of trajectories for the system operator $\sigma_{x}$. However, the calculation of the coupling energy is much more difficult. According to the results in Fig. 1(c), doubling the number of trajectories (black dotted curve and red dashed curve) makes only a few percent difference in the error. Only by increasing the number of the trajectories from $5\times 10^{6}$ (black dotted curve) to $45\times 10^{6}$ (green solid curve), we can reduce the error slightly. Therefore, the converging properties of the system operators and the bath ones are quite different.

Next, we come to the case of high temperature $\beta=1$. Compared with Fig. 1, the convergence of the results in Fig. 2 is much better. This is not surprising because the system is closer to the classical limit at high temperature. However, the convergence of the coupling energy $\langle H(t)\rangle$, shown in Fig. 2(b), is still not as good as the system operators, shown in Fig. 2(a).

III.2 Pure dephasing model with classical field control

We now consider a situation with classical field control. In addition to the original Hamiltonian in Eq. (49), we will add a classical control field Hamiltonian

[TABLE]

The coefficient $C(t)$ describes the shape of the control field, and we choose the ideal $\pi$-pulses as the control field. The separation between two pulses is $2/\omega_{0}$. Thus the coefficient $C(t)$ can be written as

[TABLE]

The derivation in Sec. II is based on time-independent Hamiltonians, but this formalism also works for time-dependent Hamiltonians (see Appendix D). The numerical results are shown in Fig. 3.

In Fig. 3, our results correctly show the effects of the control field. Compared with the results without control field in Fig. 2, the coherence time is much longer. At long times, the expectation values of the system operators become the steady values, but the coupling energy cannot become a steady value under the influence of the control pulses. The sudden changes in the Fig. 3 correspond to the control pulses. The $\langle\sigma_{z}(t)\rangle$ is not shown as it is a conserved quantity for our initial state.

III.3 Spin-Boson model

The second model we consider is the spin-boson model, which, unlike the pure dephasing mode, cannot be solved analytically. The spin-boson model is very important as it can describe light-matter interaction problems and double-well potential problems. The system Hamiltonian and coupling operator of the spin boson model are,

[TABLE]

The stochastic c-number Langevin equations of the system operators are as follows,

[TABLE]

The coupling energy can be estimated with

[TABLE]

and we assume the bath to be in the thermal state, while the system in the excited state $|e\rangle$.

The convergence of the results, as shown in Fig. 4, is comparatively good, and there is only some small fluctuations at the end of the red curve in Fig. 4(a). Similar to the pure dephasing model, the coupling energy at low temperature is hard to calculate. The coupling energy decreases monotonically when the temperature is low (Fig. 4(a)). In the high temperature case (Fig. 4(b)), the coupling energy increases first, then decreases to a steady value. As mentioned above, our calculation of the coupling energy can be unreliable for short times, so this difference should be treated carefully. The steady-state value of the interaction energy at higher temperatures is much lower than the low-temperature one. Then we add an optical pump term to the system,

[TABLE]

The peak Rabi frequency is set to $\Omega=0.5$ (in units of $\omega_{0}$), and we consider different detunings. We assume the initial state of the system its ground state. From Figs. 4(c) and 4(d), we can find that the influence of the detuning is insignificant. This means that the effects of the pump are suppressed by the bath. However, note that the pump here only differs from $S=\sigma_{x}$, which describes the coupling to the bath, by a coefficient. We will show in the next example that the situation can be different when $S=\sigma_{z}$.

III.4 Optically pumped quantum dot system

Next we consider an optically pumped quantum dot system, which has many applications in solid state quantum optics. The system considered here is a quantum dot with a pump and electron-phonon coupling, which has been studied in various experiments (e.g., qdots1 ; qdots2 ). The Hamiltonian under the rotating frame is

[TABLE]

Here, the detuning of the pump is $\delta$, and the Rabi frequency of the pump is $\Omega(t)$. The spectrum of the bath is of super-Ohmic form. We first consider the case of on-resonant pumping $\delta=0$. The parameters used here are from reference polatran3 . In addition to the population of the excited state, which is ${(1+\langle\sigma_{z}\rangle)}/{2}$, we also calculate the bath displacement induced by the system $\langle x(t)\rangle=\sum_{k}\langle(g_{k}^{*}a_{k}^{{\dagger}}(t)+g_{k}a_{k}(t))\rangle$.

From the red dashed curves in Figs. 5(a)-(c), our calculations agree well with other results polatran3 ; namely, the population oscillations for different pump Rabi frequencies, and there is bath-induced damping. The bath displacements (black solid curves) can provide more insights into the bath dynamics and coupling. When the pump intensities are moderate [Figs. 5(b) and 5(c)], the bath displacements are comparatively large and oscillate with the same frequencies of the oscillations of $\langle\sigma_{z}\rangle$. The bath displacement is still large if the pump is weak [Fig. 5(a)], but the oscillation is not resonant with the system. For very strong pumping [Fig. 5(d)], the bath displacement is very small, and there is little bath-induced damping; this is consistent with other findings polatran3 . The bath displacement here agrees with the estimation of the variational polaron transform master equation, namely, nearly zero bath response at large pump Rabi frequency. In addition, our method can also provide a clear picture of the bath response when the Rabi frequency of the pump is close to the cut-off frequency. In such a case, the bath displacement and the system population always oscillate at the same frequency, and the dissipation reaches its maximum.

Next we consider the time-dependent pump intensity polatran8 . The parameters of the bath are the same as the ones in the time-independent case. We set the detuning to be $\delta=-1.26~{}{\rm ps}^{-1}$, and the pump intensity to be $\Omega(t)=\Omega_{0}\exp(-(t/\tau)^{2})$ with $\Omega_{0}=1.28~{}{\rm ps}^{-1}$ and $\tau=20.2~{}{\rm ps}^{-1}$ corresponding to the pulse area $\Theta\equiv\sqrt{\pi}\Omega_{0}\tau=14.6\pi$. As shown in Figs. 5(e) and 5(f), the populations (red dashed curves) are inverted with negative detunings. Also, the steady values of the bath displacements are not zero in these cases. This is consistent with other works and experiments papi1 ; papi2 ; papi3 , which shows that the phonon dissipation can assist to invert the population of a pumped quantum dot when the detuning is negative.

In this section, we have calculated the expectation values of different operators in two-level systems. Although our results do not totally converge in Fig. 1(b), the coupling energy has reached its steady-state value before the results diverge. For other examples, the trajectory numbers are sufficient to obtain convergence. The numerical efficiency of our approach can change with the problem studied. Here, we summarize the numbers of trajectories used for different results in the following table:

We have shown that our method can obtain the expectation values of operators which contain bath operators. Also, our method can deal with both low and high temperature cases. However, the low temperature or the non-trivial bath parts can increase the calculation costs.

IV conclusions

We have developed a stochastic c-number Langevin equation method to access both the system information and the bath information. The problem of the nonlinear time-non-local terms in Langevin equations is avoided by formally dividing all the operators into system parts and bath parts with auxiliary stochastic fields (noise). As a Heisenberg-Langevin method, our approach can conveniently access the bath information. Such information about the bath can be quite different even when the dynamics of Pauli matrices [$\sigma_{x}(t)$, $\sigma_{y}(t)$, and $\sigma_{z}(t)$] are similar. In addition, this method is not limited to certain bath spectra or temperature in spite of the increased computing cost at low temperatures. We have also applied our method to several cases. For example, our equations work well in different well known simple models, like the pure dephasing model and the spin-boson model, and can compute both system and bath quantities. We have also reproduced some existing results in a pumped quantum dot system, and extended our model into a range of validity where typically these methods fail, In addition, the bath displacement, which is difficult to obtain with former approaches, can be calculated with our methodology. Finally, we stress that our method can be applied to problems like closed many-body systems or quantum thermodynamics.

Acknowledgements.

We are grateful to Prof. Wei-Min Zhang and Chris Gustin for their instructive comments. Z.Y.Z. is supported by the Japan Society for the Promotion of Science (JSPS) Foreign Postdoctoral Fellowship No. P17821. Y.A.Y. is supported by the National Natural Science Foundation of China (Grant No. 21373064). J.Q.Y. is partially supported by the National Key Research and Development Program of China (Grant No. 2016YFA0301200) and the National Natural Science Foundation of China (Grant No. 11774022 and No. U1801661). F.N. is supported in part by the: MURI Center for Dynamic Magneto-Optics via the Air Force Office of Scientific Research (AFOSR) (FA9550-14-1-0040), Army Research Office (ARO) (Grant No. W911NF-18-1-0358), Asian Office of Aerospace Research and Development (AOARD) (Grant No. FA2386-18-1-4045), Japan Science and Technology Agency (JST) (via the Q-LEAP program, and the CREST Grant No. JPMJCR1676), Japan Society for the Promotion of Science (JSPS) (JSPS-RFBR Grant No. 17-52-50023, and JSPS-FWO Grant No. VS. 059. 18N), the RIKEN-AIST Challenge Research Fund, and the John Templeton Foundation. S.H. thanks the Natural Sciences and Engineering Research Council of Canada for funding.

Appendix A Stochastic bath evolution operator

The solution of Eq. (27) can be obtained with the Magnus expansion,

[TABLE]

Then we will separate the annihilation operators $a_{k}$ and creation operators $a_{k}^{{\dagger}}$ in Eq. (61) with the Baker-Campbell-Hausdorff formula,

[TABLE]

Appendix B The generator of the bath operators

By taking the functional variation on the noise term in Eq. (27), we can obtain,

[TABLE]

For sufficiently long times, we can obtain the creation operator of the $k$th bath mode with Fourier transforms. Note that this approximation may introduce some error in the short-time limit. We have

[TABLE]

The annihilation operator $a_{k}$ can be easily obtained by taking the transpose conjugate of Eq. (LABEL:gc).

Appendix C The stochastic expectation value of bath identity operator $\langle I_{b}(t;z)\rangle$

Let us now calculate the average value of the stochastic bath identity operator:

[TABLE]

The operators in Eq. (LABEL:apcsuk) can be rearranged in normal order as follows,

[TABLE]

The $x^{*}_{1,s}$ is the complex conjugate of the noise term $x_{1,t}$ in Eq. (24): Then, we can calculate the trace in Eq. (LABEL:apcsuk).

[TABLE]

Now, we define the correlation function at temperature $T$ as

[TABLE]

Subsequently, the stochastic expectation value of the bath identity operator in Eq. (LABEL:apsevbu) can be expressed in a compact form,

[TABLE]

Appendix D System with classical field control

The derivations in Sec. II do not allow for a time-dependent Hamiltonian, because we have used the property that the Hamiltonian is unchanged during the time evolution. Now we provide the derivation of the time-dependent case. Consider a system with control field,

[TABLE]

Here, the $H_{\rm sys}$ is the system Hamiltonian without control, $S_{l,\rm c}$ is a system operator, and $C_{l}(t)$ is the control function. Note that the Hamiltonian in Eq. (69) is an approximate one. The original Hamiltonian should be

[TABLE]

The control field is also governed by the quantum dynamics. The contribution of the external Hamiltonian $H_{\rm ext}$ is usually assumed to be very large (classical limit). Therefore, the evolution of the operator of the external control $C_{l,\rm ext}$ is only decided approximately by $H_{\rm ext}$. If we further omit the entanglement between the system and the external field, the operator of the external field $C_{l,\rm ext}$ can be substituted with its classical expectation value, so that

[TABLE]

The effective Hamiltonian in Eq. (71) is obtained in the Schr$\rm\ddot{o}$dinger picture, so it can not be used directly in the Heisenberg picture. Now, we start from the original Hamiltonian in Eq. (70), which can be solved with our method:

[TABLE]

According to Eq. (LABEL:ssle1), the stochastic equations of the external field operator $O_{\rm ext}$ and system operators are,

[TABLE]

The influence of the system on the external field is not considered here. Then, the dynamics of the external field is just the free dynamics,

[TABLE]

with the initial state of the external field $|\psi_{\rm ext}\rangle$. Then we “trace over” the degrees of freedom of the external field in Eq. (LABEL:apssle1) by taking the expectation value:

[TABLE]

Here we have use the property $\langle\psi_{\rm ext}|\frac{\partial}{\partial t}Y_{l}(t;z)|\psi_{\rm ext}\rangle=\frac{\partial}{\partial t}Y_{l}(t;z)$ because the entanglement between the system and the external field is neglected. With Eq. (LABEL:apssle2), we can follow the way in Sec. II and obtain the time-dependent stochastic equations.

[TABLE]