Multi-State Quantum Dissipative Dynamics in Sub-Ohmic Environment : The Strong Coupling Regime

We study the dissipative quantum dynamics and the asymptotic behavior of a particle in a bistable potential interacting with a sub-Ohmic broadband environment. The reduced dynamics, in the intermediate to strong dissipation regime, is obtained beyond the two-level system approximation by using a real-time path integral approach. We find a crossover dynamic regime with damped intra-well oscillations and incoherent tunneling and a completely incoherent regime at strong damping. Moreover, a nonmonotonic behavior of the left/right well population difference is found as a function of the damping strength.


Introduction
Quantum systems are unavoidably subject to interactions with a large number of environmental degrees of freedom.Due to such interactions, an open system evolution needs to be described through the joint dynamics of system and environmental degrees of freedom, which together can be modeled as a unique isolated system.The reduced dynamics of the subsystem is then obtained by tracing out the environmental degrees of freedom.As a result, the system experiences loss of coherence, dissipation and relaxation towards equilibrium.
An effective microscopical description of the dissipation is given by the celebrated Caldeira-Leggett model [1], in which a quantum particle is subject to the linear interaction with a reservoir of N independent quantum harmonic oscillators.In the thermodynamical limit N → ∞ the reservoir is called a heat bath.The so-called spectral density function J(ω), describing the dependence of the coupling on the frequency, is assumed to be of the form J ∝ ω s e −ω/ωc , where ω c is a high-frequency cut-off.The special case s = 1 describes the so-called Ohmic dissipation.In most cases the Ohmic dissipation gives a good description of the effects exerted by the thermal bath, such as in Josephson flux qubits [2].
However, sub-Ohmic environments (s < 1) are of interest both from experimental and theoretical point of view.For example in single electron tunneling the electromagnetic environment constituted by an RC transmission line is sub-Ohmic with s = 0.5 [2].A sub-Ohmic environment with s = 0.5 can be also realized by capacitively coupling an RLC transmission line to a mesoscopic metal ring [3].Further examples of sub-Ohmic environments are found in nanomechanical devices [4].The sub-Ohmic case, for an unbiased bistable system at T = 0, is characterized by a phase transition from the delocalized to localized phase [5,6].This means that, in right/left well state representation (|R , |L ), one has σ z = 0 at t → ∞.The sub-Ohmic regime is interesting also because of the larger relative weight of the low frequency modes with respect to the Ohmic case.This feature reproduces the 1/f -noise in the s → 0 limit, which constitutes an important decoherence source in superconducting qubits, due to the presence of impurities [7].
A particle confined in a double well potential is a paradigmatic example of a bistable system, both in classical and in quantum physics.The confining potential is characterized by two minima separated by a barrier of height ∆U .The most important information of interest in such a system is the population distribution dynamics across the two minima.Classically, the only mechanism which allows for passage from one well to the other is the thermally activated jump over the barrier.On the other hand, the peculiar feature of the quantum regime is the tunneling, which allows probability amplitudes to diffuse across the barrier.Recently, the dynamics of quantum bistable systems has gained renewed interest [8][9][10][11][12].
The strong nonlinearity of the potential causes the energy levels below the potential barrier to group into well-separated doublets.If the temperature is much smaller than the inter-doublet separation (∆ inter , see Figure 1) between the first two doublets, then the system can be considered as an effective two-level system (TLS).The resulting dissipative TLS model is the so-called spin-boson (SB) model [13], as the excitations of the reservoir obey the Bose-Einstein statistics.
In the intermediate to high temperature regime, with respect to ∆ inter , and for initial preparations involving higher energies, the TLS approximation breaks down and the next-lying energy states are involved.Contrary to the SB model, which has been studied in every dissipation regime, both Ohmic and non-Ohmic [2,[13][14][15][16], dissipative multi-state systems, beyond the weak coupling regime, have been less investigated, especially for non-Ohmic environments.
The reason is that the traditional Born-Markov master equation techniques, specifically the Bloch-Redfield master equation [17][18][19], which are perturbative in the system-bath coupling, describe well the weak coupling regime.However, whenever the coupling cannot be treated as a perturbation these techniques fail.The path integral approach, being non-perturbative in the coupling with the environment, is well suited in the intermediate to strong coupling regimes.An effective approach is based on numerical ab initio evaluations [20].However, when the Hilbert space dimension of the system goes beyond two and the memory time of the reservoir becomes too large, this technique requires large computational resources.Indeed, the memory time of such non-markovian dynamics is typically very long for weak/intermediate coupling and low temperatures.Here we use an alternative path integral approach based on an integro-differential equation (see Ref. [2]) which is of the Nakajima-Zwanzig type [19] with approximate kernels.
In the TLS case the dicotomic nature of the paths in the left/right state representation makes the sum over paths possible for the free system and allows for approximate treatments of the environmental influences.On the other hand, in the general case of a M -state system, the variety of paths to be considered makes the evaluation of the kernel matrix a hard task independently of the dissipation regime considered.Therefore, further approximations on the contributing paths are needed.Following Ref. [21] we consider a temperature/coupling regime where we can neglect a whole class of paths of the reduced density matrix (RDM), namely the paths with long off-diagonal excursions [22] (clusters), due to the cut-off operated by the thermal bath.

The Model
The system is a particle of mass M , coordinate q, and momentum p, subject to a double well potential V 0 .The particle is linearly coupled to an environment of N independent quantum harmonic oscillators of frequencies ω j , which, in the thermodynamical limit N → ∞, is the so-called bosonic heat bath.
The full Hamiltonian is the sum of a free system term, a free reservoir term and a system-reservoir interaction term The Hamiltonian has a renormalization term ∝ q2 giving a spatially homogeneous dissipation, not depending on the particle's coordinate q.The result is a so-called purely dissipative bath.
The interaction of the particle with the individual degrees of freedom of the reservoir is defined by the set of constants c j and is proportional to the inverse of the reservoir's volume [2].Thus, for a macroscopic environment, the coupling with the individual oscillators is weak, which justifies the linear coupling assumption.Nevertheless, the overall influence exerted by the heat bath as a whole can be strong.
The bath spectral density function, which describes the frequency distribution of the reservoir's oscillators and their coupling with the particle, is defined by and has the dimension of mass times frequency squared.
In the general case of continuous bath the spectral density function is modeled as a power of ω, characterized by the exponent s, with an exponential cutoff at ω c The bath is said sub-Ohmic for 0 < s < 1, Ohmic for s = 1 and super-Ohmic for s > 1.The so-called damping constant γ is a measure, in the continuous limit, of the system-bath coupling.
The phonon frequency ω ph is introduced in such a way that γ has the dimension of a frequency also in the non-Ohmic case (s = 1).The exponential cut-off at high-frequency is introduced to avoid non-physical results as, for example, the divergence of the renormalization term in the Hamiltonian of Equation ( 1).The effect of the high frequency modes is taken into account by a redefinition of the particle's bare mass, which is dressed by the high-frequency bath modes [2].
Within the Caldeira-Leggett model (see Equation ( 1)), the Heisenberg equation for the particle's position operator results in the following generalized quantum Langevin equation [2,23] where and are respectively the memory-friction kernel and the bath force operator.The expectation value of the bath force operator and the bath force autocorrelation, taken with respect to the environment in the canonical respectively.
In the classical limit (h → 0) the bath force correlation function is where we used coth(βhω j /2) ∼ 2(βhω j ) −1 for h → 0. Therefore the two relations in Equation (7), in the continuum limit (N → ∞), describe a stochastic force which in turn reproduce, in the classical limit, a classical colored noise source.In other words, non-Ohmic environments correspond to classical colored noise sources [24].By comparing Equation ( 5) with Equation (2) and taking the continuous limit one finds that the Ohmic case gives, in the classical limit, a white noise source.Note that in the right hand side of the quantum Langevin equation ( 4) a term is dependent on the initial condition q(0).This term is ascribable to the factorized initial condition with the reservoir in the canonical equilibrium and vanishes at long time due to the interference of the quasi-continuum of spectral components of γ(t).

Path Integral Representation of the Reduced Dynamics
We denote by W the full density matrix, characterized by the unitary evolution where the time evolution operator is induced by the Hamiltonian (1).The RDM, in the q representation, reads We assume the following factorized initial condition With this initial condition, the dynamics of the RDM is given by the exact formal expression ρ qq (t) = dq 0 dq 0 G(q, q , t; q 0 , q 0 , t 0 )ρ q 0 q 0 (t 0 ), (10) where the propagator G is a double path integral in the left (right) coordinate q (q ) G(q, q , t; q 0 , q 0 , t 0 ) = The transition amplitudes for paths with fixed extrema (q 0 , q 0 ) → (q, q ), in this sum-over-paths, are factorized as a bare system amplitude times the Feynman-Vernon influence functional F F V .This functional accounts for the influence of the environment [25] and tends to 1 for vanishing γ.

Feynman-Vernon Influence Functional in the Discrete Variable Representaion
The exact path integral expression (11) for the propagator can be computed only in special cases, namely the free particle and the dissipative harmonic oscillator [26].In the presence of nonlinear potentials, approximate treatments exist based on spatial discretization, attained by restricting the Hilbert space of the open system [27].Specifically, based on the assumption that, given the initial condition and the damping/temperature regime, the system is not likely to visit high-energy states during its time evolution, the Hilbert space of the system is restricted to that spanned by the first M energy eigenstates.Thus the problem reduces to that of a dissipative M -state system.In the present work we consider the so called double-doublet system, in which the first M = 4 levels of the quartic potential [22] V are considered.The discrete variable representation (DVR) is then obtained by diagonalizing the position operator in this restricted Hilbert space.The continuum of position states turns into a discrete set of states localized around a grid of 4 position eigenvalues q 1 , . . ., q 4 , where In this picture, the paths of the coordinates q and q are represented by a sequence of transitions among the positions q j .In Figure 1 the potential considered in this work, the first four energy levels, and the position eigenvalues are shown.
Figure 1.Potential V 0 (see Equation ( 12)) with ∆U = 1.4 hω 0 and = 0.02 M hω 3 0 .The intra-doublet separations are ∆ 1 0.12 hω 0 and ∆ 2 0.15 hω 0 , while the inter-doublet spacing is From now on we consider the time evolution of the populations ρ ii = q i |ρ|q i , i.e., the diagonal elements of the RDM in the DVR.Taking the influence functional as F F V = exp(−Φ F V ), the influence phase in the DVR for a path with N transitions at times t 1 , . . ., t N has the following form where the stepwise functions of time ξ i and χ i are the so called charges ξ j = (q j − q j−1 ) − (q j − q j−1 ) and the indexes refer to the transition times.
The time-integrated bath correlation function , also called pair interaction because it couples the charges at different times, is related to the bath autocorrelation function by where the dimensionless quantity κ is defined by κ = (hβω c ) −1 .The function Γ(z) is the Euler gamma function and ζ(z, q) is the Hurwitz zeta function ζ(z, q) = ∞ n=0 1/(q + n) z , related to the Riemann zeta function ζ(z) by ζ(z) ≡ ζ(z, 1).

Generalized Non-interacting Blip Approximation
According to the path integral jargon, a blip is an off-diagonal configuration of an RDM path for two-level systems, while, for M -state systems an off-diagonal excursion of an RDM path is called cluster.From the expression for the phase Φ F V (Equation ( 14)) the influence functional F F V = exp(−Φ F V ) is time non-local as it couples the ξand χ-charges at every transition time.This feature prevents an exact evaluation of the path integral expression for the populations.Nevertheless a multi-state generalization of the non-interacting blip approximation (NIBA) exists, known as generalized non-interacting cluster approximation (gNICA) [21].For Ohmic friction, in the high temperature limit, the function Q(t) assumes a linearized form at short times.This implies an almost exact decoupling of the clusters.Moreover if the damping γ is sufficiently strong one can retain, in the sum over paths, the leading contributions given by RDM paths with clusters consisting of single off-diagonal transitions or blips.The gNICA scheme to the leading order can thus be called gNIBA (generalized NIBA).
In the sub-Ohmic case the decoupling of the blips is attained at shorter times, and long clusters are suppressed more effectively with respect to the Ohmic case, due to stronger cut-off operated by Q in the real part of Φ F V .These two features make the gNIBA more accurate in the sub-Ohmic than in the Ohmic regime.This can be seen from Figure 2, where the real part of Q is shown as a function of time for s = 1 and 0.5, i.e., in the Ohmic and sub-Ohmic regimes, respectively.The cut-off localizes in time the influence of the environment, which results in the factorization of the path integral expression for the propagator in Laplace space.This allows to obtain a generalized master equation for the populations in the DVR [28] which, for the double-doublet system (M = 4), reads The non-diagonal elements of the gNIBA approximated kernel matrix K are where q 2 ij = (q i − q j ) 2 , with ĤS the bare system Hamiltonian.
The diagonal elements of K are given by the probability conservation

Results
In this section, we show the dynamics of the populations in the DVR, the spatially localized representation introduced in Section 3.1.The results are obtained by numerically integrating the generalized master equation given in Equation (17).
The quantum particle is subject to the biased bistable potential in Figure 1.The barrier height (at zero bias) and the bias are ∆U = 1.4 hω 0 and = 0.02 M hω 3 0 , respectively.The frequency ω 0 ∼ ∆ inter /h is the natural oscillation frequency around the potential minima (see Equation ( 12) and Figure 1).
Due to the asymmetry of the potential, at equilibrium we expect that P ≡ P R − P L > 0, where the left and right well population are defined by P L = ρ 11 + ρ 22 and P R = ρ 33 + ρ 44 .The value of the population difference P depends on the temperature and damping strength γ.
We consider the sub-Ohmic regime with s = 0.5 .The cut-off frequency is set at ω c = 50 ω 0 and the phonon frequency at ω ph = ω 0 (see Equation ( 3)).The generalized master equation given in Equation ( 17) is numerically solved for two temperatures, namely T = 0.2, 0.3 hω/k B .In all graphs shown in this section the system is initially in the out-of-equilibrium state ρ(0) = |q 1 q 1 |, that is the particle is initially prepared as a wave packet centered at the position q 1 (see Figure 1).This initial condition involves the states belonging to the higher energy doublet, implying that the TLS approximation is not valid whatever the temperature regime is.
The energy levels structure, shown in Figure 1, implies that the bare system has multiple time scales: a fast one, of the order of the inverse of the inter-doublet separation h/∆ inter , and a slower one, of the order of the intra-doublet separation h/∆ 1 ∼ h/∆ 2 .Notice that, due to the bias, the intra-doublet separation are of the same order, while in the unbiased (symmetric potential) case we have ∆ 1 << ∆ 2 [29].
In this section, both in the text and figure, temperatures and damping strengths are expressed in units of hω 0 /k B and ω 0 , respectively.

Dynamics of the Populations
Because the coupling with the environment is in the intermediate to strong regime, the slow tunneling dynamics is incoherent.This means that the oscillations between the two wells, predicted for the free system, are completely damped out.The reason is twofold: (i) transitions among populations of distant DVR states |q i and |q j are suppressed more effectively, due to the prefactor q 2 ij in the exponential cutoff of the gNIBA kernels (Equation ( 18)); (ii) the energy scale ∆ 1 ∼ ∆ 2 of the tunneling dynamics is smaller than that of the intra-well dynamics (given by the inter-doublet separation ∆ inter ).The effective damping is therefore strong on the scale of tunneling separation and intermediate on the scale of the inter-doublet separation.As a consequence, while the intra-well dynamics displays a damped oscillatory behavior, the tunneling dynamics is overdamped.This crossover dynamical regime, found also for Ohmic damping in Ref. [29], is shown in Figure 3 where two values of the temperature, T = 0.2 and 0.3, are considered with γ fixed to the value 0.1.Both the evolutions display a transient with coherent intra-well oscillations and an incoherent inter-well (tunneling) dynamics.At higher temperature (right panel) the left well population (ρ 11 + ρ 22 ) is larger than that at the lower temperature (left panel).Temperatures and γ are given in units of hω 0 /k B and ω 0 , respectively.Time is in units of ω −1 0 .
Figure 4 shows the time evolution of the populations at the stronger damping γ = 0.3.The dynamics occurs at the transition between the crossover and the completely incoherent regimes.This transition from coherent to incoherent dynamics has been investigated for the sub-Ohmic spin-boson in Ref. [30].At both temperatures investigated, the prediction is that no full intra-well oscillation in the transient takes place.Moreover the relaxation times, at both temperatures, are larger for this value of γ with respect to that considered in Figure 3.
Note that, due to the difference in time scales between the intraand inter-well dynamics, the intra-well relaxation of the left well populations occurs before the tunneling takes place and a sort of metastability appears during the transient before the relaxation to the equilibrium takes place.A similar behavior, but of different nature, has been found in Ref. [31].There, for a symmetric two-level system in the sub-Ohmic regime, a slowly decaying quasiequilibrium state is present during the transient, due to a shifted initial preparation.This kind of shifted equilibrium initial condition is of experimental interest and is worth investigating for multi-state systems in sub-Ohmic baths in future developments.The coherent intra-well oscillations are damped and the relaxation dynamics is in the transition from the crossover to the completely incoherent relaxation regime.As well as for γ = 0.1 (Figure 3), at higher temperature (right panel) the left well population (ρ 11 + ρ 22 ) is larger than that at the lower temperature (left panel).Temperatures and γ are given in units of hω 0 /k B and ω 0 , respectively.Time is in units of ω −1 0 .
Finally, both in Figures 3 and 4, the internal states |q 2 and |q 3 , with wave packets centered away from the potential minima (see Figure 1), are more populated at the higher temperature.This effect is evident at the lower value of damping (Figure 3).The reason is that: (i) at equilibrium, the higher doublet of energy states is relatively more populated at T = 0.3 than at T = 0.2; (ii) the energy states of the higher doublet contribute more to the internal DVR states than to the external ones (|q 1 and |q 4 , centered near the minima).As a consequence ρ eq 22/33 (T = 0.3) > ρ eq 22/33 (T = 0.2).

Population Difference-Time Evolution and Equilibrium Configuration
In this section, we study the time evolution and stationary configuration of the population difference P = P R − P L .In Figure 5 we plot P as a function of time, for the two temperatures T = 0.2, 0.3 and different values of γ.
The intra-well oscillations of ρ 11 and ρ 22 , shown in Figures 3 and 4, are not visible in Figure 5.This is because we consider the overall left/right well populations.The time evolution of P is a monotonic relaxation towards an equilibrium configuration which depends on the temperature and damping strength.
All plots have qualitatively the same features but are shown to highlight how the relaxation time and the equilibrium value of P behave, at the considered temperatures, with respect to an increase of the coupling with the environment, which is quantified by γ.
Specifically, at T = 0.2, 0.3, the relaxation time grows exponentially with γ, at least around γ = 0.1 (notice that the plot is in log-scale).This behavior is similar to that predicted for Ohmic damping [21,32].
As expected, for fixed γ, P (t → ∞) is higher at the lower temperature T = 0.2, which indicates that the probability to find the particle in the left well at equilibrium is lower at this value of temperature with respect to the higher value T = 0.3.Moreover, P eq ≡ P (t → ∞) behaves nontrivially with respect to the coupling, especially at T = 0.2 (Figure 5, left panel), as shown in Figure 6 where P eq is plotted as a function of γ.The curve displays a maximum around γ = 0.15.This prediction may be ascribed to the multi-level nature of the system studied.Specifically, two competing effects are present.On the one hand, an increase of the coupling forces the relaxation towards the lower well, producing an increase of P eq .On the other hand, the higher energy doublet gets more populated and the effective damping is reduced, with a corresponding decrease of P eq .The combined effect of these two mechanisms is a possible explanation of the nonmonotonicity found.This result, together with the multiple dynamical time scales and the appearance of the so-called vibrational relaxation (intra-well motion), is another peculiar feature of the picture arising from treating the system beyond the TLS approximation. Figure 6.Population difference at equilibrium P eq = P (t → ∞) as a function of γ (in units of ω 0 ) at temperature T = 0.2 hω 0 /k B (see left panel of Figure 5).

Conclusions
In this work, we investigate the relaxation dynamics and the equilibrium configuration of a biased bistable system in a strongly dissipative sub-Ohmic environment.We consider an out of equilibrium initial condition, with the particle localized in the higher well of the bistable potential.The time evolution of the populations in a space localized representation is calculated by using a real-time path integral approach.A NIBA-like approximation scheme has been used which is suitable for the dissipation regimes considered, especially for a sub-Ohmic spectral density.
The system is studied beyond the two-level system approximation.Specifically the Hilbert space is restricted to that spanned by the states belonging to the first two energy doublets.This feature, along with the specific initial condition chosen, results in the presence of multiple time scales in the relaxation dynamics.We observe a crossover dynamical regime, characterized by damped intra-well oscillations and incoherent tunneling, at intermediate damping.At strong damping the dynamics makes a transition to the completely incoherent relaxation.
The study, performed for two different values of temperature, shows a dependence of the stationary configuration on the temperature and on the damping strength.Specifically, at fixed damping, the higher well is more populated at the higher temperature while, at fixed temperature, the left well has a minimum of the population when the damping value roughly coincides with the higher intra-doublet frequency spacing.This result, along with the observation of intraand inter-well dynamics, is a peculiar feature of treating the system beyond the two-level system approximation.

Figure 2 .
Figure 2. (Color online) Real part of the function Q (Equation (16)) vs time for s = 0.5 (sub-Ohmic) and s = 1 (Ohmic).The damping strength is γ = 0.1 ω 0 and the temperature T = 0.1 hω 0 /k B .The inset shows the behavior of Q at short times.

Figure 3 .
Figure 3.Time evolution of the four populations ρ ii in the intermediate damping regime γ = 0.1 ω 0 for sub-Ohmic damping (s = 0.5).Left panel-T = 0.2.Right panel-T = 0.3.Both the evolutions display a transient with coherent intra-well oscillations and an incoherent inter-well (tunneling) dynamics.At higher temperature (right panel) the left well population (ρ 11 + ρ 22 ) is larger than that at the lower temperature (left panel).Temperatures and γ are given in units of hω 0 /k B and ω 0 , respectively.Time is in units of ω −1 0 .

Figure 4 .
Figure 4. Time evolution of the four populations ρ ii in the strong damping regime γ = 0.3 for sub-Ohmic damping (s = 0.5).Left panel-T = 0.2.Right panel-T = 0.3.The coherent intra-well oscillations are damped and the relaxation dynamics is in the transition from the crossover to the completely incoherent relaxation regime.As well as for γ = 0.1 (Figure3), at higher temperature (right panel) the left well population (ρ 11 + ρ 22 ) is larger than that at the lower temperature (left panel).Temperatures and γ are given in units of hω 0 /k B and ω 0 , respectively.Time is in units of ω −1 0 .

Figure 5 .
Figure 5.Time evolution of the population difference P = P R − P L in the sub-Ohmic dissipation regime (s = 0.5), for different values of γ, at T = 0.2 (left panel) and T = 0.3 (right panel).Temperatures and γ are given in units of hω 0 /k B and ω 0 , respectively.Time is in units of ω −1 0 .