text:
# Semiclassical wave packet tunneling in real-time

## I Introduction

Tunneling through a potential barrier is one of the most fascinating aspects of quantum mechanics. In recent years a particular challenge has been to understand tunneling in complex systems using semiclassical methods. However, any simple description is seemingly hampered by the fact that a quantum mechanical object running towards a barrier with a typical energy smaller than the barrier height may penetrate it even though all classical trajectories with such energies are reflected completely. Thus, for static observables like e.g. tunnel splittings one works in the energy domain and calculates the energy dependent Green’s function semiclassically by switching from a real-time orbit outside the barrier to an orbit in imaginary time, i.e. with imaginary momentum, under the barrier. This technique can be traced back to the “old” WKB approximation landau  and meanwhile has been successfully extended to extract e.g. tunnel splittings also for systems with classically chaotic dynamics creag . The crucial question is then: Can semiclassical tunneling also be described in the real-time domain? This issue has turned into a fundamental challenge for our understanding of semiclassics in general and systems with explicit time dependence as e.g. in the context of driven tunneling and control of tunneling grifoni  or tunneling in the presence of chaos chaos  in particular.

The probability amplitude for a particle initially at $`q_i`$ to be at $`q_f`$ after time $`t`$ is given by Feynman’s path integral representation of the propagator as schulman

$$G(q_f,q_i,t)q_f|\mathrm{e}^{iHt/\mathrm{}}|q_i=𝒟[q]\mathrm{e}^{iS[q]/\mathrm{}}.$$
(1)
The integral sums over all paths running from $`q_i`$ to $`q_f`$ in time $`t`$ where each contribution is weighted according to its action $`S[q]=_{q_i}^{q_f}𝑑q\sqrt{2m[EV(q)]}Et`$ with potential $`V(q)`$ and energy $`E=E(q_f,q_i,t)`$. Due to the oscillating integrand in (1) tunneling appears to be the result of a complex interference pattern. In the semiclassical limit the path sum is dominated by the contributions of the stationary paths $`\delta S[q_{\mathrm{cl}}]=0`$ obeying Newton’s equation of motion and small fluctuations around them. In the last decade efficient semiclassical propagation schemes based on Gaussian wave packets have been developed, certainly the most powerful known as the Hermann-Kluk propagator (HK) kluk . However, the inclusion of deep tunneling (may be even in presence of external driving) has not been satisfactory yet. In fact, it was found that classically allowed real-time trajectories running over the barrier are not sufficient to capture strong tunneling heller2 ; heller1 . Various extensions have thus been attempted. By propagating a large number of initial wave packets tunnel splittings in a double well potential could be extracted, however, in single wave packet motion tunneling effects were absent mandel . As for barrier penetration no low energy stationary orbit exists along the real axis, it is tempting to think that one may find one in the complex plane heller1 . An individual path circumventing–in real-time–the barrier region in a complex coordinate plane is assumed, but its existence is still elusive. In grossmann  orbits run along the real axis, but the semiclassical propagation over the time interval $`t`$ is time sliced into steps over intermediate intervals. This spawning of orbits turns out to be numerically extremely expensive already for two slices and improves tunneling amplitudes only for energies not too far below the barrier top.

In the sequel we re-examine the semiclassical barrier penetration through scattering potentials starting from the exact expression (1). Our idea is this: While in the energy domain tunneling is described within a complex time plane, here, we study classical mechanics for complex energies. Our analysis reveals how tunneling is encoded in the quantum propagator (1) in terms of real-time orbits. In particular, it turns out that an individual complex “tunneling path” does not exist. Based on these results we extend the conventional HK to semiclassical wave packet dynamics in the deep tunneling regime not only through static potentials, but, for the first time, report also on accurate results for driven tunneling. The approach is shown to be very efficient for one-dimensional systems and may thus serve as a promising starting point for higher dimensional studies.

## II Complex mechanics and semiclassical approximation

We consider the motion of a particle of mass $`m=1`$ in a general one-dimensional, symmetric barrier potential $`V(q)`$ where the barrier top is located at $`q=0`$ with $`V(0)=V_0`$. $`V(q)`$ is assumed to be a smooth and analytic function of $`q`$ that can be approximated around $`q=0`$ by an inverted harmonic oscillator and for large $`|q|`$ falls off as $`V(q)V_0/[q/l]^{2k},k2`$, integer, with a typical barrier length scale $`l`$. A sufficiently high barrier is taken for granted. Typical examples include $`V_k(q)=V_0/[1+(q/l)^2]^k`$, but our results also apply to the Eckart barrier $`V_0/\mathrm{cosh}(q/l)^2`$ and the Gaussian barrier $`V_0\mathrm{exp}(q^2/l^2)`$.

Now, think of a wave packet $`\psi (q_i,0)`$ localized to the far right ($`q_i>0`$) which is propagated towards the barrier according to $`\psi (q_f,t)=𝑑q_iG(q_f,q_i,t)\psi (q_i,0)`$. We are interested in that portion of the packet that after time $`t`$ arrives on the far left ($`q_f=q_i<0`$). All real orbits connecting the two asymptotic regions run over the barrier ($`E>V_0`$) and as, for fixed end-points, $`t`$ becomes large they spend most of their time in the parabolic range around $`q=0`$. There, the marginal stability of trajectories causes the semiclassical $`G(q_i,q_i,t)`$ to die out exponentially in contrast to exact results heller2 . Classical orbits with $`E<V_0`$ coming from the far right or left reach the right or left flank of the barrier at turning points (TPs) $`q_0`$ and $`q_0`$, respectively; the long time properties of the path integral (1) are therefore governed by the dynamics in the “forbidden” range between the TPs. Mathematically, for $`E<V_0`$ no real stationary phase point to (1) obeying the proper boundary conditions exists in function space. The usual procedure is then an analytic continuation meaning here to extend classical mechanics to the complex coordinate plane.

Newton’s equation of motion, $`\ddot{q}+V^{}(q)=0`$, where $`\dot{q}=dq/dt`$ and $`V^{}=dV/dq`$, translates for complex $`q=x+iy`$ into

$$\ddot{x}+r_x=0,\ddot{y}+j_x=0.$$
(2)
Here, $`V(q)=r(x,y)+ij(x,y)`$ and the subscript $`x`$ \[$`y`$\] denotes the partial derivative with respect to $`x`$ \[$`y`$\]. We further exploited that for analytic functions $`V(q)`$ Cauchy’s relations $`r_x=j_y`$ and $`r_y=j_x`$ apply. From (2) one simply finds that the total energy $`E=ϵ_{\mathrm{re}}+iϵ_{\mathrm{im}}`$ and its real and imaginary parts

$$ϵ_{\mathrm{re}}=(\dot{x}^2\dot{y}^2)/2+r(x,y)\text{ and}ϵ_{\mathrm{im}}=\dot{x}\dot{y}+j(x,y),$$
(3)
respectively, are constants of motion. How does the corresponding classical mechanics look like? For low energies the TPs $`q_0,q_0`$ lie in the range where $`V(q)`$ can be approximated by its asymptotic behavior. Hence, we consider paths starting from large $`q_i=x_i>0`$ along the real axis with complex momentum $`p_i\dot{q}_i=\dot{x}_i+i\dot{y}_i`$. Typical trajectories are depicted in fig. 1. While basically three kinds of orbits can be distinguished, the common behavior is that as the barrier vanishes asymptotically, for large distances from the top the classical motion tends to be a free motion. If we represent trajectories in the form $`q(t)=R(t)\mathrm{e}^{i\varphi (t)}`$, for very large $`R`$ they run close to straight lines with constant $`\varphi (q_i,p_i)`$ depending merely on the initial phase space variables. Let us now discuss the types of orbits in detail.

First, we look at class (a) as it is the only one where orbits cross the line $`x=0`$ ($`\varphi =\pi /2`$) to reach the other side of the barrier. For such orbits we need initially $`\dot{x}_i<0`$ (otherwise they would run away from the barrier anyway) and for the discussion assume $`\dot{y}_i>0`$ implying $`ϵ_{\mathrm{im}}=\dot{x}_i\dot{y}_i<0`$, see Eq. (3). We further note the asymptotic form of $`V(q)=r(q)+ij(q)`$ using polar coordinates:

$$r(R,\varphi )=\frac{\mathrm{cos}(2k\varphi )}{(R/l)^{2k}},j(R,\varphi )=\frac{\mathrm{sin}(2k\varphi )}{(R/l)^{2k}}.$$
(4)
One sees immediately that a successful crossing of the dividing surface must happen with $`\dot{y}>0`$ (and $`\dot{x}<0`$ of course). Namely, at $`\varphi =\pi /2`$ the imaginary part $`j`$ vanishes so that $`ϵ_{\mathrm{im}}=\dot{x}\dot{y}=\dot{x}_i\dot{y}_i<0`$ meaning $`\dot{y}>0`$. For the required energy one derives $`|ϵ_{\mathrm{im}}|>r(x_i,0)=V(x_i)`$. Starting, however, from the imaginary axis with $`\dot{x}<0`$ and $`\dot{y}>0`$, i.e. in the direction of decreasing $`|V(q)|`$, always generates an orbit reaching the asymptotic left side of the barrier far from the real axis ($`\varphi <\pi `$). We conclude that a simple “tunneling path” connecting the asymptotic segments of the real axis on either side of the barrier via a tour through the complex plane does not exist. This is in sharp contrast to tunneling for fixed energy. There, the energy dependent Greens function $`K(q_f,q_i,E)`$ exhibits stationary phase points in imaginary time corresponding to classical paths running with energy $`E`$ in the inverted potential through the barrier range from $`q_i`$ to $`q_f`$. Here, for real-time tunneling a stationary phase path to the quantum propagator $`G(q_f,q_i,t)`$ even with complex energy cannot be found. This important result may also reflect the quite different roles “energy” and “time” play in quantum mechanics.

For our analysis the consequences are two-fold: on the one hand complex trajectories in class (a) do not play any role for a semiclassical approximation to $`G(q_f,q_i,t)`$, and on the other hand the path integral in (1) is in the low energy sector completely determined by fluctuations. To find its dominant contributions thus means to detect the dominant fluctuations; these are points in function space which lie close to orbits with $`\delta S[q]=0`$ and also obey the proper boundary conditions. Accordingly, we consider the remaining two classes of paths.

The second class (b) contains orbits with small but non-vanishing energies $`0<|ϵ_{\mathrm{im}}|<r(x_i,0)`$ which may exhibit TPs in the complex plane and always live on the same side of the barrier. Hence they are not relevant either. In the third class (c) trajectories have real total energy $`E`$, i.e. $`ϵ_{\mathrm{im}}=0`$, but start with purely imaginary momenta $`\dot{x}_i=0`$ and, as assumed, small $`ϵ_{\mathrm{re}}=\dot{y}_i^2/2+r(x_i,0)`$. These orbits display crucial features as we will explain in the following. For that purpose we focus on the limit $`ϵ_{\mathrm{re}}=0`$ and follow paths with $`x_i>0,\dot{y}_i>0`$. Writing asymptotically $`q(t)=R(t)\mathrm{e}^{i\varphi }`$ one obtains

$$\varphi (x_i,\dot{y}_i)\varphi _c^+=\pi /[2(k+1)].$$
(5)
Hence, after a transient period of time all those orbits run along the same line in the complex plane independent of their starting points $`x_i`$. And since they carry the same energy, they are also focused in phase space so that the line $`\varphi _c`$ defines a caustic. Were the trajectories optical rays, $`\varphi =\varphi _c^+`$ would be a burning line. Due to symmetry the same holds true for the complementary line $`\varphi _c^+`$ and the lines $`\pi \pm \varphi _c^+`$ on the other side of the barrier. Typically, a caustic is associated with unstable orbits and fluctuations connecting them which renders a simple Gaussian semiclassics insufficient schulman . To verify this scenario here, we consider small deviations $`\delta q=\delta x+i\delta y`$ around a certain orbit $`\overline{q}(x_i;t)`$. By linearizing the equations of motion (2) one gains $`(\delta \ddot{x},\delta \ddot{y})^T=`$$`𝑴`$ $`(\delta x,\delta y)^T`$ where $`𝑴`$ is the stability matrix evaluated along $`\overline{q}(t)`$. Along $`\varphi =\varphi _c^+`$ its diagonal elements are $`\overline{r}_{xx}>0`$ and the off-diagonal elements vanish $`\overline{j}_{xx}=0`$. Accordingly, all trajectories merging along the burning lines are unstable. Small deviations in phase space can lead from an orbit $`\overline{q}(x_i;t)`$ to another one $`\overline{q}(x_i^{};t)`$ and even allow for a turn from positive to negative momentum to run along the $`\overline{q}(x_i^{};t)`$-orbit back towards the real axis. As asymptotically paths creep along $`\varphi _c^+`$, jumps from very small positive to negative momenta require only tiny fluctuations. The reversed orbit crosses the real axis at $`x_i^{}`$ and approaches the complementary burning line $`\varphi _c^+`$ in the lower halfplane. There, a similar kind of deviation drives it to still another $`\overline{q}(x_i^{\prime \prime };t)`$ to reach again $`\varphi _c^+`$ and so forth and back. By subsequently running through these cycles between the caustics at $`\pm \varphi _c^+`$ a net-motion into the direction of the barrier top may be generated. On the left side of the barrier ($`x_i<0`$) the same kind of scenario exists and at the top $`x=0`$ the burning lines intersect (depending in detail on $`V(q)`$ within $`|q/l|1`$, see fig 1.) so that small deviations in the vicinity of the bottleneck $`x=0`$ may lead from the set of right-barrier paths to that of left-barrier paths and vice versa. This allows for a motion starting in $`x_i>0`$ to eventually reach the range on the opposite side of the barrier. So far the above discussion is restricted to class (c)-orbits with $`ϵ_{\mathrm{re}}=0`$. However, for finite but small $`ϵ_{\mathrm{re}}`$ orbits merge close to the burning lines, and we find basically the same situation. The conclusion is that two real axis paths with TPs at $`q_0`$ and $`q_0`$, respectively, are linked by a sequence of real-time complex plane orbits tied together by small fluctuations near caustic lines. Since this under-barrier-motion is not a purely stationary one obeying (2), but can be seen as nearly stationary as it follows classical orbits most of the time, it describes quasi-stationary fluctuations (QSF). The QSF allow to move from $`q_i>0`$ through the barrier range towards $`q_f<0`$ and this way dominate in absence of true stationary points, $`\delta S[q]=0`$ with $`q(0)=q_i,q(t)=q_f`$, the path integration in $`G(q_f,q_i,t)`$ between the TPs.

## III Extended semiclassical propagator

The action of the QSF can simply be approximated. For a cycle from $`x_i`$ with $`\overline{q}(x_i;t)`$ to $`x_i^{}<x_i`$ with $`\overline{q}(x_i^{};t)`$ in the interval $`\delta t`$ we find with Cauchy’s formula $`S(x_i^{},x_i,\delta t)i|W(x_i^{},x_i)|E\delta t`$ where the short action is $`W(x^{},x)=_x^{}^x𝑑q\sqrt{2m[EV(q)]}`$ and the portion from the phase space deviation along $`\varphi _c^+`$ is negligible. Accordingly, $`S(q_0,q_0,\mathrm{\Delta }t)i|W(q_0,q_0)|E\mathrm{\Delta }t`$ where $`\mathrm{\Delta }t`$ is the time interval spent between the real axis TPs $`q_0,q_0`$. Hence, one arrives at the crucial result that the real-time motion of the QSF gives rise to an imaginary part in the action which is identical to the known instanton or WKB exponent. The full action for a low energy motion from $`q_i>0`$ to $`q_i`$ now consists of two classical real axis segments from $`q_i`$ to $`q_0`$ and from $`q_0`$ to $`q_i`$, respectively, and QSF inbetween, i.e. $`S(q_i,q_i,t)2W(q_0,q_i)+i|W(q_0,q_0)|Et`$. In the semiclassical $`G(q_i,q_i,t)`$ the exponential of this action is accompanied by the contribution of Gaussian fluctuations around the real axis segments.

The most powerful representation of the semiclassical propagator is the so-called Hermann-Kluk propagator (HK) kluk . It has the advantage of being determined by an initial value problem for the classical trajectories, namely,

$$G_{\mathrm{HK}}(q_f,q_i,t)=\frac{dqdp}{2\pi \mathrm{}}h(q_f,q_i,t,p,q)R(p,q,t)\mathrm{e}^{iS(p,q,t)/\mathrm{}}$$
(6)
with the fluctuation prefactor $`R(p,q,t)`$ and an overlap factor $`h(q_f,q_i,t,p,q)=q_f|\gamma (p,q,t)\gamma (p,q)|q_i`$ where

$$x|\gamma (p,q,t)=\left(\frac{\gamma }{\pi }\right)^{1/4}\mathrm{exp}\left\{\frac{\gamma }{2}[xq(t)]^2+\frac{i}{\mathrm{}}p(t)[xq(t)]\right\}$$
(7)
is a Gaussian wave packet centered around the phase space point $`\{p(t),q(t)\}`$. In (6) one runs in the time interval $`t`$ real trajectories from $`\{p,q\}`$ to $`\{p(t),q(t)\}`$ where the contribution of each orbit is weighted according to its action and fluctuation prefactor and the Gaussian overlap of its end-points with those of the propagator.

The usual HK (6) is exact for a pure parabolic barrier and therefore captures tunneling only for harmonic fluctuations around the barrier top but fails for long times heller1 . To overcome this drawback we apply the results developed above and formally split the propagator in phase space: $`G(q_f,q_i,t)=G_>(q_f,q_i,t)+G_<(q_f,q_i,t)`$ where $`G_>`$ \[$`G_<`$\] contains orbits with $`EV_0`$ \[$`E<V_0\delta _{\mathrm{pb}}`$\] and $`G_>`$ coincides with $`G_{\mathrm{HK}}`$ (6). Accordingly, $`G_>(t)`$ describes the time evolution up to moderate times (comprising the parabolic range $`V_0\delta _{\mathrm{pb}}`$ below the top) and $`G_<(t)`$ the long time behavior. Now, while in a strict sense the complex dynamics discussed above is only valid for very low $`E`$, we assume its applicability also for somewhat larger $`E`$ and find with $`q_i`$ and $`q_f`$ on opposite sides of the barrier

$$G_<(q_f,q_i,t)=\underset{E<V_0\delta _{\mathrm{pb}}}{}\frac{dqdp}{2\pi \mathrm{}}h(q_f,q_i,t,p,q)R^<(p,q,t)\mathrm{e}^{iS^<(p,q,t)/\mathrm{}}T(q_0).$$
(8)
Here, an orbit runs from $`\{p,q\}`$ along the real axis to its TP $`\{0,q_0\}`$, jumps to $`\{0,q_0\}`$ to reach $`\{p(t),q(t)\}`$ leading to a fluctuation prefactor $`R_<`$ and action $`S_<`$. The position space jump costs $`T(q_0)=\mathrm{exp}[|W(q_0,q_0)|/\mathrm{}]`$ and eventually $`G_<`$ results from phase space averaging. In (8) we require $`|q_i|/l,|q_f|/l1`$ so that most of the dynamics is spent outside the barrier. Obviously, $`G_<`$ follows not just from switching in the integrand in $`G_>`$ to imaginary times in regions where $`E<V(q)`$.

## IV Applications

The extended HK (eHK) $`G_{\mathrm{eHK}}=G_>+G_<`$ is now employed to scattering in an Eckart barrier $`V_b(q)=V_0/\mathrm{cosh}^2(q/l)`$ that has been of wide use, e.g. as a model for the H+H2 exchange reaction. Since asymptotically $`V(q)`$ drops faster than any power of $`q`$ we have $`\varphi _c^+0`$ and burning lines stretch parallel to the real axis. In fig. 2a the correlation function $`c_{fi}(t)=\psi _f|\mathrm{exp}(iHt/\mathrm{})|\psi _i`$ between two Gaussian wave packets is depicted. Initially, $`\psi _i`$ \[$`\psi _f`$\] is centered to the far right \[far left\] with $`V_0p_i^2/2`$ so that we are indeed in a deep tunneling regime. One clearly sees the exponential drop of $`G_>`$ and the startling accuracy of the eHK over the whole time range. The most sensitive observable for a real-time treatment is the transmission probability $`P(E)`$ calculated by numerically Fourier transform $`c_{fi}(t)`$. Remarkably, we get accurate data also for very low energies (fig. 3) apart from small oscillations typical for real-time calculations tannor . In the moderate energy range $`E/V_0>0.5`$ the “real-time” $`P(E)`$ even improves the uniform WKB result. Here convergence for $`G_<`$ is achieved for roughly the same number of trajectories as in $`G_>`$ (typical number of trajectories for the set of parameters is $`510^4`$) so that in contrast to previous approaches grossmann  an extension of the eHK to two or three dimensional systems, which are of particular interest to study chaotic tunneling, seems feasible.

As an example where a dynamical approach is clearly needed we turn to an Eckart barrier driven by a periodic signal $`V=V_b+qA\mathrm{sin}(\mathrm{\Omega }t)`$ and focus on the range of non-resonant driving and weak to moderate driving amplitudes. In this case already the exact numerics is non-trivial since it is the long time tunneling behavior which is most sensitively affected by the driving and leads to a strong spreading of the wave packet. Typical results for the correlation function $`c_{fi}(t)`$ are shown in fig. 2b. Compared to the static case one sees phase shifted oscillations and a revival type of phenomenon. Semiclassically, both effects originate from an intimate interference of (fast) above-barrier-paths ($`E>V_0`$), which cross the barrier and then are back-scattered, and (slow) driven tunneling orbits ($`E<V_0`$). Even in this time-dependent case the accuracy of the eHK is quite astonishing.

## V Conclusion

To conclude, our findings reveal for the first time how tunneling is encoded in the quantum propagator in terms of classical real-time orbits. In contrast to tunneling in the energy domain, real time barrier penetration cannot be described by individual tunneling orbits. Instead, it must be seen as a diffusion along a certain set of classical paths in the complex plane. This allows for a practical approach for semiclassical wave packet dynamics even in the deep tunneling regime of static and non-resonantly driven scattering processes. Explicit examples have been given for one dimensional cases, but the efficiency of the method suggests that at least two or three dimensional cases may be feasible. The main problem then will be that from a certain TP $`q_0`$ a bunch of TPs on the other side of the barrier can be reached. This proliferation of orbits, however, seems tractable due to the exponential suppression of under barrier motion starting at $`q_0`$ and traveling over large distances. Work in this direction is in progress. Possible applications of our method and extensions of it in physics and chemistry may be e.g. mesoscopic systems in microwave fields or unimolecular reactive scattering.

Financial support by the DFG through SFB276 is gratefully acknowledged.


format: text

text:
# Corrections to Scaling are Large for Droplets in Two-Dimensional Spin Glasses

## Abstract

The energy of a droplet of linear extent $`l`$ in the droplet theory of spin glasses goes as $`l^\theta `$ for large $`l`$. It is shown by numerical studies of large droplets in two dimensional systems that this formula needs to be modified by the addition of a scaling correction $`l^\omega `$ in order to accurately describe droplet energies at the length scales currently probed in numerical simulations. Using this simple modification it is now possible to explain many results which have been found in simulations of three dimensional Ising spin glasses with the droplet model.

There are two rival theories to describe the ordered state of spin glasses: the droplet theory McM ; BM ; FH  and the replica symmetry breaking (RSB) theory of Parisi Parisi ; sgb . The droplet picture assumes that the low-temperature behavior is governed by droplet-like excitations, where excitations of linear spatial extent $`l`$ typically cost an energy of order $`l^\theta `$. Thus in the thermodynamic limit the excitations which flip a finite fraction of the spins cost an infinite amount of energy if $`\theta >0`$. These excitations are expected to be compact and their surface has a fractal dimension $`d_s<d`$, where $`d`$ is the space dimension. Furthermore it is usually assumed that the energy of different types of excitations, e.g. droplets and domain walls, induced by changing the boundary conditions, are described by the same exponent $`\theta `$. On the other hand the RSB theory assumes that there are low-energy droplet excitations in which a finite fraction of the spins in the system are reversed but which only cost a finite amount of energy in the thermodynamic limit. The surface of these excitations is space-filling so that the fractal dimension of their surface $`d_s=d`$ Marinari .

Many numerical studies have been done in an attempt to resolve the controversy, but studies in three 3D  and four dimensional systems 4D  are limited in the range of systems sizes which can be studied. In two dimensions, however, much larger systems can be studied. For this dimension, no spin-glass order exists at any finite temperature SG-kawashima1997 , but there are still puzzling results concerning droplet excitation energies about the groundstate. The energies of cross-system domain walls perturb ; HY  and droplet-like excitations KA ; picco2001  are apparently descibed by different exponents, which is not compatible with the (usual) droplet picture and might be regarded as a vestige of RSB type behavior.

In this letter we study numerically droplets in two-dimensional spin glasses similar to the droplets studied in Ref. KA, , but up to even larger sizes. We show that the numerical results mentioned above can be explained by a simple scaling correction. For larger droplets we find a behavior which is compatible with the droplet prediction, i.e. the exponent $`\theta `$ is the same for droplets and cross-system domain-walls. Smaller droplets appear to be described by an effective exponent $`\theta ^{^{}}`$. We also study how the volume of an average droplet behaves with system size. We find that when considering just smaller system sizes, droplets appear non-compact, but for larger systems, crossover to a compact behavior looks likely. Thus in two-dimensions the droplet picture seems to be entirely valid. Finally, we discuss also some recent numerical results LQ ; KM ; PY  in three dimensions, which also may be explained within the droplet picture using the correction to scaling for droplet energies introduced here.

We consider the Hamiltonian which is usually studied in numerical work of the Ising spin glass model:

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where the sites $`i`$ lie on the sites of a cubic lattice with $`N=L^d`$ sites, $`S_i=\pm 1`$, the $`J_{ij}`$ have a Gaussian distribution of zero mean and unit variance and couple nearest-neighbor sites on the lattice.

It is important to realize what the low-energy excitations look like in spin glasses. Consider a domain wall crossing a system $`M\times M\times L`$. The energy cost of this domain on the droplet picture will be of order $`M^\theta `$ and the domain wall will be fractal with dimension $`d_s`$ and with an area of order $`M^{d_s}`$. In three or more dimensions the interface may have holes through it. Because it is fractal the extent of this wandering is of order $`M`$. The wandering of the interface by an amount of order $`M`$ affects the determination of $`\theta `$, see Ref. Carter, , where it was found that the best results were obtained when $`LM`$ as then the interface is not affected by its interactions with the ends of the system. Suppose now we have two domain walls across the system. Then if their separation is large compared to $`M`$ they will be unaffected by the presence of the other. However, if they are closer together they interfere with each other and their overall energy will be greater than if the other one were absent. This is because domain walls wander to take advantage of weak bonds and since the domain walls cannot cross, the presence of the second domain wall will preclude the first domain wall from cutting through bonds it would have cut in the absence of the second wall. As a consequence of this, domain walls effectively repel each other. If they have a separation of order $`l`$, then the repulsive energy between them would be expected to vary with their separation as a power law, $`l^\omega ^{^{}}`$. No investigations of $`\omega ^{^{}}`$ seem to exist in the literature.

Consider now a droplet of linear extent $`l`$. It too will have a fractal surface described by $`d_s`$. Pictures of large droplets have now appeared KM , and a systematic investigation of them is in HKM . Their fractal nature ensures that in three or more dimensions that they have holes through them, giving them a sponge-like appearance. However, because the surface of the droplet may wander by a distance of order $`l`$ the energy of the droplet $`E`$ will be modified by its wandering and “collisions” with itself to a form which we suppose by analogy with the above is

$$E=Al^\theta +Bl^\omega .$$
(2)
The term $`Bl^\omega `$ is a scaling correction to the form of the droplet energy at large $`l`$ and $`A`$ and $`B`$ are positive constants. We will show that the scaling correction is very significant in two dimensions.

In this dimension $`\theta `$ is very accurately determined as systems of $`480^2`$ can be studied HY . The exponent associated with a single domain wall, that is $`\theta `$, is approximately $`0.29`$ Carter ; HY . However, studies where, say, $`\theta `$ is determined from the effects of thermally excited droplets (such as Monte Carlo simulations of the spin-glass susceptibility) yield an apparent $`\theta ^{^{}}`$ close to $`0.47`$ KA ; picco2001 . This discrepancy has long been a puzzle, and has prompted suggestions that perhaps different exponents describe domain wall energies and droplet energies KA . However, Eq. (2) offers a way of resolving the discrepancy. For small values of $`l`$, droplets would apparently have energies decreasing as $`l^{|\theta ^{^{}}|}`$ when their energies are affected by the correction to scaling term, but at large values of $`l`$ the decrease with $`l`$ will be slower and be as $`l^{|\theta |}`$, as expected from the conventional droplet approach.

We have tested this explicitly by studying larger sizes than before using an exact ground-state (GS) algorithm opt-phys2001 . For the special case of a planar system without magnetic field, e.g. a square lattice with periodic boundary conditions in at most one direction, there are efficient polynomial-time “matching” algorithms bieche1980 . The basic idea is to represent each realization of the disorder by its frustrated plaquettes toulouse1977 . Pairs of frustrated plaquettes are connected by paths in the lattice and the weight of a path is defined by the sum of the absolute values of the coupling constants which are crossed by the path. A ground state corresponds to the set of paths with minimum total weight, such that each frustrated plaquette is connected to exactly one other frustrated plaquette. This is called a minimum-weight perfect matching. The bonds which are crossed by paths connecting the frustrated plaquettes are unsatisfied in the ground state, and all other bonds are satisfied.

We want to mimic the droplet generation used in Ref. KA, , where a center spin was forced to reverse its orientation with respect to its orientation in the GS, the boundary spins were fixed, and a new GS under these constraints was calculated. Note that fixing the absolute orientation of a spin does not allow for a fast polynomial algorithm. Hence, a different approach must be used to study large systems.

Our generation of the droplets works in the following way. For each realization, first a ground state is calculated with free boundary conditions in both directions. Then several “hard” bonds are introduced, i.e. bonds with a high value of the absolute strength (e.g. $`|J_{ij}|=2N\times \mathrm{max}_{i,j}\{J_{ij}\}`$). If the subsystem of hard bonds does not exhibit frustration, no hard bond will be broken, if a new GS is calculated. First, all boundary spins are fixed relative to each other by hard bonds around the border chosen to be compatible with their GS orientations, i.e. the bonds between pairs of boundary spins are replaced. Second, a line of hard bonds is created which runs from the middle of (say) the left border to a pre-chosen center spin, again fixing the pair’s spins in their relative GS orientations. Next, the sign of exactly one hard bond on this line is inverted. Finally, a GS of the modified realization is calculated. With respect to the original GS, the result is a minimum energy droplet fulfilling the constraints that it contains the center spin, does not run beyond the boundary and that it has a surface which runs through the hard bond which has been inverted. The energy of the droplet is defined as the energy of the resulting configuration calculated using the original bond configuration. For each realization, this procedure is iterated over all the bonds which are located on the line from the boundary to the center, when in each case exactly one hard bond is inverted. Among all generated droplets, the one exhibiting the lowest energy is selected.

This procedure is similar, but not equivalent to the method used in Ref. KA, . The difference is that the surface of our droplets crosses the line from the boundary to the center spin exactly once. This means, it cannot cross this line more than once, i.e. it cannot “meander” across this line. But it may meander in all other directions. To reduce the influence of this constraint, we have iterated over all four lines of bonds running from the left, right, top and bottom boundary to the center spin, and selected the minimum energy droplet among all the $`2L2`$ droplets generated in this way. Hence, we have obtained droplets which may meander in any direction, but not in all four directions at the same time.

We have tested explicitly the influence of this constraint by measuring in how many directions each of the lowest energy droplets meanders. For small systems meandering does not occur, but meandering increases with system size. For the largest systems we have considered, about half the droplets do not meander at all, about 33 % meander in one direction, 14 % meander in two directions and three percent in three directions. Hence, for the system sizes studied, four-way meandering droplets probably would occur rarely, if it were allowed by the method. Hence, we believe that our results are not strongly influenced by the constraint imposed by our algorithm. Furthermore, we have also tested droplets generated in a different way, which allow for full meandering, but this approach is restricted to small sizes, see below.

We have generated droplets for system sizes in the range $`6L160`$, each time 5000 independent realizations. The average droplet energy as a function of system size is shown in Fig. 1. For small sizes up to $`L=40`$, the energy can be fitted by an algebraic function $`L^\theta ^{}`$ with $`\theta ^{}=0.45(1)`$, which is compatible with the earlier result KA  which had been interpreted in favor of the droplet exponent being different from the domain-wall exponent. For larger sizes, a clear crossover can be observed, the data is compatible with an exponent $`\theta =0.29`$. A fit of the exponent in $`[80,160]`$ yields $`\theta =0.32(1)`$, which indicates that the thermodynamic limit may not been fully reached. We have also tested the scaling ansatz Eq. (2) by fixing $`\theta =0.29`$, resulting in $`\omega =0.97(5)`$ and a very good quality goodness  $`Q=0.81`$ of the fit. Curiously, the value for $`\omega `$ is close to the value which $`\theta `$ takes in one dimension, $`1`$.

To estimate the influence of the slightly restricted meandering, we have also studied droplets, where in addition to the hard bonds around the boundary, only one additional inverted hard bond was introduced. For each system, this inverted hard bond was iterated over all bonds connecting at least one non-boundary spin. The minimum energy droplet was selected among all droplets containing the pre-selected central spin. Due to the large numerical effort, only sizes up to $`L=30`$ could be studied. When considering up to 1000 realizations, the data was compatible with an algebraic decay with $`\theta ^{}=0.47(1)`$. We have increased the statistics, such that we have treated 80000 realizations for small sizes down to 20000 realizations for $`L=30`$. Now the straight algebraic fit has a very low quality ($`Q=10^{26}`$) compared to the fit including the correction to scaling ($`\omega =2,\theta =0.29,Q=0.61`$).

For further comparison with Ref. KA, , we show in Fig. 2 the average droplet volume and surface as a function of the system size. By fitting to algebraic laws, we find a volume fractal dimension of $`D_v=1.80(1)`$ and a surface fractal dimension of $`d_s=1.10(1)`$, which are both compatible with the results obtained before. Hence, the crossover observed in the energy is not obvious when considering volume and surface and the droplets seem to be non-compact (i.e. $`D_v<2`$). We have tested, whether this result may be a consequence of a correction to scaling, by fitting $`V(L)`$ to a function of the form $`EL^2+FL^{D_v^{}}`$ resulting in $`D_v^{}=1.52(5)`$ and a much better quality of the fit; $`Q=0.59`$ as opposed to $`Q=0.002`$ for the pure algebraic fit. Hence a crossover to compact droplets seems likely. In the figure, both functions would be hard to distinguish, hence we have omitted including the correction to scaling fit. We have also tried a correction to scaling for the surface, since $`d_s`$ around 1.3 in other studies perturb , but using the correction to scaling we could change $`d_s`$ in a wide range without affecting the quality of the fit.

A crossover between large and small $`L`$ dependencies was seen in the work of Middleton M . He studied the scaling of ground-state link overlaps and observed a gradual evolution of the finite-size behavior. Middleton himself attributed this gradual evolution to averages over droplets of all sizes up to the system size $`L`$. This kind of behavior is an additional mechanism, independent of our scaling correction for the droplet energies themselves, to explain why certain averages might be slow to converge to the expected behavior for large $`L`$.

We have seen that in two dimensions that the corrections to scaling are large, so that for medium-size systems they can mask the true asymptotic behavior. A similar effect would be expected to occur in higher dimensions and indeed might already have been observed by Lamarcq et al. LQ  in their study of droplets in three dimensions. They first computed the ground state for both $`N=6^3`$ and $`N=10^3`$ spins. They then chose an arbitrary reference spin and flipped it along with a cluster containing $`V1`$ other spins connected to it. They next minimized the energy of this cluster by exchange Monte Carlo, but with the constraint that the reference spin is held fixed and the cluster was always connected and of size $`V`$. They found the largest extension (mean end-to-end distance) $`l`$ of the cluster and found that for $`V33`$ its energy $`E(l)`$ varied as $`l^\theta ^{^{}}`$, where the exponent $`\theta ^{^{}}`$ was $`0.13\pm 0.02`$. For larger $`V`$ values they found that $`E(l)`$ was increasing rather than decreasing with $`l`$ but they were unsure whether this might not be an artifact of insufficient numbers of Monte Carlo steps. Note that according to Eq. (2), when both $`\theta `$ and $`\omega `$ are small and positive, $`E(l)`$ will have a shallow minimum at some value $`l_0`$, similar to the behavior actually seen in Ref. LQ, .

Furthermore, for droplets whose linear dimensions are similar to the minimum at $`l_0`$, the apparent value of $`\theta `$ would appear to be close to zero. This may be the origin of the behavior in Refs. PY, ; KM, , where simulations of Parisi’s spin overlap distribution function and the link overlap function were only consistent with the droplet picture if $`\theta `$ was taken to be about zero rather than its domain wall value of order 0.20 in three dimensions 3D  and 0.70 in four dimensions 4D .

To summarize: there exists a natural correction to the usual scaling formula for the energy of a droplet. We have explicitly verified its validity by numerical studies of large systems in two dimensions. If a similar correction is valid in three dimensions, it can explain a number of puzzling features which have been seen in simulations. A direct check of the validity of this scaling correction needs numerical studies of larger three-dimensional systems than those currently studied.

###### Acknowledgements.

We would like to thank Alan Bray for many discussions. AKH obtained financial support from the DFG (Deutsche Forschungsgemeinschaft) under grants Ha 3169/1-1 and Zi 209/6-1. The simulations were performed at the “Paderborn Center for Parallel Computing” and the “Gesellschaft f—r Wissenschaftliche Datenverarbeitung” in Gvttingen, both in Germany.


text:
# Open environments for quantum open systems

## I Introduction

There has been much interest in quantum open systems (particularly with applications to quantum decoherence Zurek (1991, ); Giulini et al. (1996) over the past several decades.Alicki

As discussed by Alicki,Alicki  models that have been extensively studied in the literature have arisen from one of two approaches. The first starts with an a priori proposal of a mathematical form for the quantum dynamical semigroup (Alicki’s axiomatic approach), and proceeds with a study of the resulting properties of the dynamics. The second is to study the reduced dynamics of the system of interest which is part of a composite system which also includes environmental degrees of freedom (Alicki’s constructive approach). The models primarily explored for the constructive approach have a set of common features. The system of interest evolves unitarily in isolation. The environment (often an oscillator bath) evolves unitarily in isolation. The interaction between system and environment consists of a term in the composite hamiltonian, so that the composite system evolves unitarily. The initial state of the composite system is an important feature of the particular model being studied and is often but not exclusively taken to be a factor state between an initial state of the system of interest and a thermal equilibrium state for the environment.

The study of open systems is an attempt to describe how “the rest of the universe” influences the system of interest. However, most constructive models only include that part of the rest of the universe with which the rest of the universe directly interacts (through an interaction term in the composite hamiltonian). In order to account for that part of the universe with which the system of interest does not directly interact, we explore properties of models for which the environment degrees of freedom are themselves open systems. In the Section II, we motivate the quantum models by examining a classical system and environment where the environment degrees of freedom are damped oscillators relaxing into thermal equilibrium. In Section III we review some important properties of the model we will use to simulate the open nature of the environment. In Section IV we proceed to use the projection operator method in the weak coupling regime on a composite system where the environment relaxes to thermal equilibrium via a nonunitary evolution of the Lindblad form.Lindblad (1976) We illustrate the results for the particular case of linear coupling to an oscillator bath in Section V. We finish in Section VI with a discussion of our results.

## II Classical Open System

In this section we examine a classical open system with a damped environment to gain insight and anticipate properties of analogous quantum models. We are essentially following ZwanzigZwanzig (1973) with the addition of dissipation and thermalization to the dynamics of the environment.

The system is described in terms of its coordinate $`x`$. The environment consists of a bath of independent oscillators with coordinates $`q_\mu `$. The composite system equations of motion are given by

$$m\ddot{x}=\frac{U(x)}{x}+\underset{\mu }{}m_\mu \omega _\mu ^2[q_\mu a_\mu (x)]\frac{a__\mu (x)}{x}$$
(1)
for the system coordinate, and

$`m_\mu \ddot{q_\mu }`$ $`=`$ $`m_\mu \omega _\mu ^2q_\mu 2m_\mu \gamma _\mu \dot{q}_\mu `$ (2)
$`+m_\mu \omega _\mu ^2a_\mu (x)+F_\mu (t)`$
for each environment degree of freedom. These equations of motion correspond to the addition of dampening (with friction coefficients $`\eta _\mu =2m_\mu \gamma _\mu `$ ) and thermal noise terms $`F_\mu `$ to the environment portion of the closed system dynamics obtained from the composite system Lagrangian:

$$L=\frac{1}{2}m\dot{x}^2U(x)+\underset{\mu }{}\frac{m_\mu }{2}\{\dot{q}_\mu ^2\omega _\mu ^2[q_\mu a_\mu (x)]^2\}$$
(3)
Each oscillator is taken to be independent but driven to the same temperature by dissipation and white noise so that the following fluctuation-dissipation relation holds:

$$F_\mu (t)F_\nu (s)=4\gamma _\mu m_\mu k_BT\delta _{\mu \nu }\delta (ts).$$
(4)
The Kronecker delta indicates the oscillators in the bath are driven by independent noise terms. The Dirac delta function indicates that the driving forces are white noise.

With $`\mathrm{\Omega }_\mu \sqrt{\omega _\mu ^2\gamma _\mu ^2}`$ taken to be real (the bath oscillators are underdamped), the formal solutions for Eq. (2) can be written in terms of the oscillator’s initial position $`q_{\mu 0}`$ and velocity $`\dot{q}_{\mu 0}`$ as

$`q_\mu (t)=`$ $`[q_{\mu 0}\mathrm{cos}\mathrm{\Omega }_\mu t+{\displaystyle \frac{1}{\mathrm{\Omega }_\mu }}(\dot{q}_{\mu 0}+\gamma _\mu q_\mu )\mathrm{sin}\mathrm{\Omega }_\mu t]e^{\gamma _\mu t}+{\displaystyle _0^t}{\displaystyle \frac{1}{m_\mu \mathrm{\Omega }_\mu }}\mathrm{sin}\mathrm{\Omega }_\mu (ts)e^{\gamma _\mu (ts)}F_\mu (s)𝑑s`$ (5)
$`a_\mu (x(t)){\displaystyle \frac{1}{\mathrm{\Omega }_\mu }}a_\mu (x(0))[\mathrm{\Omega }_\mu \mathrm{cos}\mathrm{\Omega }_\mu t+\gamma _\mu \mathrm{sin}\mathrm{\Omega }_\mu t]e^{\gamma _\mu t}`$
$`{\displaystyle _0^t}{\displaystyle \frac{1}{m_\mu \mathrm{\Omega }_\mu }}[\mathrm{\Omega }_\mu \mathrm{cos}\mathrm{\Omega }_\mu s+\gamma _\mu \mathrm{sin}\mathrm{\Omega }_\mu s]e^{\gamma _\mu s}{\displaystyle \frac{a_\mu (x(ts))}{x}}\dot{x}(ts)𝑑s.`$
Upon substitution of Eq. (5) into Eq. (1) and rearranging, the equation of motion for the system can be written

$`m\ddot{x}`$ $`=`$ $`{\displaystyle \frac{U(x)}{x}}{\displaystyle _0^t}\eta (x(t),x(s);t,s)\dot{x}(s)𝑑s`$ (6)
$`+F_{E1}(t)+F_{E2}(t)+F_{E3}(t)`$
where the noise terms are given by

$`F_{E1}(t)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}m_\mu \omega _\mu ^2[q_{\mu 0}a_\mu (x(0))][\mathrm{cos}\mathrm{\Omega }_\mu t+{\displaystyle \frac{\mathrm{\Omega }_\mu }{\gamma _\mu }}\mathrm{sin}\mathrm{\Omega }_\mu t]e^{\gamma _\mu t},`$ (7a)
$`F_{E2}(t)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}m_\mu \omega _\mu ^2\dot{q}_{\mu 0}[\mathrm{cos}\mathrm{\Omega }_\mu t+{\displaystyle \frac{\mathrm{sin}\mathrm{\Omega }_\mu t}{\mathrm{\Omega }_\mu }}]e^{\gamma _\mu t},`$ (7b)
$`F_{E3}(t)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle \frac{\omega _\mu ^2}{\mathrm{\Omega }_\mu }}{\displaystyle _0^t}\mathrm{sin}\mathrm{\Omega }_\mu (ts)e^{\gamma _\mu (ts)}F_\mu (s)𝑑s{\displaystyle \frac{a_\mu (x(t))}{x}},`$ (7c)
and the dissipation kernel is given by

$`\eta (x(t),x(s);t,s)={\displaystyle \underset{\mu }{}}m_\mu \omega _\mu ^2\{[\mathrm{cos}\mathrm{\Omega }_\mu (ts)`$
$`+{\displaystyle \frac{\mathrm{\Omega }_\mu }{\gamma _\mu }}\mathrm{sin}\mathrm{\Omega }_\mu (ts)]e^{\gamma _\mu (ts)}{\displaystyle \frac{a_\mu (x(s))}{x}}{\displaystyle \frac{a_\mu (x(t))}{x}}\}.`$ (8)
In the past, the properties of noise terms like $`F_{E1}(t)`$ and $`F_{E2}(t)`$ have generally been extracted from the statistical distributions of the initial conditions. With the introduction of dissipation in the environment, those noise terms can readily be seen to be transient terms on the time scales of the environment (as determined by $`\gamma _\mu `$). Thus the details of the initial state of the environment are not important to the long term system dynamics.

The correlations $`F_{E3}(t)F_{E3}(\tau )`$ of the remaining noise term depend upon the correlations of the individual oscillators’ noise terms via Eq. (4), and can be written

$`F_{E3}(x,t)F_{E3}(x^{},\tau )=k_BT\eta (x,x^{};t,\tau )`$
$`+{\displaystyle \frac{1}{4\gamma _\mu \omega _\mu ^2}}[\gamma _\mu ^2\mathrm{cos}\mathrm{\Omega }_\mu (t+\tau )\omega _\mu ^2\mathrm{cos}\mathrm{\Omega }_\mu (t\tau )`$
$`\omega _\mu \gamma _\mu \mathrm{sin}\mathrm{\Omega }_\mu (t+\tau )]e^{\gamma _\mu (t+\tau )}{\displaystyle \frac{a_\mu (x))}{x}}{\displaystyle \frac{a_\mu (x^{})}{x^{}}}.`$ (9)
The first term on the right hand side of Eq. (9) is the long term correlation function of the effective noise. The second term contains a factor of $`e^{\gamma _\mu (t+\tau )}`$, and thus is transient for long timescales.

Since $`\eta (x,x^{};t,\tau )`$ is a narrow function of $`t\tau `$, we can make a Markov approximation in Eq. (6):

$$m\ddot{x}=\frac{U(x)}{x}\overline{\eta }(x,x)\dot{x}+F_s(x,t),$$
(10)
where

$$\overline{\eta }(x,x^{})=\underset{\mu }{}2m_\mu \gamma _\mu \frac{a_\mu (x))}{x}\frac{a_\mu (x^{})}{x^{}}.$$
(11)
The corresponding fluctuation-dissipation relation is

$$F_s(x,t)F_s(x^{},\tau )=2k_BT\eta (x,x^{})\delta (t\tau ).$$
(12)
The role of the spatial correlations of the noise in quantum decoherence and the lack of importance of those correlations in classical phenomena has been discussed elsewhere.Gallis (1993, 1994)

Thus, adding “fast” thermal relaxation to the environment of the system of interest lead to Markovian equations of motion with the usual fluctuation-dissipation relations. Memory effects due to the details of the initial environment state (including possible correlations with the initial system state) are erased on the environment relaxation time scales. This result provides excellent motivation for exploring similar models in the quantum mechanical domain.

## III  Properties of quantum open environment

As with the classical model, in order to account for the open nature of the environment we need to incorporate modifications to the dynamics of the environment degrees of freedom. The environment is modelled as a set of independent oscillators whose evolution is governed by a Markovian master equation:

$$\frac{\rho }{t}=L[\rho ].$$
(13)
The generator $`L`$ of the evolution is taken to be of the Lindblad formLindblad (1976):

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`L[\rho ]`$ (14)
$`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}[H,\rho ]+{\displaystyle \frac{1}{2\mathrm{}}}{\displaystyle \underset{\mu }{}}[V_\mu \rho ,V_\mu ^{}]+[V_\mu ,\rho V_\mu ^{}].`$
Formally, the solution to Eq. (14) is

$$\rho (t)=\mathrm{\Lambda }(t)[\rho ]=e^{Lt}[\rho ],$$
(15)
expectation values are defined via the trace operation, allowing the definition of the adjoint representation of the evolution operator $`\mathrm{\Lambda }^{}(t)[O]`$ via

$$\mathrm{Tr}[\rho \mathrm{\Lambda }^{}(t)[O]]=\mathrm{Tr}[\mathrm{\Lambda }(t)[\rho ]O],$$
(16)
and for the generator $`L^{}`$

$$\mathrm{Tr}[\rho L^{}[O]]=\mathrm{Tr}[L[\rho ]O].$$
(17)
The adjoint representation of $`L`$ (i.e. Heisenberg picture) is given by

$$L^{}[O]=\frac{1}{i\mathrm{}}[H,O]+\frac{1}{2\mathrm{}}\underset{\mu }{}V_\mu ^{}[O,V_\mu ]+[V_\mu ^{},O]V_\mu .$$
(18)
To model the relaxation of the environment, we will use a subset of a family of master equations that have been studied extensively in the literature Sandulescu and Scutaru (1987); Isar et al. (1994); Isar (1999); A. Isar and Scheid (1999) where the Lindblad form of the master equation can be rewritten as:

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}[H_0,\rho ]{\displaystyle \frac{i}{2\mathrm{}}}(\lambda +\mu )[q,\{p,\rho \}]`$ (19)
$`+{\displaystyle \frac{i}{2\mathrm{}}}(\lambda \mu )[p,\{q,\rho \}]{\displaystyle \frac{D_{pp}}{\mathrm{}^2}}[q,[q,\rho ]]`$
$`{\displaystyle \frac{D_{qq}}{\mathrm{}^2}}[p,[p,\rho ]]`$
$`+{\displaystyle \frac{D_{pq}}{\mathrm{}^2}}([q,[p,\rho ]]+[p,[q,\rho ]]).`$
The details of the model are determined by the specification of the diffusion coefficients $`D_{qq}`$, $`D_{pp}`$ and $`D_{pq}`$ and damping constants $`\mu `$ and $`\lambda `$, subject to the constraints:Sandulescu and Scutaru (1987)

$`D_{qq}`$ $`>`$ $`0`$
$`D_{pp}`$ $`>`$ $`0`$
$`D_{qq}D_{pp}D_{pq}^2`$ $``$ $`({\displaystyle \frac{\lambda \mathrm{}}{2}})^2.`$ (20)
The nominal hamiltonian is

$$H_0=\frac{p^2}{2m}+\frac{m\omega _0^2}{2}q^2$$
(21)
$`\stackrel{~}{\rho }`$ is taken to be the (unique) asymptotic state of $`\mathrm{\Lambda }(t)`$, that is

$$\mathrm{\Lambda }_{\mathrm{}}[\rho ]\underset{t\mathrm{}}{lim}\mathrm{\Lambda }(t)[\rho ]=\stackrel{~}{\rho }$$
(22)
for any initial environment state $`\rho `$. Since $`\stackrel{~}{\rho }`$ is a stationary state of $`L`$

$$L[\stackrel{~}{\rho }]=0.$$
(23)
The asymptotic behavior of the environment is encapsulated in $`\mathrm{\Lambda }_{\mathrm{}}`$, so that

$$\mathrm{\Lambda }_{\mathrm{}}[\rho ]=\stackrel{~}{\rho }^{(E)}$$
(24)
for any state $`\rho `$. Although $`\mathrm{\Lambda }_{\mathrm{}}`$ is defined only on the space of density operators (including pure states), we will need to extend its domain. If $`\{|e_k\}`$ is a basis for the environment Hilbert space, then

$$\mathrm{\Lambda }_{\mathrm{}}[|e_ke_k|]=\stackrel{~}{\rho }$$
(25)
for any $`k`$. Furthermore, if

$`|\psi _{nm}{\displaystyle \frac{1}{\sqrt{2}}}(|e_n+|e_m)`$
$`|\chi _{nm}{\displaystyle \frac{1}{\sqrt{2}}}(|e_n+i|e_m)`$ (26)
then

$`\mathrm{\Lambda }_{\mathrm{}}[|\psi _{nm}\psi _{nm}|]`$ $`=`$ $`\stackrel{~}{\rho }`$ (27)
$`=`$ $`{\displaystyle \frac{1}{2}}(\mathrm{\Lambda }_{\mathrm{}}[|e_ne_n|]+\mathrm{\Lambda }_{\mathrm{}}[|e_me_m|]`$
$`+\mathrm{\Lambda }_{\mathrm{}}[|e_ne_m|]+\mathrm{\Lambda }_{\mathrm{}}[|e_me_n|])`$
$`=`$ $`\stackrel{~}{\rho }+{\displaystyle \frac{1}{2}}(\mathrm{\Lambda }_{\mathrm{}}[|e_ne_m|]`$
$`+\mathrm{\Lambda }_{\mathrm{}}[|e_me_n|]),`$
and

$`\mathrm{\Lambda }_{\mathrm{}}[|\chi _{nm}\chi _{nm}|]`$ $`=`$ $`\stackrel{~}{\rho }`$ (28)
$`=`$ $`{\displaystyle \frac{1}{2}}(\mathrm{\Lambda }_{\mathrm{}}[|e_ne_n|]+\mathrm{\Lambda }_{\mathrm{}}[|e_me_m|]`$
$`+i\mathrm{\Lambda }_{\mathrm{}}[|e_ne_m|]i\mathrm{\Lambda }_{\mathrm{}}[|e_me_n|])`$
$`=`$ $`\stackrel{~}{\rho }+i{\displaystyle \frac{1}{2}}(\mathrm{\Lambda }_{\mathrm{}}[|e_ne_m|]`$
$`i\mathrm{\Lambda }_{\mathrm{}}[|e_me_n|]),`$
Which implies

$$\mathrm{\Lambda }_{\mathrm{}}[|e_ne_m|]=\delta _{nm}\stackrel{~}{\rho }.$$
(29)
For an arbitrary operator $`O`$

$`\mathrm{\Lambda }_{\mathrm{}}[O]`$ $`=`$ $`{\displaystyle \underset{nm}{}}O_{nm}\mathrm{\Lambda }_{\mathrm{}}[|e_ne_m|]`$ (30)
$`=`$ $`{\displaystyle \underset{m}{}}O_{mm}\stackrel{~}{\rho }`$
$`=`$ $`\mathrm{Tr}[O]\stackrel{~}{\rho }`$
There is in general no guarantee that an asymptotic state exists for arbitrary $`L`$, that is, not all choices for the parameters $`D_{qq}`$, $`D_{pp}`$, etc. will be appropriate to model the a system dynamically relaxing to equilibrium. For example, if the parameters satisfy

$$D_{pp}=\frac{\lambda +\mu }{2}\mathrm{}m\omega \mathrm{coth}\frac{\mathrm{}\omega }{2k_BT},$$
(31a)
$$D_{qq}=\frac{\lambda \mu }{2}\frac{\mathrm{}}{m\omega }\mathrm{coth}\frac{\mathrm{}\omega }{2k_BT},$$
(31b)
$$D_{pq}=0,$$
(31c)
then the Gibbs state is the asymptotic state.Sandulescu and Scutaru (1987); Isar et al. (1994) On the other hand, if

$$D_{qq}=\frac{\mathrm{}\lambda }{2m\mathrm{\Omega }},$$
(32a)
$$D_{pp}=\frac{\mathrm{}\lambda m\omega ^2}{2\mathrm{\Omega }},$$
(32b)
$$D_{pq}=\frac{\mathrm{}\lambda \mu }{2\mathrm{\Omega }},$$
(32c)
(where $`\mathrm{\Omega }^2=\omega ^2\mu ^2`$) there can be persistent pure states.Isar (1999); A. Isar and Scheid (1999) Since we are considering the primary effect of the openness of the environment to be effectiveness which relax the environment towards an equilibrium state, we would only consider those choices of parameters for which there is a unique asymptotic state.

Sandulescu and Scutaru have determined the time dependence and asymptotic behavior of various moments of $`p`$ and $`q`$, which will be useful in later calculations. The evolution of the first order moments is given by

$$\frac{q}{t}=(\lambda \mu )q+\frac{1}{m}p$$
(33a)
$$\frac{p}{t}=(\lambda +\mu )pm\omega ^2p,$$
(33b)
which illustrates the role of the dissipation coefficients $`\lambda `$ and $`\mu `$. With

$$S_{pq}\frac{1}{2}\{p,q\},$$
(34)
the evolution of the second order moments is given by

$$\frac{p^2}{t}=2m\omega ^2S_{pq}2(\lambda +\mu )p^2+2Dpp$$
(35a)
$$\frac{q^2}{t}=\frac{2}{m}S_{pq}2(\lambda \mu )q^2+2Dqq$$
(35b)
$$\frac{S_{pq}}{t}=\frac{1}{m}p^2m\omega ^2q^22\lambda S_{pq}+2Dpq$$
(35c)
which illustrates the role of the diffusion coefficients $`D_{pp}`$, $`D_{qq}`$ and $`D_{pq}`$.

The asymptotic first order moments are

$$p_{\mathrm{}}=0$$
(36a)
$$q_{\mathrm{}}=0.$$
(36b)
For the second order moments we have

$`p^2_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{2\lambda (\lambda ^2+\omega ^2\mu ^2)}}\{m^2\omega ^4D_{qq}`$ (37)
$`+[2\lambda (\lambda \mu )+\omega ^2]D_{pp}`$
$`2m\omega ^2[\lambda \mu ]D_{pq}\},`$
$`q^2_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{2(m\omega )^2\lambda (\lambda ^2+\omega ^2\mu ^2)}}[(m\omega )^2\omega ^2D_{qq}`$ (38)
$`+\omega ^2D_{pp}+2m\omega ^2(\lambda +\mu )D_{pq}],`$
and

$`S_{pq}_{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{p,q\}_{\mathrm{}}`$ (39)
$`=`$ $`{\displaystyle \frac{1}{2(m\lambda (\lambda ^2+\omega ^2\mu ^2)}}[(\lambda +\mu )(m\omega )^2D_{qq}`$
$`+(\lambda \mu )D_{pp}+2m(\lambda ^2\mu ^2)D_{pq}].`$
We will also need the time dependence of $`q(t)`$ in the Heisenberg picture. This can be readily extracted from Sandulescu and Scutaru ’s results for the time dependence for the first order moments.Sandulescu and Scutaru (1987) The result is

$`q_H(t)`$ $`=`$ $`\mathrm{\Lambda }(t)^{}[q]`$ (40)
$`=`$ $`q(\mathrm{cos}\mathrm{\Omega }t+{\displaystyle \frac{\mu }{\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }t)e^{\lambda t}`$
$`+`$ $`{\displaystyle \frac{p}{m\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }te^{\lambda t},`$
with

$$\mathrm{\Omega }=\sqrt{\omega ^2\mu ^2}$$
(41)
and $`\mathrm{\Omega }`$ taken to be real (the oscillators are underdamped).

The specific environment we will employ consists of a set of independent oscillators, each subject to evolution of the form Eq. (19), with position operators $`\{q_i\}`$ and associated parameters $`\{\mu _i\}`$, $`\{\lambda _i\}`$, etc. Thus the Hilbert space of the environment is actually the tensor product of the Hilbert space corresponding to each environment degree of freedom. The interaction between the environment and system is accounted for by adding an interaction term to the composite system’s hamiltonian:

$$U_I=\underset{n}{}V_n^{(S)}q_n^{(E)}.$$
(42)
We will need the correlation functions $`q_n(t)q_m`$ in the next section. Using Eq. (40), these correlations can be written in terms of the asymptotic correlations as

$`q_n(t)q_m`$ $`=`$ $`\mathrm{Tr}[\mathrm{\Lambda }^{}(t)[q_n]q_m\stackrel{~}{\rho }]`$ (43)
$`=`$ $`\delta _{nm}[q_n^2_{\mathrm{}}(\mathrm{cos}\mathrm{\Omega }t+{\displaystyle \frac{\mu }{\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }t)e^{\lambda t}`$
$`+`$ $`{\displaystyle \frac{2S_{pqn}_{\mathrm{}}i\mathrm{}}{2m\mathrm{\Omega }}}\mathrm{sin}\mathrm{\Omega }te^{\lambda t}],`$
where we have made use of the relation

$`p_nq_n`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\{p_n,q_n\}+[p_n,q_n])`$ (44)
$`=`$ $`S_{pqn}{\displaystyle \frac{i\mathrm{}}{2}}.`$
## IV  Quantum master equation from an open environment

To construct the master equation, we will use the projection operator method in the weak coupling regime, largely following the presentation of Alicki and Lendi.Alicki and Lendi (1987) The composite system (the system of interest plus its environment) is taken to evolve according to Eq. (13). The generator of this evolution is the combination of the nominal dynamics and an interaction:

$$L=L_0+L_I,$$
(45)
where the system and environment dynamics each contribute separately:

$$L_0=L_S+L_E.$$
(46)
The generators for the system and environment nominal evolutions act on on the corresponding subspaces, so that

$$L_S[A^{(S)}B^{(E)}]=L_S^{(S)}[A^{(S)}]B^{(E)},$$
(47)
and

$$L_E[A^{(S)}B^{(E)}]=A^{(S)}L_E^{(E)}[B^{(E)}]$$
(48)
for all $`A^{(S)}`$ and $`B^{(E)}`$. $`L_S^{(S)}`$ is taken to be unitary, and $`L_E^{(E)}`$ is taken to be of the Lindblad form. It is useful to note that necessarily $`L_S`$ and $`L_E`$ commute.

The projection operator is defined in terms of a partial trace:

$$P_0[O]=\mathrm{Tr}_E[O]\stackrel{~}{\rho }^{(E)}.$$
(49)
We wish to establish some important relations between $`P_0`$ and the generators $`L_S`$ and $`L_E`$.

$`L_EP_0[O]`$ $`=`$ $`L_E[\mathrm{Tr}_E[O]\stackrel{~}{\rho }^{(E)}]`$ (50)
$`=`$ $`\mathrm{Tr}_E[O]L_E^{(E)}[\stackrel{~}{\rho }^{(E)}]=0`$
for any operator $`O`$ so that $`L_EP_0=0`$. $`L_E^{(E)}`$ generates trace preserving evolution so that

$$\mathrm{Tr}_E[L_E^{(E)}[O^{(E)}]]=0.$$
(51)
for all $`O^{(E)}`$. Thus

$`P_0L_E[A^{(S)}B^{(E)}]`$ $`=`$ $`P_0[AL_E^{(E)}[B^{(E)}]`$ (52)
$`=`$ $`(\mathrm{Tr}_E[L_E^{(E)}[B^{(E)}]])A^{(S)}\stackrel{~}{\rho }^{(E)}`$
$`=`$ $`0,`$
for all $`A^{(S)}`$ and $`B^{(E)}`$ so that in general $`P_0L_E=0`$.

$`P_0L_S[A^{(S)}B^{(E)}]`$ $`=`$ $`P_0[L_S^{(S)}[A^{(S)}]B^{(E)}`$ (53)
$`=`$ $`L_S^{(S)}[A^{(S)}]\stackrel{~}{\rho }^{(E)}\mathrm{Tr}^{(E)}[B^{(E)}]`$
$`=`$ $`L_S[[A^{(S)}]\stackrel{~}{\rho }^{(E)}\mathrm{Tr}^{(E)}[B^{(E)}]]`$
$`=`$ $`L_SP_0[A^{(S)}B^{(E)}]`$
for all $`A^{(S)}`$ and $`B^{(E)}`$ so that in general $`P_0L_S=L_SP_0`$. Defining a second projection $`P_11P_0`$, it is easy to see that $`P_1L_S=L_SP_1`$. In order to focus on the reduced dynamics of the system of interest, we study the dynamics of $`P_0\rho `$ in certain approximations. Applying the projectors to the master equation for the composite system and making use of the idempotent property of projections produces the following equations:

$$\frac{P_0\rho }{t}=P_0LP_0P_0\rho +P_0LP_1P_1\rho ,$$
(54a)
$$\frac{P_1\rho }{t}=P_1LP_1P_1\rho +P_1LP_0P_0\rho .$$
(54b)
Eq. (54b) can be formally integrated and substituted into Eq. (54a), so that

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`P_0LP_0P_0\rho (t)`$ (55)
$`+P_0LP_1[e^{P_1LP_1t}P_1\rho (0)`$
$`+{\displaystyle _0^t}e^{P_1LP_1(ts)}P_1LP_0P_0\rho (s)]ds.`$
Using the relations between the projectors and generators discussed above, we have

$`P_0LP_1`$ $`=`$ $`P_0L_IP_1,`$
$`P_1LP_0`$ $`=`$ $`P_1L_IP_0,`$
$`P_0LP_0`$ $`=`$ $`L_SP_0+P_0L_IP_0,`$ (56)
so that Eq. (55) can be written

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`(L_S+P_0L_IP_0)P_0\rho (t)+P_0L_IP_1e^{Lt}P_1\rho (0)`$ (57)
$`+{\displaystyle _0^t}P_0L_IP_1e^{L(ts)}P_1L_IP_0P_0\rho (s)𝑑s.`$
This result is exact, within the constraints placed on the model so far.

The interaction, as specified by Eq. (42), can be written

$`L_I[O]`$ $`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}[U_I,O]`$ (58)
$`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \underset{n}{}}[V_n^{(S)}q_n^{(E)},O].`$
This form, along with the first order moments of the environment asymptotic state, simplifies the first term of the right hand side of Eq. (57). Specifically, for any operator $`O`$,

$`P_0L_IP_0[O]`$ $`=`$ $`\mathrm{Tr}_E[{\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \underset{n}{}}[V_n^{(S)}q_n^{(E)},\mathrm{Tr}_E[O]\stackrel{~}{\rho }^{(E)}]]\stackrel{~}{\rho }^{(E)}`$ (59)
$`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \underset{n}{}}\mathrm{Tr}^{(E)}[q_n^{(E)}\stackrel{~}{\rho }^{(E)}][V_n^{(S)},\mathrm{Tr}_E[O]]\stackrel{~}{\rho }^{(E)}={\displaystyle \frac{1}{i\mathrm{}}}{\displaystyle \underset{n}{}}q_n^{(E)}_{\mathrm{}}[V_n^{(S)},\mathrm{Tr}_E[O]]\stackrel{~}{\rho }^{(E)}.`$
From Eq. (36) each term in the sum has a factor $`q_n^{(E)}_{\mathrm{}}=0`$ so that

$$P_0L_IP_0=0.$$
(60)
Furthermore

$$P_0L_IP_1=P_0L_I(1P_0)=P_0L_I$$
(61)
and

$$P_1L_IP_0=(1P_0)L_IP_0=L_IP_0,$$
(62)
so that Eq. (57) can be written

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`L_SP_0\rho (t)`$ (63)
$`+`$ $`P_0L_IP_1e^{Lt}P_1\rho (0)`$
$`+`$ $`{\displaystyle _0^t}P_0L_Ie^{L(ts)}L_IP_0P_0\rho (s)ds].`$
The second term on the right hand side of Eq. (63)

$$P_0L_IP_1e^{Lt}P_1\rho (0)$$
(64)
represents a transient term which depends upon the initial state of the composite system. For most constructive models the initial state is taken to be a factored state of an arbitrary system state and the environment in its asymptotic state:

$$\rho (0)=\rho ^{(S)}(0)\stackrel{~}{\rho }^{(E)}.$$
(65)
For this type of initial condition $`\rho (0)=P_0\rho (0)`$ so that

$$P_1\rho (0)=0.$$
(66)
However, because we have added relaxation to the dynamics of the environment, this type of factoring assumption is not necessary. We will be taking the weak coupling limit and so we will show that Eq. (64) is (approximately) zero to the lowest nonvanishing order of the interaction in the remaining terms of Eq. (57), which turns out to be second order. Formally we can write

$$e^{Lt}=e^{L_0t}+_0^te^{L(ts)}L_Ie^{L_0s}𝑑s.$$
(67)
Substituting Eq. (67) into Eq. (64) and keeping only second order, we have the approximation

$`P_0L_IP_1e^{Lt}P_1\rho (0)P_0L_IP_1e^{L_0t}P_1\rho (0)`$
$`+P_0L_IP_1{\displaystyle _0^t}e^{L_0(ts)}L_Ie^{L_0s}𝑑sP_1\rho (0).`$ (68)
We are interested in timescales $`t`$ which are assumed to be much longer than the relaxation timescales of the environment, so that

$`e^{L_0t}`$ $`=`$ $`e^{L_St}e^{L_Et}`$ (69)
$``$ $`e^{L_St}\mathrm{\Lambda }_E\mathrm{}.`$
Similarly, for the integral from $`t`$ to $`s`$, either $`t`$ or $`ts`$ (or both) is large compared to the relaxation timescales of the environment, so that either

$$e^{L_0(ts)}e^{L_S(ts)}\mathrm{\Lambda }_E\mathrm{}$$
(70)
or

$$e^{L_0s}e^{L_Ss}\mathrm{\Lambda }_E\mathrm{}$$
(71)
or both. From Eq. (30), for an operator $`O`$ decomposed in terms of the basis for the environment $`\{|e_k\}`$ and the basis for the system $`\{|\varphi _\nu \}`$ we have

$`\mathrm{\Lambda }_E\mathrm{}[O]`$ $`=`$ $`O_{\alpha m,\beta n}\mathrm{\Lambda }_E\mathrm{}[|\varphi _\alpha |e_m\varphi _\alpha |e_m|]`$ (72)
$`=`$ $`O_{\alpha m,\beta n}|\varphi _\alpha \varphi _\alpha |\mathrm{\Lambda }_E\mathrm{}^{(E)}[|e_me_m|]`$
$`=`$ $`O_{\alpha m,\beta n}|\varphi _\alpha \varphi _\alpha |\delta _{mn}\stackrel{~}{\rho }^{(E)}`$
$`=`$ $`\mathrm{Tr}_E[O]\stackrel{~}{\rho }^{(E)}=P_0O.`$
When Eq. (72) is applied to Eq. (68) with either Eq. (72) or Eq. (72) applying, then all terms in Eq. (68) end up with factors of either $`P_0P_1`$ or $`P_1P_0`$, both of which are $`0`$. Using this result and keeping only up to second order in the interaction term, the master equation Eq. (63) becomes

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`L_SP_0\rho (t)`$ (73)
$`+{\displaystyle _0^t}P_0L_Ie^{L_0(ts)}L_IP_0P_0\rho (s)𝑑s.`$
We can rewrite Eq. (73)using Eq. (61) and Eq. (62) to get

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`L_SP_0\rho (t)`$ (74)
$`+`$ $`{\displaystyle _0^t}P_0L_IP_1e^{L_0(ts)}P_1L_IP_0P_0\rho (s)𝑑s.`$
In the integrand in Eq. (74) we see the factor $`P_1e^{L_0(ts)}P_1`$ which (following the discussion above) will be zero for $`ts`$ longer than environment relaxation timescales so that the primary contribution to the integral is for $`st`$. This yields the simplest Markovian master equation we will extract from the model:

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`L_SP_0\rho (t)`$ (75)
$`+P_0L_IL_IP_0\rho (t).`$
Using the particular form of the interaction given by Eq. (58) and noting that the independence of the oscillators comprising the environment implies

$$q_nq_m_{\mathrm{}}=\delta _{nm}q_n^2_{\mathrm{}}$$
(76)
allows us to rewrite Eq. (58) in a final form:

$`{\displaystyle \frac{P_0\rho }{t}}`$ $`=`$ $`L_SP_0\rho (t)`$ (77)
$`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle \underset{n}{}}q_n^2_{\mathrm{}}[V_n,[V_n,P_0\rho (t)]].`$
Eq. (77) represents the added effect of environment induced noise on the system’s dynamics which is responsible for phenomena such as quantum decoherence. However, additional effects such as dissipation are not present and will require a more careful handling of the Markov approximations.

To reconsider the Markov approximations, we return to Eq. (73). The unitarity of the isolated system’s evolution implies that

$$e^{L_St}[AB]=e^{L_St}[A]e^{L_St}[B].$$
(78)
The naive Markov approximation is introduced into Eq. (73) by examining

$`e^{L_0(ts)}`$ $`L_I`$ $`P_0\rho (s)=e^{L_E(ts)}e^{L_S(ts)}L_IP_0\rho (s)`$ (79)
$`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}e^{L_E(ts)}e^{L_S(ts)}[U_I,P_0\rho (s)]`$
$`=`$ $`{\displaystyle \frac{1}{i\mathrm{}}}e^{L_E(ts)}[(e^{L_S(ts)}U_I),(e^{L_S(ts)}P_0\rho (s))]`$
$``$ $`{\displaystyle \frac{1}{i\mathrm{}}}e^{L_E(ts)}[(e^{L_S(ts)}U_I),P_0\rho (t)]`$
The Markovian master equation can be written as

$$\frac{P_0\rho }{t}=L_SP_0\rho (t)+KP_0\rho (t)$$
(80)
where

$$K[P_0\rho ]=_0^{\mathrm{}}P_0L_Ie^{L_0s}L_IP_0P_0\rho 𝑑s,$$
(81)
which becomes for our particular model:

$`K[P_0\rho (t)]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle _0^{\mathrm{}}}P_0[U_I,e^{L_E(s)}[(e^{L_S(s)}U_I),P_0\rho (t)]ds`$ (82)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle _0^{\mathrm{}}}𝑑s\{{\displaystyle \underset{nm}{}}\mathrm{Tr}_E[[V_n^{(S)}q_n^{(E)},e^{L_E(s)}[(e^{L_S(s)}(V_m^{(S)}q_m^{(E)})),P_0\rho (t)]]]\stackrel{~}{\rho }^{(E)}\}.`$
Since $`L_E`$ acts only on the environment, we can make use of the cyclic property of the trace over the environment to rewrite Eq. (82) as

$`K[P_0\rho (t)]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle _0^{\mathrm{}}}ds\{{\displaystyle \underset{nm}{}}\mathrm{Tr}_E[[V_n^{(S)}q_n^{(E)}(s),`$ (83)
$`[V_m^{(S)}(s)q_m^{(E)},P_0\rho (t)]]]\stackrel{~}{\rho }^{(E)}\},`$
where

$$V_m^{(S)}(s)e^{L_S^{(S)}(s)}V_m^{(S)},$$
(84a)
$$q_n^{(E)}(s)\mathrm{\Lambda }_E^{(E)}(s)q_n^{(E)}.$$
(84b)
In terms of the correlation functions specified by Eq. (43) we can write

$`K[P_0\rho (t)]=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle _0^{\mathrm{}}}ds{\displaystyle \underset{n}{}}\{`$ (85)
$`q_n(s)q_nV_nV_n(s)P_0\rho (t)`$
$``$ $`q_nq_n(s)V_nP_0\rho (t)V_n(s)`$
$``$ $`q_n(s)q_nV_n(s)P_0\rho (t)V_n`$
$`+`$ $`q_nq_n(s)P_0\rho (t)V_n(s)V_n\},`$
where $`V_nV_n^{(S)}𝕀^{(E)}`$ and $`𝕀^{(E)}`$ is simply the identify operator on the environment subspace. We note that here $`V_n(s)=e^{L_Ss}V_n`$ is not the Heisenberg evolved operator, that is $`V_n(s)e^{L_S^{}s}V_n`$.

While Eq. (83) is similar to previous results, the main difference comes from the nature of $`L_E^{(E)}(t)`$. Previous authors have taken $`L_E^{(E)}(t)`$ to necessarily generate unitary evolution in order to use properties analogous to Eq. (78) which does not generally apply to nonunitary evolution.

Since Eq. (83) has the same mathematical structure as results obtained without dissipative effects in the environment, we expect similar shortcomings. In particular, to insure complete positivity for the reduced dynamics, we will now apply an averaging process (sometimes refereed to as the Rotating Wave Approximation),Davies (1974); Alicki and Lendi (1987) defined by:

$$\overline{K}=\underset{a\mathrm{}}{lim}\frac{1}{a}_0^ae^{L_0\tau }Ke^{L_0\tau }𝑑\tau .$$
(86)
Since $`L_EP_0`$

$$e^{L_E\tau }P_0=P_0,$$
(87)
so that with Eq. (46) and Eq. (81) we can write

$$\overline{K}=\underset{a\mathrm{}}{lim}\frac{1}{a}_0^ae^{L_S\tau }Ke^{L_S\tau }𝑑\tau .$$
(88)
Using Eq. (78) and Eq. (82),Eq. (88) becomes

$`\overline{K}[P_0\rho (t)]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}\underset{a\mathrm{}}{lim}{\displaystyle \frac{1}{a}}{\displaystyle _0^a}d\tau {\displaystyle _0^{\mathrm{}}}ds{\displaystyle \underset{n}{}}\{`$ (89)
$`q_n(s)q_nV_n(\tau )V_n(s\tau )P_0\rho (t)`$
$``$ $`q_nq_n(s)V_n(\tau )P_0\rho (t)V_n(s\tau )`$
$``$ $`q_n(s)q_nV_n(s\tau )P_0\rho (t)V_n(\tau )`$
$`+`$ $`q_nq_n(s)P_0\rho (t)V_n(s\tau )V_n(\tau )\}.`$
The operators $`\{V_n\}`$ can be decomposed in terms of the energy eigenstates of the system hamiltonian $`H_{(S)}`$:

$`V_n`$ $`=`$ $`{\displaystyle \underset{\mu \nu }{}}|\mu \mu |V_n|\nu \nu |`$ (90)
$`=`$ $`{\displaystyle \underset{\mathrm{\Delta }\omega }{}}{\displaystyle \underset{ϵ_\mu ϵ_\nu =\mathrm{}\mathrm{\Delta }\omega }{}}|\mu \mu |V_n|\nu \nu |`$
$`=`$ $`{\displaystyle \underset{\mathrm{\Delta }\omega }{}}V_{n,\mathrm{\Delta }\omega }.`$
It is clear that with this decomposition

$$V_{n,\mathrm{\Delta }\omega }=V_{n,\mathrm{\Delta }\omega }^{}.$$
(91)
The time dependence of the operators $`\{V_n(s)\}`$ becomes

$$V_n(s)=\underset{\mathrm{\Delta }\omega }{}e^{i\mathrm{\Delta }\omega s}V_{n,\mathrm{\Delta }\omega }.$$
(92)
When Eq. (92) is substituted into Eq. (89), there will be oscillating terms which the integral over $`\tau `$ will cancel, via

$$\underset{a\mathrm{}}{lim}\frac{1}{a}_0^a𝑑\tau e^{i(\mathrm{\Delta }\omega ^{}+\mathrm{\Delta }\omega )}=\delta _{\mathrm{\Delta }\omega ^{},\mathrm{\Delta }\omega }$$
(93)
so that Eq. (89) becomes

$`\overline{K}[P_0\rho (t)]`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle _0^{\mathrm{}}}ds{\displaystyle \underset{n,\mathrm{\Delta }\omega }{}}e^{i\mathrm{\Delta }\omega s}\{`$ (94)
$`q_n(s)q_nV_{n,\mathrm{\Delta }\omega }V_{n,\mathrm{\Delta }\omega }^{}P_0\rho (t)`$
$`q_nq_n(s)V_{n,\mathrm{\Delta }\omega }P_0\rho (t)V_{n,\mathrm{\Delta }\omega }^{}`$
$`q_n(s)q_nV_{n,\mathrm{\Delta }\omega }^{}P_0\rho (t)V_{n,\mathrm{\Delta }\omega }`$
$`+q_nq_n(s)P_0\rho (t)V_{n,\mathrm{\Delta }\omega }^{}V_{n,\mathrm{\Delta }\omega }\}.`$
With the definition

$$h_{n,\mathrm{\Delta }\omega }+iS_{n,\mathrm{\Delta }\omega }_0^{\mathrm{}}𝑑se^{i\mathrm{\Delta }\omega s}q_n(s)q_n,$$
(95)
and using Eq. (91) we can write Eq. (94) as

$`\overline{K}[P_0\rho (t)]`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[{\displaystyle \underset{n,\mathrm{\Delta }\omega }{}}{\displaystyle \frac{1}{\mathrm{}}}S_{n,\mathrm{\Delta }\omega }V_{n,\mathrm{\Delta }\omega }V_{n,\mathrm{\Delta }\omega }^{},P_0\rho (t)]`$ (96)
$`+{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle \underset{n,\mathrm{\Delta }\omega }{}}h_{n,\mathrm{\Delta }\omega }([V_{n,\mathrm{\Delta }\omega }^{}P_0\rho (t),V_{n,\mathrm{\Delta }\omega }]`$
$`+[V_{n,\mathrm{\Delta }\omega }^{},P_0\rho (t)V_{n,\mathrm{\Delta }\omega }]).`$
The first term in the right hand side of Eq. (96) is simply an additional hamiltonian term. The remaining terms (those with $`h_{n,\mathrm{\Delta }\omega }`$) are in the Lindblad form *if* $`h_{n,\mathrm{\Delta }\omega }`$ is positive. Using Eq. (43) and the moments given in Eq. (38) and Eq. (39) in Eq. (95) we find

$$h_{n,\mathrm{\Delta }\omega }=\frac{[(\lambda _n+\mu _n)^2+\mathrm{\Delta }\omega ^2]m_n^2D_{qqn}+D_{ppn}+m_n\lambda _n\mathrm{}\mathrm{\Delta }\omega +2(\lambda _n+\mu _n)m_nD_{pqn}}{[\lambda _n^2+(\mathrm{\Omega }_n+\mathrm{\Delta }\omega )^2][\lambda _n^2+(\mathrm{\Omega }_n\mathrm{\Delta }\omega )^2]}.$$
(97)
and

$$S_{n,\mathrm{\Delta }\omega }=\frac{C_{0n}+C_{1n}\mathrm{\Delta }\omega +C_{2n}\mathrm{\Delta }\omega ^2+C_{3n}\mathrm{\Delta }\omega ^3}{2\lambda m^2[\lambda _n^2+(\mathrm{\Omega }_n+\mathrm{\Delta }\omega )^2][\lambda _n^2+(\mathrm{\Omega }_n\mathrm{\Delta }\omega )^2]}.$$
(98)
where

$$C_{0n}=\mathrm{}m_n\lambda _n(\lambda _n^2+\mathrm{\Omega }_n^2)^2,$$
(99a)
$$C_{1n}=[(\mathrm{\Omega }_n^23\lambda _n^2)(\mu +\lambda )^2+(\lambda _n^2+\mathrm{\Omega }_n^2)^2]m_n^2D_{qqn}+(\mathrm{\Omega }_n^23\lambda _n^2)(2(\mu _n+\lambda _n)m_nD_{pqn}+D_{ppn}),$$
(99b)
$$C_{2n}=\mathrm{}\lambda m_n(\lambda _n^2+\mathrm{\Omega }_n^2),$$
(99c)
$$C_{3n}=m_n^2D_{qqn}[\lambda _n^2+\mathrm{\Omega }_n^2+(\mu _n+\lambda _n)^2]2D_{pqn}m_n(\mu _n+\lambda _n)D_{ppn}.,$$
(99d)
The denominator in the right hand side of Eq. (97) is the product of two sums of squares, and hence is guaranteed to be positive. The numerator is quadratic in $`\mathrm{\Delta }\omega `$ of the form

$`y`$ $`=`$ $`a\mathrm{\Delta }\omega ^2+b\mathrm{\Delta }\omega +c`$ (100)
$`=`$ $`(m_n^2D_{qqn})\mathrm{\Delta }\omega ^2m_n\lambda _n\mathrm{}\mathrm{\Delta }\omega `$
$`+m_n^2(\lambda _n+\mu _n)^2D_{qqn}+D_{ppn}+`$
$`2(\lambda _n+\mu _n)m_nD_{pqn}.`$
If $`y`$ is positive, then $`h_{n,\mathrm{\Delta }\omega }`$ is positive as well. The coefficient of the quadratic term is positive, thus $`y`$ has positive concavity. If $`y=0`$ has no real roots, then $`y`$ must be positive, which can be tested by the condition $`4acb^2>0`$. Upon some rearrangement, we can write

$`4acb^2`$ $`=`$ $`4m_n^2\{[m_nD_{qqn}(\lambda _n+\mu _n)+D_{pqn}]^2`$ (101)
$`+D_{qqn}D_{ppn}D_{pqn}^2({\displaystyle \frac{\lambda _n\mathrm{}}{2}})^2\}.`$
The first term inside the braces is a square and hence positive, while the remaining terms satisfy Eq. (20) and so the resulting expression is necessarily positive. Thus it is sufficient that our open system model for the oscillators is of the Lindblad form (as discussed in the previous section) to guarantee that the rotating wave approximation of the weak coupling limit generates a Lindblad form for the evolution of the system.

## V  Example: Oscillator linearly coupled to bath

In this section we illustrate our results with a test model consisting of an oscillator linearly coupled to a damped oscillator bath. To simplify notation we take $`\rho `$ to be the reduced density operator for the system. The system has a nominal hamiltonian given by

$$H_S=\frac{1}{2}m_S\omega _S^2Q^2+\frac{P^2}{2m_S},$$
(102)
and the interaction with the environment is given by

$$U_I=\underset{n}{}C_nQq_n.$$
(103)
With this choice, we can write the operators $`V_{n,\mathrm{\Delta }\omega }`$ in terms of the system creation and annihilation operators $`a^{}`$ and $`a`$

$$V_{n,\mathrm{\Delta }\omega }=C_n\sqrt{\frac{\mathrm{}}{2m_S\omega _S}}(a^{}\delta _{\mathrm{\Delta }\omega ,\omega _S}+a\delta _{\mathrm{\Delta }\omega ,\omega _S}).$$
(104)
The contribution to the hamiltonian through $`\overline{K}`$, as it appears in Eq. (96), becomes

$`\mathrm{\Delta }H`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{n,\mathrm{\Delta }\omega }{}}(S_{n,\mathrm{\Delta }\omega }V_{n,\mathrm{\Delta }\omega }V_{n,\mathrm{\Delta }\omega })^{}`$ (105)
$`=`$ $`\mathrm{}\delta \omega _S(aa^{}{\displaystyle \frac{1}{2}})+\mathrm{\Delta }E,`$
where

$$\delta \omega _S=\underset{n}{}\frac{C_n^2}{2m_S\omega _S\mathrm{}}(S_{n,\omega _S}+S_{n,\omega _S})$$
(106)
and

$$\delta E=\underset{n}{}\frac{C_n^2}{2m_S\omega _S}(S_{n,\omega _S}S_{n,\omega _S}).$$
(107)
Thus $`\delta \omega _S`$ is simply a frequency shift for the system and $`\delta E`$ is a C-number shift in the energy. The remaining contributions to $`\overline{K}`$ essentially are of the form of Eq. (19) with $`\mu =0`$, that is

$`\overline{K}[\rho ]`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\mathrm{}\delta \omega _S(aa^{}{\displaystyle \frac{1}{2}})+\mathrm{\Delta }E,\rho ]`$ (108)
$`{\displaystyle \frac{i}{2\mathrm{}}}\lambda ([q,\{p,\rho \}][p,\{q,\rho \}])`$
$`{\displaystyle \frac{D_{pp}}{\mathrm{}^2}}[q,[q,\rho ]]{\displaystyle \frac{D_{qq}}{\mathrm{}^2}}[p,[p,\rho ]].`$
The dissipation and diffusion coefficients are given by

$`\lambda `$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{C_n^2}{2\mathrm{}m_S\omega _S}}(h_{n,\omega _S}h_{n,\omega _S})`$
$`Dpp`$ $`=`$ $`(m_S\omega _S)^2Dqq`$ (109)
$`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{C_n^2}{4}}(h_{n,\omega _S}+h_{n,\omega _S})`$
If the environment degrees of freedom are thermalized (i.e. driven asymptotically to a Gibbs state) then Eq. (31) holds and the dissipation and diffusion coefficients are given by

$`\lambda `$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{C_n^2\lambda _n}{m_S^2(\lambda _n^2+(\mathrm{\Omega }_n\omega _S)^2)(\lambda _n^2+(\mathrm{\Omega }_n+\omega _S)^2)}}`$
$`D_{pp}`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{C_n^2\mathrm{}\mathrm{coth}(\frac{\mathrm{}\omega _n}{2k_BT})[(\lambda _n+\mu _n)(\lambda _n^2+\mathrm{\Omega }_n^2)+(\lambda _n\mu _n)\omega _S^2]}{m_S\omega _n(\lambda _n^2+(\mathrm{\Omega }_n\omega _S)^2)(\lambda _n^2+(\mathrm{\Omega }_n+\omega _S)^2)}}`$ (110)
and the frequency shift is given by

$$\delta \omega _S=\underset{n}{}\frac{C_n^2(\lambda _n^2+\mathrm{\Omega }_n^2\omega _S^2)}{2m_S^2\omega _S(\lambda _n^2+(\mathrm{\Omega }_n\omega _S)^2)(\lambda _n^2+(\mathrm{\Omega }_n+\omega _S)^2)}$$
(111)
From Eq. (33), the dissipative time scales of the environment are determined by $`\lambda _n\mu _n`$ and $`\lambda _n+\mu _n`$. If dissipation is weak, then

$$\lambda _n\mu _n,\lambda _n\mu _n\omega _n$$
(112)
and

$$\mathrm{\Omega }_n\omega _n.$$
(113)
The dissipation and diffusion coefficients each have a factor of

$`{\displaystyle \frac{1}{(\lambda _n^2+(\mathrm{\Omega }_n\omega _S)^2)}}`$ (114)
$``$ $`{\displaystyle \frac{1}{(\omega _n\omega _S)^2}}`$
from which we see that for a weakly damped environment, the greatest effect is from the oscillators in the environment with frequencies close to the system’s frequency. With this approximation, we can further simplify the dissipation and diffusion coefficients with:

$`\lambda `$ $``$ $`{\displaystyle \underset{\omega _n\omega _S}{}}{\displaystyle \frac{C_n^2}{4m_S^2\omega _S^2\lambda _n}}`$
$`Dpp`$ $``$ $`{\displaystyle \frac{\mathrm{}m_S\omega _S}{2}}\mathrm{coth}{\displaystyle \frac{\mathrm{}\omega _S}{2k_BT}}{\displaystyle \underset{\omega _n\omega _S}{}}{\displaystyle \frac{C_n^2}{4m_S^2\omega _S^2\lambda _n}}.`$ (115)
Thus the parameters of $`\overline{K}`$ satisfy Eq. (31) and the system is driven towards thermal equilibrium by an effective evolution of the Lindblad form.

## VI  Conclusions

We have constructed open environment models to account for the environment’s interaction with the “rest of the universe”. In Section II, our brief investigation of a classical model provides some insight and expectations for the development of a quantum mechanical model. The review the open system model in Section III served to establish many properties used in our derivation of the effective master equation in Section IV. Although Section IV largely follows previous work, the introduction of nonunitary evolution for the environment provided some novel aspects to the master equation derivation. We were able to show that an effective master equation of the Lindblad form could be obtained for a rich family of dissipative environment models. Our illustration of the resulting master equation with a bilinear environment-system interaction provided an demonstration in which an environment dynamically driven toward thermal equilibrium can naturally result in the dynamical thermalization of the system of interest.


text:
# 1 Introduction

## 1 Introduction

The rheological properties of soft materials, such as colloidal dispersions, presumably originate in a number of physical mechanisms, like shear–induced phase transitions, direct potential and hydrodynamic interactions, advection of fluctuations, and shear banding or localization among others; see e. g. the collection of papers in . At higher particle concentrations, the non–linear rheology depends on how steady shearing interferes with solidification during glass formation. Recently, we developed a theory for the non-linear rheology of dense colloidal suspensions aimed at this point . It describes how the structural dynamics is fluidized by advection of density fluctuations, while hydrodynamic interactions, non–linear flow profiles and ordering phenomena are neglected. Computer simulation studies can ensure that the latter processes are absent and thus provide crucial tests of the presented scenario. In this contribution, theoretical calculations will be compared to Brownian dynamics simulations of hard spheres by Strating  — without adjustable parameter in principle —, and to molecular dynamics simulations of a sheared binary Lennard–Jones mixture by Berthier and Barrat .

## 2 Theory

### 2.1 General aspects

A system of Brownian particles is studied in a prescribed steady shear solvent flow with constant velocity gradient and shear rate $`\dot{\gamma }`$. The equation of motion for the temporal evolution of the many–particle distribution function is known , and has been solved for hard spherical particles at low densities . This model constitutes a first microscopic approach to real dense colloidal suspensions, and may serve as a model sheared glassy fluid . It considers the “Brownian part” of the viscosity only, which, in Stokesian Dynamics simulations, Foss and Brady found to dominate compared to the hydrodynamic one for small shear rates

While the (approximate) approach developed in  gives general steady state quantities (like the shear distorted static structure factor) and their time–dependent fluctuations close to glassy arrest, we will concentrate on the thermodynamic shear stress $`\sigma (\dot{\gamma })`$ and the connected shear viscosity $`\eta (\dot{\gamma })=\sigma /\dot{\gamma }+\eta _{\mathrm{}}`$; here $`\eta _{\mathrm{}}`$ is the viscosity of the background solvent. The equations of motion exhibit a glass transition bifurcation, around which asymptotic expansions capture the transition from shear–thinning fluid flow to solid–like yielding. With the separation parameter $`\epsilon `$ denoting the (relative) distance from the transition, and $`t_0`$ a time scale obtained by matching onto microscopic transient motion, the following behaviors of $`\sigma `$ in the “structural window” have been established

$$\sigma =\sigma (\dot{\gamma }t_0,\epsilon )\{\begin{array}{cc}\dot{\gamma }t_0|\epsilon |^\gamma c_1\hfill & \epsilon <0\hfill \\ c_2\left(1+c_3|\dot{\gamma }t_0|^m\right)\hfill & |\epsilon ||\dot{\gamma }t_0|^{\frac{2a}{1+a}}\hfill \\ c_2(1+c_4\sqrt{\epsilon })\hfill & \epsilon >|\dot{\gamma }t_0|^{\frac{2a}{1+a}}\hfill \end{array},$$
(1)
where the $`c_i`$ are positive material–dependent parameters (for the exponents $`\gamma `$ & $`a`$ see e.g. ). The “structural window”, here, is defined as the double regime $`|\epsilon |1`$ and $`|\dot{\gamma }t_0|1`$, where the slowing–down of the structural dynamics dominates the steady state stress. While the divergence of the Newtonian viscosity $`\eta _0=t_0|\epsilon |^\gamma c_1`$ (first line of Eq. (1)) upon approaching the transition, applies to the linear–response regime of a fluid ($`\epsilon <0`$), and is known from mode coupling theory (see the references in ), the novel predictions close to and above ($`\epsilon 0`$) the transition describe the universal non–linear response of glasses to steady shearing with rate $`\dot{\gamma }`$. Importantly, a “dynamic yield stress” $`\sigma ^+(\epsilon )=\sigma (\dot{\gamma }0+,\epsilon 0)`$ is obtained because a finite stress has to be overcome in order to force the glass to yield even for vanishingly small shear rate; $`\sigma ^+`$ is connected to the constants $`c_2`$ and $`c_4`$ in Eq. (1). While the yield stress varies strongly with distance to the transition deep in the glass, at fixed parameters close to the transition, the stress increases from $`\sigma ^+`$ with a power–law in $`\dot{\gamma }`$, where the material–dependent exponent $`m`$ lies around 0.15 in the models studied below. The given asymptotes are only the leading orders for $`\epsilon 0`$ and $`\dot{\gamma }t_00`$, while corrections can be obtained systematically , or are included in model calculations to be presented below.

The dominance of the structural dynamics in determining $`\sigma (\dot{\gamma }t_0,\epsilon )`$ entails that all exponents or constants are functions of the equilibrium structure factor $`S_q`$ alone, except for the time scale $`t_0`$ which matches to shorter non–structural dynamics. Thus, hydrodynamic interactions or inertial terms only influence the value of $`t_0`$, which ideally could be determined from an analysis of the intermediate scattering functions of the system . This result arises because the small–shear rate rheology of glassy suspensions is dominated by steric hindrance (the “cage–effect”) which is not qualitatively affected by the properties of the solvent around the particles. It is in agreement with the findings in Stokesian dynamics simulations  that shear thinning is dominated by a decrease of the Brownian part of the viscosity. The elimination of particle forces in favour of the quiescent–state structure factor $`S_q`$ is an approximation of unknown quality in the present situation, but in part motivated by the consideration of small shear rates.

### 2.2 Models and simplifications

The equations of motion, from which $`\sigma `$ in Eq. (1) follows uniquely for a given $`S_q`$, have not been solved yet. Two approximate models were presented and discussed in  and shall be used in the following. While the schematic $`F_{12}^{(\dot{\gamma })}`$–model only incorporates the competition of two effects (divergent structural relaxation times with increasing $`\epsilon `$ and loss of memory induced by shearing), the semi–microscopic ISHSM combines a semi–quantitative description of a quiescent hard sphere colloidal dispersion  with an isotropically–averaged shear advection of density fluctuations. Both models only depend on two parameters which map onto $`\epsilon `$ and $`\dot{\gamma }`$ introduced in Eq. (1), and thus can be viewed as minimal models for the described scenario.

A problem when analysing data using both models arises from the ratio $`c_2/c_1`$ in Eq. (1), which has a simple physical meaning. It gives the ratio of yield stress to transverse elastic constant (viz. shear modulus $`G_{\mathrm{}}`$) of the glass at the critical point, $`c_2/c_1=\widehat{c}_1\sigma _c^+/G_{\mathrm{}}^c`$, where the numerical constant $`\widehat{c}_1=`$ 1.0 (1) for the ISHSM ($`F_{12}^{(\dot{\gamma })}`$). This ratio can be interpreted as a critical yield strain. Both models underestimate the effect of shearing leading to $`\sigma _c^+/G_{\mathrm{}}^c=`$ 0.33 (0.34) for the ISHSM ($`F_{12}^{(\dot{\gamma })}`$), while experiments give values around 0.05 indicating smaller strains are necessary for yielding . While the schematic $`F_{12}^{(\dot{\gamma })}`$–model is not meant to quantitatively capture such ratios, this error in the ISHSM presumably arises from the oversimplified handling of the shear–induced anisotropy of density fluctuations. The ISHSM treats all directions equivalent to the vorticity direction that is perpendicular to the flow plane. Perhaps unsurprisingly this underestimates the effects of shearing. We correct for this error in an ad–hoc fashion by rescaling the shear–rate $`\dot{\gamma }`$ when considering $`\eta (\dot{\gamma })`$. For the $`F_{12}^{(\dot{\gamma })}`$–model this procedure is rigorously equivalent to an adjustment of the ratio $`\sigma _c^+/G_{\mathrm{}}^c`$.

## 3 Results and comparison with simulation data

Before applying Eq. (1) to a solution of colloidal hard spheres at packing fraction $`\varphi `$, at first the latter’s critical value, $`\varphi _c`$, entering in $`\epsilon =C(\varphi _c\varphi )/\varphi _c`$ (with $`C=1.5`$ ) needs to be determined. This is done by testing whether the divergence of the quiescent viscosity (and corresponding structural relaxation time) for $`\epsilon 0`$ is observed. The inset of fig. 1 shows viscosities from experiments  and from Brownian dynamics simulations . Also included are self–diffusion coefficients from , which are predicted to vanish with $`D|\epsilon |^\gamma `$. Replotting the data with the calculated $`\gamma =2.62`$ , fits to the data above $`\varphi 0.50`$ give $`\varphi _c=`$ 0.57 , 0.58 & 0.60 for $`D`$, $`\eta `$ (two outliers neglected) and simulations, respectively. Interestingly, the two experimental data sets provide rather close estimates for $`\varphi _c`$ and indicate a strong coupling of diffusion and viscosity, $`D\eta 0.4\eta _{\mathrm{}}D_0`$ for $`\epsilon 0`$ (neglecting the difference in $`\varphi _c`$), with $`D_0`$ the dilute single particle diffusion coefficient. The numerical factors is about half the predicted value . We speculate that the discrepancy of the extrapolation of the simulation results arises in part because the data are not fully in the asymptotic regime.

With the quantitative knowledge of $`\epsilon `$, only the matching time $`t_0`$ is required to analyse the steady state viscosities in the structural window using the ISHSM. We chose to obtain it via the full fitting procedure which consists in matching by eye the numerical solutions to the non–Newtonian viscosity data. In this way, $`t_0`$ is mainly determined by the increase of the Newtonian viscosity, because $`\eta t_0\sigma /(\dot{\gamma }t_0)`$ holds and $`\sigma /(\dot{\gamma }t_0)`$ becomes independent of $`t_0`$ in the fluid for vanishing shear rate. The main panel of fig. 1 shows $`\eta `$ from the Brownian dynamics simulations as function of the dimensionless Peclet number Pe$`{}_{0}{}^{}=\dot{\gamma }d^2/D_0`$, which measures the effect of shearing relative to the time a single particle diffuses its diameter $`d`$. The fits by eye using the ISHSM are included for packing fractions close to the transition, $`\varphi 0.50`$. From the fits, mainly from the divergence of $`\eta _0`$ given in Eq. (1), the matching time $`t_0=0.019d^3\eta _{\mathrm{}}/k_BT`$ is estimated, and inclusion of corrections to asymptotic behavior in the ISHSM–fits shifts the glass transition packing fraction closer to the other determinations; $`\varphi _c=0.59`$ follows from the $`\epsilon `$ used in fig. 1. Note that the solvent viscosity is included in the theoretical curves, $`\eta =\eta _{\mathrm{}}+t_0\sigma /(\dot{\gamma }t_0)`$111The ISHSM calculations provide $`\sigma (\epsilon ,\dot{\gamma }t_0,`$Pe$`{}_{0}{}^{})`$ for all values of $`\epsilon `$ and Pe0, while Eq. (1) captures the asymptotic behavior for $`\epsilon 0`$, $`\dot{\gamma }t_00`$ & Pe$`{}_{0}{}^{}0`$. Because we aim at describing the proximity of the glass transition, we match the parameters of Eq. (1) ($`\varphi _c`$ and $`t_0`$) there. Without matching the ISHSM gives $`\varphi _c=0.52`$ & $`t_0=0.025d^3\eta _{\mathrm{}}/k_BT`$. . In the shear–thinning region, the viscosity diminuishes and approaches a behavior like $`\eta \sigma _c^+/\dot{\gamma }`$ with strong corrections, though, masking the power–law . Because of the overestimate of $`\sigma _c^+`$ in the ISHSM, this decrease would set in at too high $`\dot{\gamma }`$ values only. In order to correct for the quantitative error, the theoretical curves are plotted versus rescaled shear rate, $`\dot{\gamma }0.25`$; i.e. $`\sigma /(\dot{\gamma }t_0)=f_\eta (\dot{\gamma }0.25)`$. With this ad–hoc correction, satisfactory agreement of theory and simulation results is seen for Pe$`{}_{0}{}^{}1`$, where the steady state viscosity varies over two orders on variation of shear rate and packing fraction. For larger Peclet numbers, the data presumably lies outside the structural window where Eq. (1) applies. Motivated by numerical findings in , we speculate that the enhanced, $`\varphi `$–dependent steady state viscosities around Pe$`{}_{0}{}^{}=`$ 1 – 10 in the Brownian dynamics simulations originate in the hard core repulsion. If so, hydrodynamic interactions which prevent particles from close contact could appreciably affect $`\eta `$ in this region.

A second set of steady state shear stresses and viscosities is provided by recent large scale molecular dynamics simulations of a sheared simple liquid (a binary Lennard–Jones mixture) above and below its glass transition temperature . Because kinetic parameters do not enter the theoretical predictions, and as linear flow profiles were obtained in the simulations, the universal predictions of our approach can be compared again. Figure 2 shows stationary shear stresses from the simulations and fits by eye using the $`F_{12}^{(\dot{\gamma })}`$–model as specified in . The model provides a relaxation time $`\tau `$ as function of $`\epsilon `$ and a dimensionless shear rate, denoted as Peclet number Pe0, which are mapped onto the data as specified in the figure caption. The data nicely span the glass transition temperature, $`T_c0.435`$ already known , and are well compatible with a transition from a shear–thinning fluid to a yielding glass with finite yield stresses at and below the transition temperature.

## 4 Conclusions and outlook

We presented results of a microscopic theory of the nonlinear rheology of colloidal fluids and glasses under steady shear , and compared them with simulation and experimental data. This brought out the existence of a universal transition between shear-thinning fluid flow, with diverging viscosity upon increasing the interactions, and solid yielding, with a yield stress that is finite at (and beyond) the glass point. Numerical calculations could explain simulation results over up to 3 decades in variation in shear rate and viscosity. A quantitative analysis of larger data sets is required in order to determine the theoretical parameters for both simulations more accurately than the estimates found here.

The approach we outlined should be improved with respect to the handling of shear–induced anisotropies, and stress–induced effects. The latter may lead to shear thickening behaviour that, for many colloidal materials, occurs at higher flow rates than those addressed here. This avenue will be explored in a future paper  on a version of the schematic model which is modified to include explicit stress- (as well as strain–rate–) dependence.

We thank J.–L. Barrat and L. Berthier for discussions. M.F. was supported by the DFG, grant Fu 309/3.