{'fraction_non_alphanumeric': 0.05009367530202209, 'fraction_numerical': 0.014199883713418179, 'mean_word_length': 4.147517541817698, 'pattern_counts': {'":': 1, '<': 17, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum and statistical mechanics in open systems: theory and examples

Using the system-bath model Hamiltonian this thesis covers the equilibrium and out of equilibrium properties of quantum open systems. Topics included are the calculation of thermodynamical quantities of open systems, derivation of quantum master equations, phase space and numerical methods and Linear and non Linear Response Theory. Applications are the transport in periodic potentials and the dynamics of spins.


Chapter 1 introduction
In physics the "system" is the small part of the universe on which we have some control. In other words, it is the part of "interest" that we measure and/or manipulate. Knowing the system Hamiltonian, H s , we can calculate its time evolution and find the future state of the system, as well as telling its past (how it reached the current state). The system of interest can be, for instance, an electron, an atom, a piece of metal, a cup of coffee, or the solar system.
Everyday experience teaches us that any system interacts with its environment. For example, we all know that a coffee thermos in a cold winter morning will eventually cool down, no matter how sophisticated the flask is. Nobody is astonished by this: we cannot fully isolate the coffee from the outside. Admittedly, in many cases the effects of the environment can be neglected, and the dynamics is then determined by the system Hamiltonian alone. In other cases, however, we must take into account that our system is not isolated, but coupled to its surroundings [1,2].
Either because we do not know the evolution of the environment (on which we have not control), or because we only look at the system dynamics, we do not have complete info on the whole (system plus bath). Think about classical mechanics, for example; one point in the system's phase space may correspond to several configurations of the total system. If we cannot fully determine the current state, we are unable to reconstruct the past or predict the future evolution from H s alone. On the other hand, the total phase space of system plus environment is very large. As a result the Poincare recurrence time (to return to the initial state) becomes infinite for most practical purposes: we never find the coffee warming up after some time (irreversibility).
In this context the formalism of open systems is important in that we need not change the current theoretical frameworks to account for irreversibility. In addition, this formalism provides some clues for the difficulty in finding quantum superpositions in everyday macroscopic life. The theory saves face.
In the following pages we study the effects of the environment on the dynamics of the system of interest. Although we start surveying the classical limit, our work will focus on the quantum domain (non-relativistic). A constant throughout the manuscript will be to relate the results obtained with their classical limits (appropriately defined).

CHAPTER 1. INTRODUCTION
Although the formalism of the first chapters is general, the theory will be applied later on to systems with one degree of freedom. Examples from solid-state physics where the calculations can be applied are the phase difference across a Josephson junction, or collective excitations in lattices, like solitons or fluxons [1,3]. In these systems openness is important, while they can be approximately described by a few collective degrees of freedom. We will also tackle spin systems, like quantum superparamagnets or ensembles of 1/2-spins (two-level systems) [4,5].
As for the spirit of this manuscript, we have chosen not to repeat here the details of the calculations included in the articles we have written in the last three years. The reader interested in those details can consult our publications, which are cited as PAPER I, II,. . . In any case the thesis is quite self-contained, and can be followed without recourse to the articles (at least that was our aim).
The text is organized as follows. We begin reviewing the classical limit in Chap. 2. We introduce the bath-of-oscillators formalism for the environment, obtaining reduced dynamical equations (Langevin & Fokker-Planck equations). Then we quantize and proceed to derive dynamical equations for the reduced density matrix (quantum master equations; Chap. 3). Chapter 4 is devoted to the thermal equilibrium properties of open quantum systems. Then in Chapter 5 we introduce the phase-space formalism and apply it to quantum master equations, giving a natural link to the classical equations. This first part is general, and we close it by introducing in Chap. 6 the tools we will employ later on in explicit calculations: linear-response theory & the continued-fraction technique to solve master equations. The last two chapters (the longest) are devoted to the application of the general formalism to two problems of interest. In chapter 7 we study the Bloch electron (a particle in a periodic potential) coupled to a dissipative bath. Finally, in Chap. 8 we study the equilibrium and dynamical properties of quantum superparamagnets. 12 Chapter 2 classical theory: irreversibility and quantization In this chapter we place the problem in a historical context, beginning with the 19th century discovery of Brownian motion. Then we introduce the bath-of-oscillators formalism and address in this frame the problems of irreversibility and quantization. We conclude introducing the specific systems we will study in this thesis.

history of Brownian motion [6, Chap. 1]
In 1827 Father Brown observes under his microscope the irregular motions executed by pollen grains when suspended in a fluid (figure 2.1). Such a motion disagreed with physical laws known at the time, where particles followed smooth and predictable trajectories according to Newton laws (we left aside chaos, unknown then). Due to his background, he was a botanist, Brown thought that the animated motion observed might be a manifestation of life. Soon he would be discarding such an explanation. 13 It wasn't until 1905 when Einstein, in his annus mirabilis, explains in a satisfactory way the phenomenon, rightly called Brownian motion. 1 Einstein's approach is based on two main points: (i) the irregular motion is caused by the impact of the fluid's molecules on the colloidal grain (ii) the motion of those molecules is so intricate that the grain motion can only be described statistically. From these assumptions an equation for W (x, t) could be derived, giving the probability of the pollen grain being around x at time t (a sort of Liouville equation, but in configuration space).
This is a diffusion equation where D is related with the standard deviation of the grain position x(t) 2 − x(0) 2 = 2D t.
A few years later Langevin proposed an explanation, a model, which using his own words was "infinitely more simple" than that of Einstein & Smoluchowsky. He assumed that the force due to the fluid can be decomposed into: (i') a systematic part (viscous force) and (ii') a fluctuating force, accounting for the unpredictability of the collisions. Incorporating these forces into Newton's equation, Langevin wrote: Here γ is the viscous damping parameter and ξ(t) the fluctuating force, on which he only assumed that it can take random values both positive and negative. This was the first example of a stochastic differential equation [6,7]. Except for the short time regime (inertial), the treatments of Einstein and Langevin (more complete) lead to equivalent results. Both equation (2.1) and (2.2) describe the pollen grain dynamics. The Liouville and Newton time evolutions, respectively, get modified due to the interaction with the fluid, in such a way that the effective equations for the grain display irreversible dynamics.

problems with irreversibility and quantization
At this point we could raise the following question. Assuming that, though certainly complex, both the pollen grain and the fluid molecules obey Newton's reversible laws, how comes one can arrive at irreversible equations? This is the irreversibility paradox [2]: To this, a second complain was added: 2 1 In 1906 Smoluchowsky's paper is published, putting forward the same explanation of Brownian motion. We would like to dispel an urban legend attributing the name to the irregular and unpredictable motions on stage of the late James Brown. 2 These objections were first raised against Boltzmann's kinetic equation [8].
• Zermelo's paradox : If we start from a configuration with many pollen grains at the center of the fluid, we know that after some time they will spread throughout the volume, according to the diffusion equation (2.1). That is, we will never find that at some later time all grains return to the center, the initial configuration. This is at odds with Poincare's theorem, which asserts that the trajectories of a bounded system (in phase space) will pass arbitrarily close to the initial state, after a time called the recurrence time.
Finally, along with the irreversibility/recurrence issues, we find the problem of quantization. How could we quantize an equation that is irreversible, or stochastic? In general phenomenological quantization schemes (say, the counterpart of plugging −γẋ in Newton equation) can lead to violations of basic issues (like the very normalization of the state, or commutation relations). 3 With the formalism introduced in next section we will try to answer both questions.

bath of oscillators formalism
Here we discuss the rigorization of the diffusion and Langevin equations. This will allow us to address the issues of irreversibility and quantization.
As physicists we feel comfortable when we see a Hamiltonian. Therefore we formally write: This total Hamiltonian is made of the system Hamiltonian H s , the environment H b (the bath), and the interaction between them H sb . Naturally, we will have to give content to each part of H tot . This looks like a though task though. On top of the problem of modelling the system itself (think of nuclear physics), now we have to model everything in its surroundings. Fortunately, in many cases of interest one has the following simplifying conditions: (i) The environment "feels" only weakly the presence of the system. Loosely speaking, the bath is big, a macroscopic object (but not necessarily classical); the liquid carrier in Brown's set up. The fluid's properties change very little when a pollen grain, or a handful, are suspended on it. Naturally, the converse does not hold. The grain is indeed affected by the bath.
(ii) The environment dynamics is faster that the characteristic time scales of the system. For example, due to their microscopic nature the fluid molecules move around so fast that they look like random to us.
Under these conditions the bath can be modelized as a set of harmonic oscillators linearly coupled to the system [10,11]: Here (x, p) are the system canonical variables, not necessarily coordinates and velocities. Note that we can extract the interaction as ∼ Q α F (x, p). Thus, the coupling is linear in the bath coordinates Q α , whereas it can be arbitrary in the system variables F (x, p). 4 When extracting the coupling ∼ Q α F (x, p), one also gets a term ∝ F 2 (x, p) which only depends on the system variables. This counter-term cancels the renormalization of H s produced by a purely linear ∼ Q α F (generally speaking, an arbitrary coupling induces fluctuations, dissipation and renormalizations). So what we are doing is assuming that H s is the Hamiltonian to which we have access experimentally, where all relevant renormalizations have already been included (see [10] and what follows).
A further case where Hamiltonian (2.4) can be used is: (iii) The bath oscillators provides an "exact" description of the environment. Think about the phonons in a solid, or the normal modes of the electromagnetic field, the photons, or less evident cases like electron-hole excitations.
Well, once we have introduced our working Hamiltonian, we would like to obtain dynamical equations for the (sub)system degrees of freedom. To this end we consider a dynamical variable A(x, p) depending only on the system, and write Hamilton's equations for it: dA/dt = {A, H tot } with { , } the Poisson bracket and the total Hamiltonian (2.4). Writing similarly the equations for the bath variables, Q α and P α , we have [12] dA The first line includes the contribution of the system dynamics in the absence of interaction (first term), the part from the renormalization terms (middle), and the coupling to the bath coordinates (last). The second line is simply the equations of motion for the bath oscillators (as forced by the system).
Here we start to appreciate the advantages of model (2.4), as we know how to solve Eqs. (2.7). They are just scores of non-interacting and forced oscillators, with F (x, p) playing the role of a driving. This was made possible by the coupling being linear in the bath coordinates. Now, from undergraduate physics we know that the solution of a harmonic oscillator forced by F (τ ) reads: (2.8) 4 The most general Hamiltonian linear in the bath would be Htot = Hs(x, p) + 1 2 Fα(x, p) " 2 (2.5) where one accounts for a possible coupling with the oscillator momenta, and for different F and G for each mode. However, to illustrate the formalism is enough with (2.4). The general (2.5) brings little new, while it makes the calculations cumbersome. 16

BATH OF OSCILLATORS FORMALISM
where Q h α = Q α (t 0 ) cos[ω α (t − t 0 )] + P α (t 0 )/ω α sin[ω α (t − t 0 )] is the general solution of the homogeneous equation (free oscillators, without F ), while the last term is a particular solution of the inhomogeneous equation.
Plugging now the solution (2.8) in the system equation (2.6), and integrating it by parts, one finds (β ≡ 1/k B T ) We see that the renormalization term in (2.6) was cancelled by the term − α u 2 α 2ω 2 α {A, F 2 } arising from the integration by parts. This was indeed the reason behind including ∼ F 2 in the starting Hamiltonian (2.4). It ensures that the renormalization (if physical) was already included in H s , and not counted twice; then H s is the Hamiltonian (the levels) to which we have access in experiments. As for the other two terms in (2.9), we are going to see they describe fluctuations and relaxation with We have extracted β from C(t − τ ) so that it matches the correlator of f (t) [Eq. (2.11) below], in correspondence with the quantum definition (chapter 3). From the equations above we begin to grasp the structure of a Langevin equation like (2.2). We would have the "random force" f (t), which is just the free evolution of the bath oscillators (see Fig. 2.2), while the memory integral in Eq. (2.9) would correspond to the viscous damping. The later originates from the inhomogeneous solution in (2.8), and brings back the dynamics at previous times (a sort of back-reaction on the system of its previous action on the bath). As for the fluctuations, recall that we do not have full control on the bath, whereas f (t) includes the initial conditions Q α (t 0 ), P α (t 0 ) for each oscillator. It seems reasonable to draw them from an equilibrium state at the initial time. 5 Being harmonic oscillators, the distribution of positions and momenta would be Gaussian, so we just need to specify the first two moments of f (t) Here we see the link mentioned between the integral kernel and the force correlator. The latter can be written in terms of the bath spectral density: in the following way [13] (2.14) Now we are ready to give a couple of examples to give content and fix ideas.

Brownian particle
Let us consider the Hamiltonian of mechanics H s = p 2 /2m + V (x), and F (x, p) = −x (i.e., the coupling is linear in the system coordinate as well). Assume now that the correlator C(τ ) is a function strongly peaked around τ = 0, expressing that the system loses memory quickly (Markovian limit). Quantitatively C(τ ) = 2mk B T γ δ(τ ), corresponding to J(ω) = mγ ω, with γ a measure of the coupling strength. 6 Setting first A = x and then A = p in the equation of motion (2.9) one gets This equation reduces to the one proposed by Langevin when H s = p 2 /2m (that is, for the free particle). A stochastic equation with Gaussian f (t) and δ-correlation f (t)f (τ ) ∝ δ(t − τ ), can be converted into an equation for the probability distribution W (Fokker-Planck equation; see e.g. [2,7]). For the mechanical problem (2.15) one obtains the Klein-Kramers equation [14]: This Fokker-Planck equation for a particle in a potential is a generalization of the diffusion equation (2.1). 7 6 The spectral density (2.14) is customarily assumed of the form J(ω) ∼ ω α . If α = 1 the bath is called Ohmic (Markovian case) and the cases α ≷ 1 are called super/sub-Ohmic. The choice α = 1 if often done for simplicity; it is hard to solve non-local (integro-differential) equations! But if required by the physics of the problem, one cannot escape working with α = 1, as we will do later on. 7 Actually, Einstein Eq. (2.1) follows when neglecting inertial effects, the so called overdamped regime. And it was for a free particle; in a potential the overdamped equation for the distribution W (x, t) = classical Brownian spin/dipole Fortunately the equation of motion (2.9) is valid for any Poisson bracket, so we can handle problems more general than the Brownian particle. Given that in this thesis we will study simple spin Hamiltonians, we are going to discuss first the classical case H s = H s (m x , m y , m z ) with the angular-momentum Poisson algebra Here θ and φ are the polar and azimuthal angles on the sphere. The constant g is the gyromagnetic ratio related with the spin number by S = µ B /g, with µ B Bohr's magnetron.
In this case we set A = m i in (2.9) and we get (Ohmic bath) 8 This equation has the structure of the phenomenological equation proposed by Landau-Lifshitz in 1935. Here (∂/∂m) i ≡ ∂/∂m i is the gradient in m-space, while the noise and the damping include: We have used an Ohmic J(ω) = λ ω, with λ the coupling strength; the analogue of the damping γ in particle problems (footnote 6). We will use two different letters for particles and spins following the standard notation in the literature (besides, γ has dimensions of frequency 1/t whereas λ is dimensionless). As we did for the Brownian particle, we can also write a Fokker-Planck equation associated to the spin Langevin equation (2.20) where now (∂/∂m)· is the divergence operator in m-space. Let us illustrate this with a simple case: coupling F ∼ S i linear in the spin variables plus axially symmetric H s (θ), as in H s ∝ −K cos 2 θ − B z cos θ. Under these conditions ∂ φ B ef = 0 and we can restrict ourselves to axially symmetric solutions ∂ φ W = 0. Then, introducing z = cos θ and 2τ N = β/(gλ), equation (2.22) reduces to [15]: (2.17) 8 Note the useful relation: valid for any function of the spin variables mi. It follows from the Poisson bracket (2.18) and the definition of the vector product ×. This is nothing but Debye's equation for dipolar relaxation (possibly the first Fokker-Planck equation describing rotational Brownian motion).

revisiting the paradoxes of irreversibility
Now we are in a position permitting to address the two objections of section 2.2: • First, let us tackle the time-reversal objection due to Loschmidt. It is clear that the Hamilton equations (2.6) and (2.7) are reversible. However, the generalized Langevin equation (2.9) is completely equivalent to them. What happens is that one thinks of reversing the velocities of the system only. Therefore, we have a "practical" irreversibility: although the whole system is reversible, we cannot change the sign of the velocities of all degrees of freedom, including those of the bath (on which in practice we have little control).
• Second, Zermelo's paradox. The fact that both the effective equations (and the experiments) show that the pollen grains are never brought back together, is just a consequence of the total phase space being very large. Again, we cannot overlook the bath. Indeed one can estimate the recurrence time and find that it diverges as we let the number of bath oscillators tend to infinity [8,16]. It is only for finite "baths" that such time remains finite. But in practical cases the recurrence time can be several times the estimated age of the universe. For practical purposes we can be sure that the grains will never return to the initial configuration.

quantization
We have just seen how to handle the irreversibility problem (get along with it, would be more correct). Let us move on to the issue of quantizing a system with Langevin or Fokker-Planck dynamics. The idea is to consider the Hamiltonian (2.4) as the problem to be quantized. Then we would just need to transform the variables (x, p) & (Q α , P α ) into Hilbertian operators and proceed by the book. Actually, the main derivation of Sec. 2.3 can be followed analogously in the quantum case; "simply" replace Poisson brackets { , } by commutators [ , ], with A an operator in the Heisenberg picture. This is how one derives quantum Langevin equations [17]. However, these are not as useful/generic as their classical counterparts; now both x and the force f (t) are operators, and calculations to-the-end exist only in simple cases [18,19] On top of this, there are no simple recipes (as those in the classical case), to convert a Langevin equation into an equation for the distribution (Fokker-Planck like).
To bypass these problems one usually resorts to quantum master equations. Recall that the bath is not under control, so it is out of question any preparation of pure states (of system plus bath). Then one needs the density matrix formulation (more general) to compute averages A , our main interest. If A acts only on the Hilbert space of our system, we can introduce the reduced density matrix ̺. Then using the properties of partial tracing, we get A just as a trace over the system: Here s + b indicates the original trace over the whole system+bath using the total density matrix ̺ tot . Much as in the classical case the distribution W (x, p) involved the system variables only, now we have ̺ playing this role. In this sense the equations of motion for ̺ will be the analogue of the Fokker-Planck equations. Those equations (next chapter) generalize the Von Neumann equation for the closed evolution: This simply parallels the fact that the Hamiltonian part of the Fokker-Planck was the Liouville equation

systems studied in this thesis
We are going to deal with both particle problems and spin Hamiltonians (translational and rotational problems, respectively). In the final applications we will consider Hamiltonians with one degree of freedom (but which could describe effectively composite systems like a Josephson junction or a molecular magnet).

particle systems
We will consider the standard Hamiltonian with [x, p] = i . In particular we will address the problem of a quantum particle running in a periodic potential V (x) = V (x + L) plus a constant force (this is related with the Wannier-Stark problem [20]). We will study the dynamics in a cosine potential ( Fig. 2.3): and in ratchet potentials (lacking inversion symmetry) In figure 2.3 we have plotted the energy bands for these potentials (see PAPER I, Sec. 7.3). We have introduced a dimensionless parameter with S 0 some relevant/characteristic action; for instance S 0 = E 0 /ω 0 with E 0 the barrier height and ω 0 the oscillation frequency at the bottom of the wells. The parameter K tells us how classical (or quantum) the system is. Letting K → ∞ gives the classical limit; then the bands tend to a continuum (for K = 200 the levels of Fig. 2.3 are already very close).

spin systems
In parallel with mechanical problems, we will consider simple spin Hamiltonians ( Fig. 2.4) Here [S i , S k ] = iǫ ijk S k , with S 2 |m = S(S +1)|m and S z |m = |m . The above is an effective spin Hamiltonian [21]. The term −B z S z is the Zeeman coupling to the external field, while D is the anisotropy constant, of spin-orbit origin. When D = 0, Hamiltonian (2.29) reduces to the often studied isotropic spin. The cases D ≷ 0 correspond to anisotropic spins, with "easy axes" and "easy plane" anisotropy, respectively. Later, in chapter 8 we will focus on D ≥ 0 in (2.29), which is the minimal model for superparamagnets [22]. The classical limit will be taken letting S = (µ B /g)/ tend to infinity. Note the analogy between K for periodic potentials and S here. When taking the limit, it will be convenient to "normalize" the operators as H s = −DS 2 (S z /S) 2 − B z S(S z /S). Then we let S grow keeping constant DS 2 and B z S, the amount of anisotropy and Zeeman energies. Mathematically, more and more levels are introduced, towards a continuum, keeping constant the barriers' height.

dynamics of open quantum systems
In this chapter we address the dynamics of open quantum systems. We first discuss the notion of quantum decoherence, through a solvable example. To handle more complex systems one can resort to quantum master equations, the quantum analogue of the Fokker-Planck equations of the previous chapter. Here we introduce those equations and give some examples. We discuss the approximations involved in their derivation, as well as other invoked in solving master equations.

quantum decoherence
In the previous chapter we discussed how irreversibility could emerge from reversible equations. Now we would like to see how the classical world we are familiar with arises out of the microscopic quantum world it is made of. A proposed explanation for such emergence resorts to the coupling with the environment, in a irreversible process called decoherence (see [23] for an introduction to the trendy topic). It can be seen as a consequence of having the dynamics of very many degrees of freedom but only looking at some reduced dynamics. The notion can be introduced with a simple example, due to Zurek [24], which is fully solvable.
The question one aims to answer is: if everything around us, including ourselves, is made of atoms obeying the laws of quantum mechanics, how comes we do not find superpositions of states in the macroscopic world? The most famous/shocking of those superpositions is due to Schrödinger: Or others equally exotic like |Ψ ∝ |here +|there , i.e., the possibility of having delocalized objects. To address these points let's go through the promised simple example. 1 As usual, we will consider that the system interacts with many degrees of freedom, so that the total quantum Hamiltonian is H tot = H s + H sb + H b . For simplicity, we assume that both system and bath have discrete and non-degenerate spectra. In the basis of eigenstates of H s and H b one can write Here |e α are the eigenstates of the bath oscillators, and |s j those of the system. In this basis, the H sb in the bath-of-oscillators model (2.4) can be written as To proceed we assume that σ ii ′ αα ′ ≪ γ iα , meaning that decoherence times are much shorter than relaxation processes. Then which allows to solve analytically the problem. If σ ii ′ αα ′ ≪ γ iα we would have to resort to a quantum master equation (next section), and solve it numerically, most likely. Let us assume now that our system is in a superposition like (3.1) and let us study the dynamical evolution. Since we have to specify the state of the bath too, we assume for simplicity that it is a pure state (though everything can be repeated for the realistic case of a mixed bath). Then, the initial state reads: Now with the H sb considered, the time evolution at time t results in the entangled state: 2 Calculating now the reduced density matrix (2.24), that is, ̺ = Tr b (|Ψ(t) Ψ(t)|) one finds in the {|s j } basis
Here a i a * j e − i (δ i −δ j )t is just the coherent rotation due to the evolution under H s , while z ij (t) is due to the interaction with the environment [ω ij α = (γ iα − γ jα )/ ]. In general ω ij α = 0 and hence z ij (t) results to be the sum of many oscillating factors with different frequencies. Then at long enough times z ij (t) → 0, so that e −iω ij α t can be seen as a random phase. Therefore |b α | 2 e −iω ij α t becomes a random vector in the complex plane. The problem can thus be mapped onto the random walk in 2D, where each move is |b α | 2 e −iω ij α t and |b α | 2 = 1 is kept fixed (see Fig. 3.1). Then z ij (t) is a random variable with a probability distribution given by P (z) ∝ e −N z 2 , where N is the number of bath oscillators. Therefore at sufficiently long times, and considering the bath to be macroscopic (N → ∞), one writes z ij (t) ∼ = 0 . (3.9) Or in other words, the off-diagonal elements drop to zero and the original superposition j α j |s j is "destroyed". In the sense that the reduced density matrix is indistinguishable from a statistical mixure with |a i | 2 as diagonal elements.
The time t d , from which we can take z ij (t) ∼ = 0, is called the decoherence time. The loose of coherence is a consequence of looking only at the system. In fact, everything follows from the unitary evolution of the total quantum system; that is, the coherence is still out there, but now is diffused/shared between a large number of degrees of freedom. In a sense, this decoherence follows from the "non-locality" of quantum mechanics, as the state (3.6) is fully entangled between system and bath.
Finally, is worth to recall that every density matrix is always diagonalizable (as it is Hermitian). Or conversely, given a diagonal density matrix, it can be seen as non-diagonal by some other basis. Then, the concept of decoherence seems to be "basis dependent".
What happens, decoherentists say, is that the interaction with the bath "determines" in which basis ̺ will be diagonal; the so called pointer basis [23]. It seems that Nature chose us to be localized objects, as well as being alive or death. To give some numbers, the decoherence time for a superposition of dust grains with a diameter 10 −3 cm, set 1 cm apart is estimated to be around 10 −36 seconds [25]. That is, for most practical purposes, and for most macroscopic objects, quantum superpositions cannot be observed. 3

derivation of master equations
It will not be always possible to express ̺(t) analytically, as in the previous section. Our goal here is to derive differential equations for ̺. In contrast to the classical case [Eq. (2.9)], we do not know how to obtain in general an exact equation for the model (2.4). Therefore, we will be happy with an approximate equation. The main assumption involved will be to consider that the system-bath interaction is weak.
The most usual techniques in deriving master equations are those of projection operators [2,28], path integral [9,29,30,31], or cumulant expansions [7,32]. For our purposes, it will suffice with a simpler derivation, based on plain quantum-mechanical perturbation theory. The approximations invoked will be clearly seen as we proceed.
We will start working with ̺ tot ; we would trace later on to obtain the effective dynamics for ̺ [Eq. (2.24)]. The evolution of ̺ tot is given by where we have assumed that H tot is not time dependent. This is not essential though; we could consider that H tot depends on t, and then the evolution operator would be a time-ordered exponential U (t; In what follows it will be convenient to split V = H sb from H = H s + H b , and assume that the system-bath interaction is weak. That is, we will treat V = H sb as a perturbation (we shall quantify later what we mean by "weak"; be patient!). We want to calculate U (t; t 0 ) perturbatively in V = H sb ; to that end, we will use the following Kubo-type identity [33, p. 148]: This equation 4 can be iterated by expressing e −(i/ )(H+V )s in the integrand in just the same way, and so on. To second order in V one has: The decoherence approach does not exclude the existence of those superpositions. It only asserts that they vanish at too short a time to be observable [26]. For a discussion of the state of the play in the search of quantum superpositions in macroscopic objects, see Leggett's insightful review [27]. On the other hand, decoherence does not seem to solve the "measurement problem". It is true that it provides a diagonal ̺, with its elements giving the outcome probabilities of each eigenstate. But it does not furnish the mechanism by which, is a given measurement, the system would "collapse" in the corresponding subspace. 4 The proof of (3.11) follows multiplying both sides by e + i H(t−t 0 ) and differentiating with respect to t.

DERIVATION OF MASTER EQUATIONS
where V (s) is the Heisenberg evolution under H (interaction picture), that is, V (s) = e +(i/ )H s V e −(i/ )H s . Taking the Hermitian conjugate of (3.12) we get a similar expression for U + (t; t 0 ). Both results can be inserted now in (3.10), giving the formal time evolution of ̺ tot (t) = U ̺ tot (t 0 )U + . As we seek for a differential equation, we can differentiate ̺ tot with respect to t. Keeping only terms to order V 2 , one gets: where the operators without argument are assumed in the Schrödinger picture, whereas V (s − t) is the formal Heisenberg evolution under H. The tilde in ̺ tot means unperturbed evolution under H, that is, Now we have to trace over the bath to get the evolution equation for ̺, Eq. (2.24). It is convenient to write the system-bath coupling as (F = F + , E = E + ) with F acting only on the system's Hilbert space and E over the bath's. 5 Now we will assume that at the initial time t 0 the density matrix ̺ tot is given by the following tensor product (decoupled initial conditions): Besides, ̺ b would be the equilibrium density matrix of the bath in the absence of the system, i.e., That is, the bath is considered at thermal equilibrium before "meeting" the system [see paragraph before (2.11) and (2.12)]. The factorized form (3.15) is difficult of justify, but it is an assumption required to proceed further with the derivation. To be more reassured, think that in this work we are interested in the long time behavior, where the system would have "forgotten" any initial conditions. After all this preparations, let us trace over the bath's degrees of freedom, Tr b (· · · ). Thanks to (3.15), we know that ̺ tot (t) = ̺ s (t) ⊗ ̺ b (t) because H s and H b act on different Hilbert spaces ([H s , H b ] = 0). Besides, ̺ b ∝ e −βH b → d̺ b /dt = 0 and ̺ s (t) = e −(i/ )Hs(t−t 0 ) ̺ s (t 0 )e +(i/ )Hs(t−t 0 ) . Therefore, for the first order term of (3.13) we have where · · · b ≡ Tr b (· · · ̺ b ). In most cases of interest the bath is such that E b = 0, and the above term drops. For the first term in the integrand of (3.13) we have: 5 The more general form P j Fj ⊗ Ej is easily handled adding summations here and there.

CHAPTER 3. DYNAMICS OF OPEN QUANTUM SYSTEMS
where is the correlator of the bath operator. This shows a clear analogy with the classical case [see Eq. (2.14)], so we have used the same letter to denote it. Notice that to derive (3.17), the factorization (3.15) has been essential; then ̺ tot (t) = ̺ s (t) ⊗ ̺ b (t). As the coupling H sb is factorized too, Eq. (3.14), we were able to take the parts acting on the system out of the trace action. Proceeding in the same way with the term [V (s − t), ̺ tot ]V in (3.13), we finally arrive at with the following relaxation term Note that we have replaced ̺ s by ̺ inside R, as they differ in second order terms, while R is already of second order in V . Now we can be more specific on what we mean by weak interaction. The equation obtained, Eq. (3.19), was based on the expansion (3.12). There we assumed that V was a perturbation on H, so that it made sense to retain only the first few terms in the expansion. Well, but (3.12) would only be valid at short times, V t/ ≪ 1, because letting |t − t 0 | → ∞ the integral could become very large and break down the perturbation theory [34]. So, is it that the master equation only makes sense at short times? This would reduce considerably its applicability (let alone the description at short times being biased by the decoupled initial condition). Fortunately, the fact is that we implicitly assume the existence of correlation time τ c obeying

Heisenberg evolution
In much the same way as we obtained an equation governing the evolution of the reduced density matrix, we can obtain the corresponding equation for any operator A within the Heisenberg picture. The Heisenberg evolution for A is given by A(t) = U + (t, t 0 )A(t 0 )U (t, t 0 ). We can use the same approximate expressions (3.12) for U and U + and then differentiate with 3.4. APPLICATION TO THE BATH-OF-OSCILLATORS MODEL respect to t. Keeping terms to second order in V , we would get an equation dA/dt = (· · · ) analogue to (3.13). Assuming then (i) the factorized initial conditions (3.15), (ii) that the operator A only acts on the system's variables A = A ⊗ I b (I b is the identity in the Hilbert space of the bath), (iii) taking into account that A(t) = Tr b (A(t)̺ b ), and finally (iv) taking the partial trace Tr b (· · · ̺ b ), one arrives at with the following relaxation term Here the operators without time argument are understood to be evaluated at t. Note that R in the Heisenberg equation is structurally different (simpler if you want) from the R obtained for the density matrix ̺, Eq. (3.20). Remember that now we are in the Heisenberg picture, so that ̺ tot (t) = ̺ tot (0) = ̺ s ⊗ ̺ b ; we also exploited this to trace over the bath above. Finally, from the A solving the above equation one can get the desired average as A = Tr s (A̺ s ) .

application to the bath-of-oscillators model
Let us now particularize the master equation (3.19) to the bath of oscillators model (2.4). Note first that the Hamiltonian (2.4) can be rewritten as: where the counter-term 1 2 F α |u α | 2 /ω 2 α has been absorbed into H s [see the discussion after (2.4) and (2.9)]. In this form, the model has the structure H sb = F ⊗ E for which we derived the master equation (3.19). (The case of Heisenberg evolution (3.23) is analogous.) In the quantum case it is convenient to rewrite (3.26) as: where b + α , b α are bosonic operators (creation and destruction) with the new constants c α = /(2ω α ) u α . We recognize the last term in (3.27) as the standard quantum oscillator Hamiltonian ωb + α b α . Let us now particularize the equations of the previous section. The correlator (3.18) follows from the time evolution of the bosonic operators b α ∝ e iωαt as C(τ ) = α |u α | 2 n α e +iωατ + (n α + 1) e −iωατ where n α = 1/(e β ωα − 1) is the Bose distribution. Introducing now the spectral density (2.13) as in the classical case: we can write the correlator as This renders the classical result (2.14) in the limit → 0, letting coth(β ω/2) → 2/(β ω).
Here we have added to R the term coming from the counter-term in (3.26), with γ(0) = α u 2 α /2ω 2 α which can be expressed in terms of the so-called damping kernel γ(τ ) = ∞ 0 dωJ(ω)/ω cos(ωτ ). That term cancels with iγ(0)[F 2 , ̺] arising after integration by parts in R. Therefore, as discussed in the classical case, all renormalizations on the system have been accounted for from the outset.
Let us close with a few comments. Equation (3.30) is sometimes criticized for not being of Lindblad type, and hence for not ensuring positivity in the time evolution of ̺. 6 In appendix D. 3 of PAPER IV we addressed this point from a physical rather than a rigorous point of view. Let us summarize our position. Equation (3.30) was derived under several approximations, the most important being γ τ c ≪ 1; outside its validity range the equation would violate positivity and what not. Besides, the derivation assumed decoupled initial conditions at t 0 = −∞ [see Eq. (3.15)]. Therefore, taking arbitrary initial conditions at t = 0 (which are the only cases of reported violations of positivity), it seems natural that the approach may produce nonsense.
So then, can we only study the evolutions under those conditions? Yes and no. Imagine we want to study the evolution from an equilibrium state at t = 0 of system plus bath. If the total system is in some sense "ergodic" it does not matter the initial condition we used at t 0 = −∞ (for instance, factorized, or other); at t = 0 we will have our Boltzmann  "preparing" the initial conditions dynamics studied H s Figure 3.2: Scheme of the protocol discussed in the text to work consistently with the master equations. We assume that bath and system first met a long long time ago: long enough so as the decoupled initial conditions no longer matter. Then, from t = 0 onwards we manipulate H s with fields, etc., to study the response of the system. initial condition. From that time onwards, we can manipulate the system with fields, etc.
(that is, modifying via H s ) and study the corresponding response (see Fig. 3.2). And all this because we are respecting a continuous process coming all the way from t 0 = −∞. 7 Therefore, if kept within its range of validity (weak coupling), the equation would preserve positivity, hermiticity, etc. to the approximation considered. 8

a few examples
Within the oscillator bath model, there are some cases where we can write the master equation (3.30) explicitly.

when the eigenstates of H s are known
If we know the eigenvalues ǫ n and eigenvectors |n of the system Hamiltonian, H s |n = ǫ n |n , we can write explicitly the equation for the components of the density matrix ̺ nm in that basis. To this end, we sandwich both sides of (3.19) as n|∂ t ̺|m , and employ: where ∆ nm = ǫ n −ǫ m . After some algebra (done with care in appendix D.1 of PAPER IV), one obtains with the explicit relaxation term Here the relaxation rates are with the part γ(0) = α u 2 α /ω α coming from the counter-term [see after Eq. (3.31)]. When the correlator C(τ ) is given by Eq. (3.29), one can prove the important relations The first is the detailed balance condition. The second is a consequence of having accounted for the renormalization (and will be exploited later on). The equation obtained is of the Redfield type, like those introduced for the study of magnetic resonance [36].

RWA
The "rotating wave approximation" is common in many branches of physics. Roughly speaking, the RWA consists in dropping the terms which "rotate" much faster than the characteristic frequencies of the problem addressed. As an illustration imagine we have terms like e i(ω 1 +ω 2 )t and e i(ω 1 −ω 2 )t (with ω 1 , ω 2 > 0). Then one argues that e i(ω 1 +ω 2 )t oscillates much faster than e i(ω 1 −ω 2 )t and can be replaced by its average (zero). Let us see how the RWA is implemented in our problems. Let us jump into the "rotating frame": ̺ nm (t) = ̺ nm e −i/ ∆nmt where the densitymatrix equation (3.34) reads: The crudest approximation assumes that the terms e i/ (∆nm−∆ n ′ m ′ )t with ∆ nm −∆ n ′ m ′ = 0 revolve in time much faster than the characteristic relaxation scales of the system, so they are averaged to zero. That is, we discard all terms with ∆ nm − ∆ n ′ m ′ = 0. If we have a system spectrum without degenerations or regularities (see below), that is, where all ∆'s are different, 36]. Therefore, in such a RWA the master equation splits into diagonal and off-diagonal elements as follows where we have returned to the unrotated basis. 34

A FEW EXAMPLES
It is not difficult to see that R nn ′ nn ′ is real: mmnn (as ̺ is Hermitian). But in the RWA one only retains the real parts of R nnmm , arguing that the imaginary parts reflect an extra renormalization. Notice that the equations for the off-diagonal terms are now uncoupled, and their solution is simply Therefore, we just need to solve the equation for the diagonal terms. Then analytical advances are sometimes possible, and in the worst case, it will be computationally "cheaper". Furthermore, equation (3.39) is of Lindblad type [28] (by some chance), and the equilibrium solution ∂ t ̺ = 0 results to be the canonical equilibrium ∝ e −βHs (see next chapter), what makes it more "user friendly". As for our discussion of decoherence, given that the W (∆ nm ) > 0 (see PAPER IV), the R nnmm give directly the inverses of the decoherence times. In this case the basis where the density matrix becomes diagonal, the pointer basis, is that of eigenstates of H s (this is related with the equilibrium solution being ∝ e −βHs ). Finally, the equation for the diagonal terms has the same form as the Pauli master equation, but there the off-diagonal dynamics is neglected.
The RWA is many times difficult to justify (for a illustration of the limitations see [37]). For instance, we could have systems with ∆ nm − ∆ n ′ ,m ′ very small. Or it could also happen that ∆ nm − ∆ n ′ ,m ′ = 0 in more cases than those considered, for example, if there is a degeneracy, or if H s has an equispaced spectrum ∆ nn±1 = cte. In the latter case, we would have ∆ nm − ∆ n ′ ,m ′ = 0 also when n − m = n ′ − m ′ for any shift p in the indexes. Then, we are not allowed to drop the R nn+p,mm+p terms in (3.34) and we write We will refer to this equation as "improved RWA", and it will be invoked in the last chapter.

semiclassical or high temperature equation
This is another instance where the relaxation term can be written in a simple form. Let us start having a look at the formula (3.29) for the correlator. Although the frequency integral extends to infinity, we have J(ω > ω c ) = 0 with ω c the bath cutoff. 9 Then, at high enough temperatures, β ω c ≤ 1 one can approximate coth(β ω/2) ∼ = 2/(β ω). Now let us consider a low frequency part J(ω) ∼ mγ ω, which corresponds to the Markovian case discussed in the classical limit (see Sec. 2.3.1). Under these conditions, the correlator (2.14) is left as Therefore the relaxation term reduces to This is the semiclassical or high-temperature equation. It has the advantage of not having an integral in R (time local) and hence is easier to handle. If we now particularize the above R to the case of a mechanical Hamiltonian H s = p 2 /2m + V (x), with the particle linearly coupled to the bath coordinates, F = −x, we get with [A, B] + = AB +BA the anticommutator. This is nothing but the "famous" Caldeira-Leggett equation, originally derived with the path integral formalism [9]. Let us close with a few remarks on the validity range as compared to the original equation (3.30). First, we used J(ω) = mγ ω, with γ having units of t −1 and being related with the coupling strength [see Eq. (2.13)]. On the other hand, for such J(ω) the bath correlation time is τ c = β (see [28,Sec. 3.6]) and hence the validity range from Eq. (3.22) would be β γ ≪ 1. But on top of this restriction, Caldeira-Leggett put the β ω c ≤ 1 condition. Then, one fears that a priori the validity range would be quite restrictive. However, in the case of the harmonic oscillator, H s = p 2 /2m + 1 2 mΩ 2 x 2 , exact quantum master equations have been obtained [38,30]. And for this case Unruh and Zurek [39] checked empirically that, except at short times, the restrictions on the hightemperature equation could effectively be relaxed to ω c /Ω ≪ e k B T / γ . Although limited to the harmonic oscillator, this is anyway good news on the applicability range of the Caldeira-Leggett master equation.

summary
We are aware that this chapter has been terse. Let us summarize the main points.
First we have discussed the concept of quantum decoherence. The decoherence program is an example of how quantum theory can save face. In a way similar to how reversible dynamics could account for the irreversibility found in nature, by invoking the interaction with the environment, quantum mechanics manages to explain the "non-observation" of certain odd superpositions in the macroscopic world (for the good of cats and other creatures). The reason, we are told, is the extremely fast lost of coherence when we only look at our system, through the notion of partial tracing [Eq. (2.24)].
Then we discussed the description of the dynamics of open systems by means of a master equation for the reduced density matrix ̺. Those master equations are valid in the range γ τ c ≪ 1 with γ quantifying the system-bath coupling strength and τ c the bath correlation time. The prime difficulty with such equations is in dealing with the time integral in the relaxation term R [Eq. (3.20)]. We saw how this can be handled without problems when the spectrum of H s is known, getting a Redfield-type equation [Eq. (3.35)].
Finally, we discussed a couple of approximations widely used in the specialized literature, like the rotating wave approximation [Eq. (3.39)] and the high-temperature approximationà la Caldeira-Leggett [Eq. (3.45)]. Both equations have validity ranges, in principle, more restrictive than γ τ c ≪ 1, although when they can be used result in a great simplification.

equilibrium properties
Up to this point we were not concerned with the equilibrium properties in open systems. The equations discussed, classical or quantal, describe an irreversible dynamics in which the system in contact with the reservoir exchanges energy (chapters 2 and 3). Waiting long enough the system will equilibrate to the bath, reaching a stationary regime. Which is the state reached is the question that we will try to answer in this chapter.
We will start revisiting the classical case. Then we introduce the thermodynamic perturbation theory formalism, and expand the partition function in powers of the systembath coupling, getting explicitly the partition function Z to second order. After this, we will check that the master equations derived in the previous chapter are thermodynamically consistent. That is, that they have as stationary solution the equilibrium distribution obtained here to the same order in the bath coupling. We conclude discussing how to compute equilibrium quantities in open systems, and illustrate this with the example of the harmonic oscillator.

the classical limit
The best way to address the equilibrium properties starting from the dynamics is to resort to the Fokker-Planck equations. The necessary condition for W eq to be a stationary solution is ∂ t W eq = L cl W eq = 0, where L cl is a compact writing of the Fokker-Planck differential operator. If we focus on the particle case [Eq. (2.16)], we can easily check that L cl e −βHs = 0. 1 Furthermore, it happens that any initial condition, in the limit t → ∞, tends to the same solution [14]. Therefore, the system relaxes to the equilibrium Gibbs-Boltzmann distribution W eq ∝ e −βHs . The dynamics of classical open systems places us in a comfortable theoretical framework where the Fokker-Planck equations capture the approach to equilibrium.
Physically, what we do is to shift everything by a distance ∝ x. As the integrals extend from −∞ to +∞ their limits do not change. Compare with the analogy of the Bohr-Van Leeuwen particle in a magnetic field (final comment in Sec. 4.2.1).
Let us reexamine this from another point of view. Although in an effective/reduced way, the Fokker-Planck or Langevin equations describe the total dynamics of system plus bath (2.4) (see the discussion on irreversibility in Sec. 2.4). But the total system is a macroscopic entity, and hence we can apply the Gibbsian hypothesis. That is, at sufficiently long times the averages can be computed with the Gibbs distribution for the whole system As we are discussing the reduced distribution, we integrate out the bath, getting or in other words, Z tot = Z b Z s , with Z b = dP α dQ α e −βH b and Z s = dpdx e −βHs . This result has relied on the integrals over the bath being Gaussian (see Fig. 4.1) so that the change of variables Q α + u α F/ω 2 α → Q α gives W eq = e −βHs /Z s . This is the reason behing (4.2) giving the reduced distribution equal to that of the system without couping. Now we understand that e −βHs /Z s is stationary solution of (2.16), since the Fokker-Planck equation is nothing but the evolution of the total system with the bath variables integrated out.
We remark that, although it results natural to accept that the system relaxes to the distribution e −βHs , from the point of view of (4.2) this is a "happy" coincidence. We mean that the result is "model dependent" and for other models (coupling nonlinear in the bath coordinates, or other kinds of baths), the reduced distribution could be different from e −βHs . Otherwise, every system with several degrees of freedom would factor into the product of reduced distributions for each degree of freedom.
In fact, this result is so particular that, as we will soon see, does not hold for the same Hamiltonian in the quantum case, and the system relaxes to a distribution different from the canonical one. What is the meaning behind this? In the books of statistical mechanics one assumes that the body is macroscopic and that the interaction with the bath is somehow mediated by its "surface". As the surface L 2 to volume L 3 ratio goes to zero in large systems ∝ 1/L, the interaction with the reservoir H sb is put aside, leaving the Gibbs distribution [40]. But for the open systems we are treating here we cannot neglect H sb . System plus bath, on the other hand, constitute a macroscopic object. And here it is where we can assume Gibbs. But then nothing ensures that the reduced distribution is going to be ∝ e −βHs . We insist again, the result (4.2) is quite particular.

the quantum case
We now address the quantum problem. We will try to answer to the following questions • Is the reduced distribution the canonical one for the system without coupling?
To which we have already advanced the answer (previous paragraph), and, if it is not so, how is that distribution?
• Is the reduced density matrix a stationary solution of the master equations discussed in chapter 3? That is L q ̺ eq with L q a compact way of writing the master operator in (3.30).

thermodynamic perturbation theory
In principle we do not have a practical expression for e −βHtot (much as in the dynamics we did not know an exact differential equation for the reduced density matrix). Therefore, we will assume again that the system-bath coupling is weak. This will allow us to perform the equilibrium calculation to the same perturbative order as we did in chapter 3 for the master equation. Besides, with this calculation we will gain some insight in the way the first corrections to the uncoupled distribution emerge. In order to obtain an approximate expression for e −βHtot we will start again from the Kubo identity (3.11). As in Sec. 3.2 we will identify H ≡ H s + H b and H sb ≡ V and expand to the second order: This is nothing but the evolution (3.12) in imaginary time, with V (−i σ) = e σHs V e −σHs . But note that here the expansion is in powers of βV . We now trace with respect to the reservoir and apply the result to the bosonic bath model (as in section 3.4), getting Here Z s = Tr s (e −βHs ) and the "rest" R which is a byproduct of the counter-term in (2.4). The result above is already properly normalized, Tr(̺) = 1, and Z (2) is essentially the second order correction to Z tot = Tr(e −βHtot ): The superindexes (2) tell that this result is valid to second order. The first order terms disappear as they did in the dynamics since in the bosonic bath (3.27) we have E b = 0. We have thus answered the first question, showing that the reduced distribution, contrary to the classical counterpart, does not agree with the canonical result at zero coupling.
We close this section with a few more remarks on the result obtained. First notice that (4.10) That is, the equilibrium quantities, as computed from the partition function, will depend on the bath through the imaginary parts of the coefficients R nn ′ mm ′ [Eq. (3.36)]. The real parts are associated to relaxation and decoherence times (section 6.2), that is, to how fast we reach the stationary solution. But now we are in the statics, and we do not care about the speed of the approach, but about the final state instead; it then makes sense that the equilibrium results do not depend on the real parts of R nn ′ ,mm ′ .
It is instructive to see that in the limit β → 0, both Z (2) and the term R go to zero, recovering the classical result of damping-independent equilibrium properties in a Gaussian bath. Consistently, in the high temperature limit the quantum corrections to the canonical result disappear too (Z tot = Z s Z b ). 2 We conclude with an analogy to a famous problem: the absence of orbital magnetism in classical statistical mechanics. Let us start from the Hamiltonian of a charged particle in a magnetic field and compute the partition function Z = d q d p e −βH . As the momentum integrals are Gaussian, on making the variable change p i = p i − eA i , the vector potential disappears from the final result, producing a zero magnetization M z ∝ ∂ ln Z/∂B z ≡ 0. This is the content of the Bohr-Van Leeuwen theorem.
In our case we have Gaussian integrals too, but it is u α what disappears from Z (instead of A), so that the equilibrium quantities become independent of the system-bath couplings. Quantum mechanically (we are thinking in the obtainment of Z from the density matrix), the Gaussian argument does not hold. Correspondingly, the equilibrium quantities do depend on H sb , much as the charged particle developed a non-zero magnetization (Landau).

stationary solution of the master equation
Let us address now the second question, that is, let us see if the damping-dependent equilibrium distribution (4.7) is a stationary solution of the quantum master equation (3.30). We would not hide that the answer is in the affirmative. It should be so if we did things properly, since [e −βHtot , H tot ] = 0, while the master equations are just the Von Neumann evolution of the total system with some tracing over the bath degrees of freedom (remember the classical discussion in Sect. 4.1). What we aim is to prove here is that the equality still holds when e −βHtot and H tot are calculated order by order. Let us see this. 3 As all our calculations are perturbative (both in the dynamics and in equilibrium), we will check L q ̺ eq ? = 0 order by order. Taking into account the decomposition (4.7) of the density matrix ̺ = ̺ (0) eq + ̺ (2) eq we obtain the following consistency conditions On writing R = R nm |n m| we have [H s , R] = ∆ nm r nm |n m| with ∆ nm the energy difference between the n and m eigenstates of H s . Therefore, the equality (4.13) follows for the diagonal elements (n = m), given that R[̺ (0) nn ]|n n| = 0 (we are omitting the marker "eq" in ̺ to alleviate notation). Finally, is not difficult to prove the relation (4.14) showing that the second condition is obeyed too. So we have just managed to prove the thermodynamic consistency of our master equations. This is also relevant for the linear-algebra eigenvalue problem associated to the master equation. Solving it, one of the eigenvalues will be zero, corresponding to the stationary solution (4.7), as we have just checked. As for the rest, we expect they having negative real parts, ensuring convergence at long times to the stationary solution. Positive real parts would imply a divergence with time which would not be physically acceptable.
Well, at this point we have already answered our two questions. That is, the stationary solution of the master equation (3.30) is the equilibrium result (4.7), and this does not agree in general with the canonical distribution at zero coupling. In other words, although by now it feels quite natural, the equilibrium properties of quantum open systems do depend on the damping strength. When we started this chapter this would have struck us, due to the intuition from the classical case and the Gibbsian prejudice. Indeed, in many works it is still stated that the system would approach the canonical e −βHs . Here we have argued that in general this is not the case. 4

computing thermodynamic quantities
Once reached this point, we could use our formula (4.9) for Z and compute thermodynamic quantities. To this end, it will prove convenient to redefine the free energy as follows (all other quantities will follow from this) Here F tot = −T ln(Z tot ) with Z tot = Tr(e −βHtot ) the partition function of the whole system, while F b = −T ln(Z b ) with Z b would be the partition function of the bath uncoupled from the system Z b = Tr b (e −βH b ). That is, we are redefining the zero of free energy, removing the contribution of the bare bath. Besides, in this way our definition above for Z = Z tot /Z b allows obtaining this F by the custom rules The idea will be to use now for Z the perturbative result (4.9). From the definition of F we can derive energies and entropies in the usual way Thus, our convention gives U = U tot − U b and S = S tot − S b since all quantities follow from F by linear operations. Then, in processes not involving the bath, we are just redefining the zero of all thermodynamical quantities. But even in the bath is involved, as we have merely extracted its bare contribution, it can be easily restored to the final results.
Let us illustrate what these definitions would give in a calculation of the specific heat C. Imagine we want to find the C of our system, which as usual in this work, is an open system and interacts, for example, with the phonons of the lattice where it "lives". To that end one applies a controlled heat pulse with a resistor, and measures the temperature change C tot = ∆Q/∆T . To isolate the contribution of the system from other sources (phonons), one can repeat the measurement in the absence of "system" (for instance, if it consists of magnetic ions, we measure a compound without spins, but having the same lattice). This would give the bath specific heat C b , which experimentalists routinely substract from the original C. But the quantity C = C tot − C b is just what directly comes from our redefinition of the energy, when computed in the standard way C = −∂ T U.
Instead of getting the thermodynamical energy as in (4.17), we could have assumed that the relevant system energy is H s = Tr(̺H s ). But in general this would not agree with our U = H s . 5 The bare average of the system Hamiltonian H s misses the important contribution of the system-bath coupling. And we do not mean, as above, the additive part of the bath alone, but the work done to "insert" the system in the bath. Indeed, the use of H s to discuss thermodynamical principles instead of U can produce "violations" of the second law (see the lucid discussion of Ford et al. [42,43,44,45]).
On the other hand, for observables not depending on H sb -we have in mind the magnetization of a paramagnet where the coupling does not depend on the magnetic field-we have That is, in this case both the definition from the free energy, or directly from the reduced density matrix ̺, provide us with two equivalent ways of computing the magnetization (that we will use in chapter 8).

example: the harmonic oscillator
Let us apply the thermodynamic formalism discussed to the example of a harmonic oscillator H s = p 2 /2m + 1 2 mΩ 2 x 2 , bilinearly coupled to the bath F = −x. This exactly solvable model has been thoroughly studied in the literature [1,18,46]. That allows us to compare our perturbative results with exact ones. As F = −x can be expressed as ∝ (a + + a) and the harmonic oscillator's spectrum is equispaced ∆ nn±1 = ∓ Ω, we have with ǫ n = Ω(n + 1/2). Doing the geometric sums we arrive at with Φ ± ≡ Φ( Ω) ± Φ(− Ω) and the Φ function from (3.36). Now we can obtain F directly via (4.16), and by derivation U, S, C, etc, as in (4.17). To illustrate, we obtain the following specific heat in the high-temperature range Ω/k B T ≪ 1 while in the regime of low temperatures Ω/k B T ≫ 1: Both results agree with the expansion to first order in the damping of the exact results by Ingold and Hanggi [47]. 6 Notice that in the low temperature regime, the dissipative corrections give C ∝ T . This is is contrast with the celebrated γ = 0 result where C drops exponentially at low temperature, due to the gap between the ground state and the first excited state. To make sense of this, we repeat again that now the total system is our system plus the bath, which removes the gap (the reservoir has infinitely many oscillators with a quasi-continuous distribution of ω α 's).
Another way of looking at this is resorting to the definition of density of states That is, the inverse Laplace transform of the partition function. For the central harmonic oscillator, Hanke and Zwerger [49] obtained d(ǫ) exactly, finding that between the ground state and the first excited state there exists a finite density of states (that is, there is no gap, see Fig. 4.2) 6 We have used an Ohmic bath spectral distribution J(ω) ∼ mγ ω regularized with an exponential cutoff ωc. Then we know the explicit form of Φ: in terms of the ψ function [48] and the bath cutoff ωc.

final comments
Performing the RWA in the master equation one gets Eq. (3.39), whose stationary solution coincides with the canonical. 7 That is, the equilibrium distribution at zero coupling.
Possibly this result has contributed to the "belief" that the system always equilibrates towards e −βHs in the quantum case as well. But in general this is not so, and here we have studied the actual form of the equilibrium distribution, and its relation with the stationary solutions of the quantum master equations.
In the low damping range, the corrections computed here to the thermodynamical quantities, though sensible, are probably small and not easy to detect experimentally. In a way, this contrasts with the situation in dynamics, where weak coupling can produce sizable effects via relaxation and decoherence. Anyway, from the conceptual point of view the rigorous treatment of the equilibrium properties contributes to the understanding of those dissipative contributions, absent in the classical case.

working in phase space
We conclude the general discussion of the formalism of open systems with the reformulation of the master equations for ̺ in phase space. The phase space approach, together with canonical quantization and path integral are three equivalent ways of quantizing [50]. In this formalism the theory "lives" in the space defined by the canonical conjugate variables of Hamiltonian dynamics (phase-space). Given that quantum and classical theory share now the same mathematical space, the formalism provides an ideal framework to study the transition from quantum to classical (or the other way around). In this frame we will be able to relate explicitly the classical and quantum theories of dissipation (chapters 2 and 3).
We will motivate this chapter as Wigner did, in 1932, in his seminal article [51]. We will introduce historically the development of the theory for spinless particles. Afterwards we will provide a more modern view, following the postulates of Stratonovich to map Hilbert operators and functions in phase space. We then address the problem of systems with spin. In both cases (particles and spins) we transform the master equations into phase space. Then the corresponding Fokker-Planck equations are obtained in the limit " → 0".
Keeping in mind the spirit of Caldeira and Leggett, who posed the theory of dissipation as a quantization problem [10], this chapter closes the circle, getting now the classical limit from the quantum case. The formalism discussed here will be very useful in Chap. 7, where we will solve quantum master equations in phase space by adapting techniques developed to solve classical Fokker-Planck equations.

the problem
In classical mechanics the dynamics is given by the Liouville equation: in the one-dimensional particle case (V = R 2 ) while v = (θ, φ) for classical "spins" (V = S 2 with S 2 the sphere in three dimensions).

49
In quantum mechanics, the Liouville equation is substituted by Von Neumann equation: where now H and ̺ are Hermitian operators in some Hilbert space.
In spite of the formal similarities { , } ↔ [ , ] these two equations seem difficult to relate. We mean, it is not straightforward to take some appropriate limit in Von Neumann equation and arrive at Liouville's. They live in different mathematical spaces! In this sense it turns to be involved to study the classical-to-quantum transition, or the quantum-toclassical, in a continuous way.
To bypass this problem we resort to the phase-space formalism, in which both classical and quantum mechanics "live" in the same space. Formulating the quantum theory as closely as possible to the classical case would allow to study their crossover, and also see how the corrections to the classical case are built; as well as studying the problems from another point of view. In particular, one can export notions and tools of the classical case to the quantum world.

a constructive tale
We are going to review the construction of the phase-space formulation of quantum mechanics from a historical perspective.

Wigner's function (1932)
In that year Wigner introduces the following "distribution" [51]: where ψ is the quantum-mechanical wave function. We have added the subindex to W ψ in order to distinguish it from the classical probability distributions (note that W ψ = W * ψ ). Next introducing position "eigenstates" |x , one defines the transform of any operator A as follows With these definitions the expectation value of A in the state |ψ can be written as an "average" over phase space: Thus we see that W ψ plays the role of a distribution. Advancing results to fix ideas, W A (x, p) is the classical dynamical variable A(x, p) corresponding to the operator A.

A CONSTRUCTIVE TALE
The transformation properties of W ψ are as follows [52]. Applying ψ(x) → ψ(x + r) We thus see that W ψ is real, normalized, respects the symmetries, and the averages are computed like in classical mechanics (same for W ̺ ). 1 But the analogy is not fully complete, some differences should show up! W ψ is not positive, and besides W ψ ≤ π/ . That is, the value at any point (x, p) of the quantum distribution W ψ is bounded. This precludes, for instance, that one could write W ψ (x, p) = δ(x−x 1 )δ(p−p 1 ) for a distribution localized at (x 1 , p 1 ). As Moyal pointed out in 1949, this reflects the uncertainty principle, forbidding well defined x and p "simultaneously" [53].

the 40s: Groenewold & Moyal
In 1946 Groenewold expresses the transform of the product of two operators as follows where the arrows indicate on which W the derivatives act upon.
In modern terminology the transform of the product defines an inner product in V, called star product: Then the algebra of the product of operators (composition) is transformed in the ⋆-product in V, inheriting its properties. Namely, it will be associative, W C(BA) = W (CB)A , but in general non-commutative, W AB = W BA . Von Neumann equation gets transformed in phase-space as Introducing now H s = p 2 /2m + V (x) and using Groenewold's result (5.7) we obtain Von Neumann dynamics in phase space One can proceed analogously with ̺ (mixed states) introducing the analogue to W ψ (5.1): So the calculation of averages (5.3) now reads: To abbreviate we removed the subindex ̺ for the density matrix transform, and written V (2s+1) ≡ ∂ 2s+1 x V . Besides, we have extracted the s = 0 term from the sum (V ′ ∂ p ) gathering all -dependent terms. In this way the formal limit → 0 recovers the classical Poisson bracket: With the glasses of phase-space, the equation for the quantum dynamics is built from the classical one simply adding more terms to the Poisson bracket. 2 Note that for the harmonic oscillator V = mΩ 2 x 2 /2 we have V (2s+1) = 0 when s ≥ 1, so Von Neumann's equation in phase space is identical to the Liouville equation. What is going on? Quantum and classical mechanics are the same for the harmonic oscillator? No, what happens is that in the classical and quantum cases the admissible solutions are different. Recall that W could take negative values and in addition is bounded [55].
In 1949 Moyal rounds up the work noting that Wigner's distribution (5.2) is just the inverse of Weyl's prescription for symmetric quantization ordering [40], that is: Besides, he computes the uncertainty relations in phase space [53]. Therefore, reached this point -calculation of observables, ⋆-product, and dynamics-the formalism in phase space was completed, constituting an alternative to the standard canonical quantization in Hilbert space. From the 50s onwards many applications were developed of this formalism. Many of them exploited its analogy with the classical formalism, adapting its schemes and tools to the quantum world [52,56,57,58].

Bopp (1961)
In 1961 Bopp introduces an operationally better way of computing the star product where B(Q, P) acts on the functions in V, with the same functional form a B(x,p) but withx &p being replaced by the differential operators Q and P. As an aside, note that these operators obey in V the same commutation rules as their Hilbert analoguesx and p, that is, [Q, P] = i . 3 Let us apply Bopp's result to Von Neumann dynamics. To this end we use ̺ = ̺ + and H = H + , whence (H̺) + = ̺H, together with W A = W * A + for adjointing. Then,

OTHER SYMMETRIES AND STRATONOVICH'S POSTULATES
W H̺ = W * ̺H , and we obtain ∂ t W = 2Im H(Q, P) W . (5.14) The Hamiltonian H is now a function of the Bopp operators Q and P which act on the functions in R 2 . The equation above is a compact way of writing the overwhelming (5.10).

Caldeira-Leggett equation in phase space (1983)
Let us apply the formalism to the dynamics of open systems. This will prove useful in Chapter 7. We want to transform into phase space the master equation ( (5.17) We have written the equation is such a way that the first three terms correspond to the classical Klein-Kramers equation (2.16). And the last term contains the purely quantum terms coming from the closed evolution. It is important to remark that the terms (5.15) and (5.16) are the diffusion and dissipation terms in the classical case. That is, the Caldeira-Leggett equation is nothing but the Von Neumann evolution augmented with the classical irreversible terms. This identification has been possible thanks to the phase space formalism. Who

other symmetries and Stratonovich's postulates
We have seen how the phase-space formalism developed from the introduction of the Wigner distribution. That is, everything followed from the map between operators and dynamical variables in phase space [Eqs. (5.2) and (5.4)]. But the map is restricted in the sense of "only" being applicable to "particle systems", we mean, to the algebra [x,p] = i , corresponding to the Poisson bracket {x, p} = 1. Then, how could we proceed if we want to tackle problems with other commutation rules? We have in mind a phase-space description for systems with angular-momentum type algebra [S i , S j ] = i ǫ ijk S k , corresponding to the classical Poisson bracket {m i , m j } = gǫ ijk m k of Eq. (2.18) (g is the gyromagnetic ratio, and ǫ ijk the completely anti-symmetric unit tensor).
The answer was given by Stratonovich [59], completed by Varilly & Gracia-Bondia [60], and later generalized by Briff and Mann for any finite Lie group [61]. In brief, the idea is to find a map between Hilbert operators and phase-space functions by imposing the "physically" reasonable properties obeyed by the Wigner function (section 5.2.1).

the Stratonovich postulates (1956)
Mathematically one starts with a physical system having a given "symmetry group" ([x,p] = i , or [S i , S j ] = i ǫ ijk S k , etc) corresponding to a classical phase space V. 4 Then one seeks for a map between the Hilbert space and V obeying the following Stratonovich postulates:

(i) linearity
A → W A is linear and bijective.
(ii) "reality" where G is our symmetry group and A g := U(g −1 )A U(g) with U the group generators in Hilbert space.
In [60,61] it was proved that these properties determine uniquely the map. The map is bijective (i) and hence invertible. Condition (ii) ensures that W A is real if A = A + [see after (5.2)]. Postulate (iii) is simply the normalization condition [see below (5.3)]. Linearity plus (iv) ensure that the averages are computed as in Eq. (5.3). Finally (v) expresses respect to the basic symmetries of the system. In the case of particles, these postulates produce the Wigner function discussed above.

the case of spin systems
Let us apply the formalism to spin problems. Here the dynamics is governed by a spin Hamiltonian H(S x , S y , S z ) with In the standard basis S z |m = m|m and S 2 |m = S(S + 1)|m , for m = −S, . . . , S, the spin quantum number. 5 In this case, following Ref. [60] we have the map where T lm = (2l + 1)/(2S + 1) × j,j ′ S, j; S, j ′ |S, j ′ |S, j ′ S, j| are irreducible tensors, the S, j; S, j ′ |S, j ′ are Clebsch-Gordan coefficients and Y lm spherical harmonics [34,48]. Klimov and Espinoza [62] obtained a differential form for the spin ⋆-product which reads: Here N S = √ 2S + 1 and L 2 is the angular momentum operator on the sphere In formula (5.20) the operators S ± and G inside the brackets act only inside the bracket, while G −1 on the left acts on everything to its right. Equations (5.19) and (5.20) are, respectively, the spin analogues of the Wigner map and Groenewold ⋆-product for particles.

Bopp operators for spins
At first sight the ⋆-product (5.20) looks of little use in explicit calculations. Then, in the spirit of Bopp, we would like to find a simpler representation for W B ⋆ W A . To this end we explicitly calculated in PAPER VI the transformation W S i A , i = x, y, z, finding where Hereη i (L 2 )Y lm = η i (l)Y lm for i = 1, 2, and the numerical coefficients are The vector m(θ, φ) = (sin θ cos φ, sin θ sin φ, cos θ) was introduced in (2.18) while L is 6 Finally, invoking the associativity of the star product, for an arbitrary operator B(S x , S y , S z ) we have where B(S x , S y , S z ) is the same function of its arguments as B(S x , S y , S z ), but now the S i are "classical" differential operators on the sphere. To abbreviate, instead of writing B(S x , S y , S z ) we simply write B(S).
Thus the S i in the sphere play the same role as Q = x + i( /2)∂ p and P = p − i( /2)∂ x in R 2 . For example, in analogy with the Q and P commutation, we have here [S i , S j ] = iǫ ijk S k . That is, the S i obey the same commutation rules as their S i partners (for the same reason as discussed in Sec. 5.2.3). To conclude the analogy with the particle problem, the dynamics can also be written in the following compact form [cf. Eq. (5.14)] The following formulas are of assistance

an example
Let us apply the formalism to the spin Hamiltonian studied in this work The Poisson bracket (2.18) produces in this case the differential equation: In the quantum case we employ (5.30) and (5.24) getting: The first thing one notices is that for D = 0 (isotropic spin) the differential equation is the same in the classical and quantum cases:

application to spin master equations
As we did in the particle case (section 5.2.4) let us apply the phase-space formalism to open systems. We will consider the semiclassical regime, and eventually take the limit S → ∞. Therefore, it will suffice to use the semiclassical master equation (3.44) which has been written with the conventions used in this chapter (recall footnote 5).
In order to transform the master equation to phase space, we will resort to our spin Bopp operators (5.24). We need to compute [cf. (5.15) As we did for the Caldeira-Leggett particle case, we will consider bilinear coupling, that is F = α i S i , linear in the system variables. In contrast with the mechanical case, here [H s , F ] depends always of H s , so that we cannot write the dissipative terms in a generic way. 8 Considering the simplest case H s = B i S i (isotropic spin), we obtain: There exists an analogy with classical spin problems. There a bilinear coupling induces multiplicative noise in the Langevin equation, in contrast with the particle case, where the noise enters additively for bilinear coupling. 10 Formally, this is due to the angular Poisson bracket {m i , m j } = g ǫ ijk m k which, even for F linear in m i , gives a non-constant {A, F } in the Hamilton equation (2.6) for any dynamical variable A [15]. Quantum mechanically, it is the angular-momentum commutation relations which render [F, [H s , F ]] non constant, which eventually produces the term M in the master equation.

other distributions
What we have told here is not the whole story. Canonical quantization comes associated with the ordering problem. In path integrals, similarly, it shows up in the choice of the evaluation point [63]. Analogously, in phase space this results in the ambiguity in the map between Hilbert operators and V functions -we can construct other distributions, each associated with a given ordering. In the 60s, Cahill and Glauber [64] constructed distributions (called P and Q functions) in terms of the creation and destruction operators of the harmonic oscillator. Those distributions differ from the Wigner distribution (5.1) in that they correspond to normal and anti-normal ordering, respectively [52]. For normal ordering (P-function) [cf. Eq. (5.12)]: while for the Q-function P (xp+px)/2 = xp + 2i . That is, the phase-space function xp corresponds to the operatorxp in normal ordering and topx for anti-normal. These distributions can also be generated from the Stratonovich postulates of Sec 5.3.1. Then the map gets parametrized by a superindex σ and the transformed function is denoted by W

(σ)
A ; for σ = ±1 one has the normal and anti-normal cases, while σ = 0 corresponds to the Weyl symmetric ordering discussed throughout. Postulate (iv), the trace property, is then modified as Therefore the symmetric function (σ = 0) is special, unique, in the sense that to compute the averages we only resort to functions with σ = 0. Any other σ requires transformed functions with σ and −σ.
In the time before Varilly & Gracia-Bondia, the spin distributions followed the approach of Cahill and Glauber. The construction was done, instead of using a and a + , in terms of Bloch states or spin coherent states [65] (those were the 70s). With this approach Takahashi and Shibata found the spin Bopp operators corresponding to the orderings σ = ±1 [66]; an important result they called the product theorem. In PAPER VI, under the modern postulated-based approach, we generalized the formalism to any ordering (here we have described the σ = 0 case only).

summary
In this chapter we have surveyed the construction of the phase-space approach to quantum mechanics. The particle Bopp operators, and our generalization to spin problems, have allowed us to transform easily the master equations into phase space (sections 5.2.4 and 5.4.3). With the equations so transformed we have recovered, in the classical limit, the corresponding Fokker-Planck equations, and we have seen how the first quantum corrections are generated. In this way we have connected the classical and quantum theories of dissipation of chapters 2 and 3. On the other hand, as we will see in Chap. 7, transforming the Caldeira-Leggett equation (5.17) to phase space with be instrumental for its numerical solution.

methods
This chapter is a transition between the general treatment hitherto discussed and the chapters to come, devoted to specific applications. Here we explain what we intend to calculate and how to carry it out practically.
We present first the linear response theory for the obtainment of dynamical susceptibilities. Then we use the example of the two-level system to illustrate how those susceptibilities can give us relaxation and decoherence times. Finally, we introduce the continued-fraction method to solve master equations and the practical computation of response functions.

linear response theory (general formalism)
We start recalling the protocol we discussed in chapter 3 for the study of the dynamics using master equations (section 3.4, fig. 3.2). During −∞ < t < 0 we let the system thermalize to the bath. At t = 0 we add some field, or change some parameter in the system Hamiltonian H s , and study the dynamics. Thus master equations are a natural frame to address problems where the system is perturbed, and study its adjustment to the new conditions.
In order to not modifying the intrinsic nature of the system studied, the probe can be made small enough [67]. But then the dynamics can be studied within perturbation theory. Linear response theory (LRT), along with dealing with such adjustments, provides the link with a typical experimental situation -the measurement of the dynamical susceptibilities.
So hands at work. For −∞ < t < 0 we can write the dynamics formally as where L 0 is a compact way of writing the master equation with H s = H 0 for t < 0. 1 When t = 0 the previous dynamics would have brought the system to the stationary solution (recall chapter 4) Then at t = 0 we modify H s as H s = H 0 + δ · H 1 , with δ accounting for the size of the perturbation. To first order in δ we can write: and the formal solution of the master equation is: We now ask about the evolution of the average of some observable A = Tr(A̺) and introduce: where A 0 = Tr(A̺ 0 ) is the average in the absence of perturbation. Then ∆A(t) measures how far we are from the reference unperturbed evolution.

The generic form of the change is H
The response function R follows from (6.4) in the form Therefore our task will be to calculate R(t), which is independent of f (t) within linear response theory. We will carry this out in two typical situations: the relaxation experiment, and the AC response.

relaxation experiment
Let us consider a step function perturbation f (t) = Θ(−t), so that f (t) = 1 for t < 0 while f (t) = 0 when t > 0. Note that we have reversed the "set up", saying that in the interval −∞ < t < 0 the system was "perturbed" with a constant δ, and at t = 0 the perturbation is removed (Fig. 6.1). 2 In this case A 0 = A (∞), since the "unperturbed" evolution takes place for t > 0, and we know that the system will eventually reach the stationary solution. Then (6.5) is written as That is, ∆ r A(t) describes how the equilibrium is approached as t → ∞. Hence the name "relaxation experiment" and the subindex "r" used.

AC response (frequency domain)
According to the experimentalists, it is often easier to measure the response to an oscillating probe than to conduct a relaxation experiment. In is convenient to introduce then the Fournier transform of the response function R: At t = 0 we switch on an oscillating probe f (t) = e iωt (see Fig. 6.1) and the response follows as Besides linear response theory teaches us that χ(ω) of the AC experiment is related with the response ∆A r (t) in a relaxation experiment [68, App. C]: Here χ eq = χ(ω = 0) is the response to an infinitely slow field, so that χ eq = ∆A r (0)/δ (6.12)

Bloch equations
To fix ideas and give content to the formalism, let us consider the solvable problem of the two-state system. In the basis of eigenstates |1 and |2 we can always write where σ z is the Pauli z matrix and the energy levels are ǫ 1 = −B z /2 and ǫ 2 = B z /2. The notation reflects that we cannot help thinking of a 1 2 spin with a Zeeman Hamiltonian. In this way we bear in mind a physical realization of the two-level model.
The above H s is coupled to a bosonic bath as in the model (3.27). The dynamics will be given in general by the master equation (3.34): As we set the problem ourselves, we make our life easier assuming that the rotating wave approximation can be made (section 3.5.2). Then the above system of ODEs simplifies; in particular the off-diagonal elements decouple and evolve trivially. For the two-state problem the remainder equations for the diagonal [Eq. (3.39)] are called the Bloch equations [36]: Here P 1|2 = W 12 |F 12 | 2 , with the matrix element of the coupling F ( σ) and the rate W , providing detailed balance P 1|2 = e −βB P 2|1 [Eq. (3.37)].
As mentioned the non-diagonal elements simply obey (̺ 12 = ̺ * 21 ) and its conjugate for ̺ 21 (t). Therefore one introduces the time constant T 2 ≡ 1/Γ d , measuring the speed at which the off-diagonal elements go to zero -the decoherence time.
We are left with solving the Bloch equations for the diagonal elements. Using that 1 = Tr(̺) = ̺ 11 + ̺ 22 and introducing the population difference M z ≡ ̺ 11 − ̺ 22 , we have with the relaxation constant Γ r = (P 1|2 + P 2|1 )/2. (The notation again tells that we have in mind the spin problem, with the magnetization M z ∝ σ z .) Following the literature one defines the relaxation time T 1 ≡ 1/Γ r , which measures how fast the diagonal terms equilibrate to their stationary values. And these, as we pointed out in 4.5, agree with the Boltzmann distribution M z (∞) ∝ e −βǫ 1 − e −βǫ 2 . 3 This solvable example, plus LRT, will illustrate how one can obtain T 1 and T 2 from the susceptibility curves. Given that the susceptibility is an experimentally accessible quantity, this provides a way of measuring relaxation and decoherence times.

LRT and T 1
Let us reconsider the relaxation experiment of section 6.1.1 in a two-level system. For t < 0 the system was perturbed by an extra constant δB z : At t = 0 we remove the probe, δB z = 0, and follow the evolution. Immediately before t = 0 the system would have reached the stationary state in the presence of B z + δB z . Under the RWA, we repeat again, the stationary solution is the canonical one: ̺ 11 (0) = e β(Bz+δBz ) , ̺ 22 (0) = e −β(Bz +δBz) and ̺ 12 (0) = ̺ 21 (0) = 0. The same at t = ∞, but now with δB z = 0. As δB z is suitably small we can expand to first order getting M z (0) = M z (∞) + δB z · ∂ B M z . All the quantities here are equilibrium ones; by definition ∂ B M z = χ eq is the equilibrium susceptibility, so we can write ∆M z (t) = δB z · χ eq z e −Γrt (6.22) We can plug the above relaxation function into Eq. (6.11), getting the dynamical susceptibility χ z (ω) as 4 We have written the result for χ z (ω) in two equivalent ways found in the literature. These functional forms are called after Debye's seminal work on dielectric relaxation. Let us split into the real and imaginary parts which have been plotted in Fig. 6.2. In the zero frequency limit the real part tends to the equilibrium susceptibility. Besides, the maximum in the imaginary part is located at ω = 1/T 1 . In an ac experiment, where one applies a sinusoidal field and records the response, one can directly obtain χ eq and T 1 . Another usual form of representation is the Cole-Cole plot, with the imaginary part plotted vs. the real part (Fig 6.2, right). For Debye relaxation one has Im = (1 − Re)Re, in evident notations; that is, a semicircle with radius 1/2 located at (0, 1/2). Note the marked points: the zeroes of the imaginary part are at ω = 0 and ∞, corresponding to 1 (maximum) and 0 (minimum) of the real part. The maximum of the imaginary part, Im = 1/2, is reached at ω = 1/T 1 , taking half the value of the maximum real part. This plot is used as indication of non-Debye behavior if deviations from a semicircle are found. 4 Here we use The arrows mark special points in the spectrum, namely, the peak location of the imaginary part, which gives T 1 , and the equilibrium susceptibility approached in the zero frequency limit. Right: Cole-Cole representation (Im vs. Re).

LRT and T 2
Let us now apply the LRT formalism to a problem with a transverse probing field (see also PAPER IV). At negative times the Hamiltonian of our two-level system is then where σ x is the corresponding Pauli matrix (we could use σ y instead). At t = 0 we switch δB x off. As before, for the stationary solutions we will use the canonical form. Then to first order in δB x we have ̺ 11 (0) = ̺ 11 (∞) = e βBz , ̺ 22 (0) = ̺ 22 (∞) = e −βBz , for the diagonal elements, while the coherences (off-diagonals) read ̺ 12 (0) = ̺ 21 (0) = δB x M z /B z and ̺ 12 (∞) = ̺ 21 (∞) = 0. Now we monitorize σ x (the operator coupled to the perturbation) through its average M x ≡ ̺ 12 − ̺ 21 . This relaxes from its equilibrium value ∝ χ eq x = M z /B z as follows From this relaxation curve we can get the AC susceptibility χ x (ω) [Eqs. (6.11) and (6.23)]: This is different from a simple Debye form, and we are going to see that includes absorption. Let us extract the real and imaginary parts These are plotted in Fig. 6.3, where we see again that the equilibrium result is recovered as ω → 0. As for the imaginary part, it consists of two Lorentzian-type curves, centred around the level difference ±B z and having width Γ d . In this case the perturbation did not commute with the base Hamiltonian (6.26), inducing transitions between its levels. Spectroscopists call T 2 (= 1/Γ d ) the life-time of the levels; for us it is the decoherence time (the time to erase the off-diagonals). Letting Γ d → 0 (closed system) the absorption lines go over two Dirac deltas at ±B z . Their broadening in an open system gives T 2 directly. 5

beyond the 2-state model
The above example was a simple case where we could do everything analytically. In general, we will have to deal with quantum master equations like (3.30). But we saw that if H s has a discrete and finite spectrum, we can write a density-matrix equation for the ̺ nm as in (6.14). This is a linear system with N ×N equations. Diagonalizing the system's matrix (numerically most likely), we obtain N × N eigenvalues and eigenvectors. Then, with the appropriate initial conditions we can write where we have split the contribution of the zero eigenvalue (giving the stationary solution). With these ̺ nm (t) we can construct the relaxation curve ∆A(t) for any observable A of interest [Eq. (6.8)]. Plugging it into Eq. (6.11), and doing the integrals with (6.23), we The arrows mark these special points.
obtain the corresponding dynamic susceptibility The a i ∝ n,m c (n,m) i n|A|m involve the matrix elements of A (PAPER III). Real Λ i (the Γ r in the two-level system) contribute Debye-type relaxation terms to the susceptibility. On the other hand, complex eigenvalues (like iΓ d ± B z in the 2-level system) would contribute with absorption profiles.
Let us consider now an observable such that the a i of the complex eigenvalues are all zero (the observable does not see those modes). Then the total susceptibility is made up Debye's; an example with two of them is plotted in Fig. 6.4. Each peak on the imaginary part is located at the corresponding Λ i .
In a multiexponential case like this, one tries to characterize the response with a few effective parameters. The integral relaxation time τ int is defined as the area below the relaxation curve: This quantity governs the low frequency behavior of the susceptibility. On the other hand, the high-frequency properties are better characterized by the initial (t = 0) slope of the relaxation curve. This defines the effective time τ ef Both time constants coincide in a problem with a single exponential decay. We will see examples with several relaxation/absorption peaks in Chapters 7 and 8.

continued-fraction methods
It is not always possible, or practical, to diagonalize our master equations (large number of levels N , continuous spectrum, etc. . . ). In some cases we can resort to the continuedfraction method [14]. This method is a relative of the solution of kinetic equations by expansion into complete sets (à la Grad) [69]. It has been fruitfully exploited to solve classical Fokker-Planck equations [14,70] and quantum master equations [71,72,73,74].
The ideas is to expand the unknown solution ̺ of our dynamical equation ∂ t ̺ = L̺ in some basis {u j } (all this applies to the classical or quantum W ) Plugging this into ∂ t ̺ = L̺ one gets a system of coupled differential equations for the expansion coefficients, that is:Ċ Here I is the range of the coupling. If , the system of differential equations is finally transformed into a set of algebraic equations for the C i , with recurrences between the indexes Compared with (6.35), the coefficient of the central term was simply augmented to (Q i,i − s), while −C i (0) appears as an inhomogeneous term.

solving recurrence relations: case I=1
Let us assume that, with a smart basis choice, we have ended up with a recurrence like (6.36) with I = 1 (first non-trivial case). Then we introduce some simplified notation: f i ≡ C i (0) for the source terms ("forcing"), and the constants Then the recurrence (6.36) would take the "canonical" form: At this point we could diagonalize the matrix associated to the system (6.37). Then, instead of diagonalizing L directly, we would have made an analytical effort leading to 69 a tri-diagonal matrix; the numerical work would be more efficient. But we can take an analytical step further. If we introduce the following ansatz [14], in equation (6.37), we easily get the following relations for the "ladder" coefficients S i and the shifts a j (associated to the inhomogeneous term) Now we can explain the name "continued fraction" given to this method; iterating the denominators with the same formula one could write: This would provide a "closed form" for the coefficients C i [s]. We could say that this is a semi-analytical technique, where we end up evaluating (numerically) continued fractions. Eventually, inverse Laplace transformation would give the time evolution of the coefficients C i , and hence of the full distribution ̺. 6

Brinkman hierarchy: matrix recurrence relations
The best way to learn how to convert a kinetic equation into a suitable recurrence relation is through an example. Let us consider the Klein-Kramers equation (2.16), which in appropriate dimless units reads We follow Risken [14] and begin expanding in a basis [cf. Eq. (6.34)] W (x, p) = w(x, p) c n (x, t)ψ n (p) (6.41) with the pre-factor w extracted for later convenience. Inserting this expansion in the Fokker-Planck equation we get the following "equations of motion" for the coefficients c n : Here L is the differential operator on the right-hand side of the Klein-Kramers equation, whileL is some similarity transformed version of it:L f = w −1 L(w f ).
Following Risken we put w = e −(p 2 /2+V (x))/2 , which results in That is, the dissipative part acts formally as the Schrödinger Hamiltonian of a harmonic oscillator (on the p dependence). This cries for introducing "ladder" differential operators a = ∂ p + p/2 and a + = ∂ p − p/2 obeying [a, a + ] = 1 and giving On the other hand, the closed part (Poisson bracket) gets transformed as with the shifted derivatives D ± = ∂ x ∓ V ′ /2. Finally, if we choose for the ψ n the Hermite functions, ψ n (p) ∝ H n (p/ √ 2), the term a + a becomes diagonal, giving Q nn , while a and a + produce Q nn±1 [Eq. (6.42)]. Putting all this together one gets the so-called Brickmann hierarchy [75,76] having the form of a 3-term recurrence. However, as the D ± still contain operators on x, we expand on a x-basis {u α (x)}, obtaining a recurrence relation, now without operators, but involving two indexes:ċ α n = Q α,α+β n,n+m c α+β n+m . In order to convert this into a one-index recurrence (the problem we know how to solve), we introduce the following vectors and matrices: Here L is some truncation index for the x-basis, and the matrices Q are Finally, we could Laplace transform as in (6.36), getting an algebraic recurrence like (6.37), with the Q's and C's given by the vectors C and matrices Q. The resulting recurrence could be solved like the scalar case of subsection 6.4.1. The only difference is that the fraction bars in (6.40) have to be understood as matrix inversions (A/B → B −1 A). In the statics no Laplace transform is required, and Risken solved the problem studying parameter ranges out of the reach of other techniques [14].

the case I > 1
What happens if the recurrence involves 5 or more terms? Then we build appropriate matrices and vectors (by mounting the ones above) and the recurrence can always be cast into the "canonical" 3-term form [14]. The method will be useful as long as I is not too large. We will solve problems where I = 1 or 2.

stationary solutions, LRT & continued fractions
Let us see how to apply the continued-fraction method to obtain stationary solutions and dynamical susceptibilities. The approach will be used later in chapters 7 and 8.

obtaining stationary solutions
When we want to obtain the stationary solution L̺ = 0 we do not need to resort to Laplace's transform. Indeed, plugging the expansion for the coefficients of the stationary solution. This is an algebraic recurrence relation like (6.37), with f i = 0, solvable simply as discussed in Sect. 6.4.

application to LRT
Let us have another look at the AC experiment of section 6.1.2. Recall that from t = 0 on we apply a sinusoidal probe f (t) = e iωt ; then the evolution operator becomes L = L 0 + δ · e iωt L 1 , while the solution can be written as ̺ = ̺ 0 + δ · e iωt ̺ 1 . Compare with Eqs. (6.3) and (6.4); here we have explicitly extracted the t dependence e iωt from ̺ and L. Therefore, both ̺ 1 and L 1 are time independent (of course, ̺ 0 and L 0 too). Using now̺ = iω̺ 1 e iωt and equating terms with and without oscillating factors we get 0 = L 0 ̺ 0 (6.50) Inserting the expansion ̺ 0 = i C i u i in the first equation we get the recurrence (6.49), which we can solve by continued fractions, getting the C i . Next, we insert the solution for ̺ 0 in the second equation and expand ̺ 1 as in (6.34), finding a recurrence like (6.49), but now with f i = 0. The inhomogeneous part comes from L 1 ̺ 0 , thus justifying the carrying along of the forcing term from the generic recurrence (6.37). In summary, with two recurrence relations we can get ̺ to first order in δ. With this ̺ we obtain the AC response of the averages of interest [Eq. (6.10)].
In this chapter we have considered the system's response to first order in the probing field/force. To go beyond LRT, one would need the response to higher orders [cf. Eq. (6.4)] We will see in chapter 8 how to proceed along this line.
In the second part of the chapter we have briefly introduced the continued-fraction method to solve master equations. This is a semi-analytic technique, and for the cases we are interested in (stationary and first few order responses) we will see that it can be very efficient. Besides, it can be used when the spectrum is continuous (Chap. 7). In addition, one gets the full ̺, not only the observables.
It has some drawbacks too. It is very specific of the problem addressed, requiring a specific pre-analysis in each case (finding a good basis {u j } to expand ̺, obtainment of the required matrix elements, etc.). Besides, it can become numerically unstable in certain parameter ranges (typically very low damping and temperatures [14]). applications I: brownian particle in a periodic potential In this chapter we apply all the preceding open-system formalism to the problem of a quantum particle in a periodic potential V (x) = V (x + L). That is, we study the Bloch problem taking into account the coupling to a dissipative bath. Besides the substrate potential, the system can be acted upon by a force −x · f (t), allowing the study of transport properties (the current). The external force breaks periodicity and the band picture, so that the problem becomes a delicate one (continuous unbounded spectrum, Bloch oscillations, Wannier-Stark ladders [20] . . . ). On top of this the particle interacts with the bath, further complicating matters.
The chapter begins with the solving of the Caldeira-Leggett equation (3.45) in phase space [Eq. (5.17)] by means of continued fractions. In the classical limit this equation, which goes over the Klein-Kramers equation (2.16), was solved in this way by Risken (section 6.4.2). Here we discuss the extension of that work to the quantum regime (PAPER I).
In the paradigmatic example of a cosine potential V (x) = cos(x) we will see how the classical and quantum distributions W (x, p) can differ. After this we address transport problems by adding a driving −x · f (t). We compute the current and the mobility in the cosine potential as well as in potentials lacking space-inversion symmetry: ratchet potentials of the type V (x) = sin(x) + r/2 sin(2x). The phase-space approach will allow us to connect smoothly with the classical results.

solving the Caldeira-Leggett master equation by continued fractions
We start writing the Caldeira-Leggett equation in phase space as [cf. Eq. (5.17)]: The first two terms are simply the classical Poisson bracket. The third one includes the dissipation and diffusion originating from the coupling to the bath. Being a high temperature semiclassical equation (Sect. 3.5.3), those terms are equal to the classical ones. The sum in the last term incorporates the quantum mechanical "corrections" to the closed Hamiltonian dynamics (Moyal terms, section 5.2.2). Equation (7.1) was made dimless by means of a characteristic length x 0 (e.g., the period of the potential) and a characteristic energy E 0 (the potential amplitude). Besides the parameters have been thermally rescaled introducing ω T = (k B T /mx 0 ) 1/2 . Thus the momentum p in (7.1) was made dimless and thermalized (the potential V → V /k B T too), while the other parameters read Here ω 0 and S 0 are the characteristic frequency and action [ω 0 = (E 0 /mx 2 0 ) 1/2 and S 0 = E 0 /ω 0 ]. The temperature has been absorbed in the definition of γ T and the Moyal coefficients κ (s) . In this way the form (7.2) shows that the quantum terms become less important the higher K and/or the temperature are. Indeed the parameter K tells how quantum the system is (K ∝ 1/ , Chap. 2, Sect. 2.6). Thus the limit K → ∞ recovers the Klein-Kramers equation from the Caldeira-Leggett master equation in phase space. From now on we will simplify further the notation by omitting the subindex T en γ. 1

the classical limit (K → ∞)
As we stated in chapter 5, one advantage of the phase-space formalism is that by varying one parameter the master equations go over their classical Fokker-Planck counterparts. In our problem, letting K → ∞ in (7.1) we recover the Klein-Kramers equation (2.16), in dimensionless form In section 6.4.2 we transformed this equation into a recurrence form, the so called Brickmann hierarchy (6.46). Recall that what we did was to expand the distribution in Hermite polynomials for p while the x-basis was left unspecified, W (t, x, p) ∝ c α n (t)u α (x)ψ n (p). This led to a 2-index recurrenceċ α n = Q α,α+β n,n+m c α+β n+m which was converted into a recurrence in one index by defining the vectors (6.47) and the matrix coefficients (6.48) In this form, the recurrence could be solved by continued fractions (Sect. 6.4). The xbasis choice is usually dictated by the space symmetries of the problem including boundary conditions.
1 Note the relation between K and the thermal de Broglie wavelength. The definition λ dB = ∆x = /∆p with ∆p = √ M kBT (thermal uncertainty from equipartition) gives

application to periodic potentials
If the substrate potential is periodic, one Fourier expands as V ′ (x) = q V ′ q e iqx . Then it is natural to choose plane waves for the x basis To get matrix elements A αβ = dx u α A u β , we use orthogonality From these we can get the matrix elements of Q, which have the structure Q ± α,α+β ∝ [iα δ α,α+β ± V ′ β ]. Hence, the number of diagonals with non-vanishing elements equals the number of harmonics of V (only one for the cosine potential).
Hitherto we always had in mind solving the recurrence (7.5) in the index n. But we could also have chosen to transform the 2-index recurrenceċ α n = Q α,α+β n,n+m c α+β n+m by introducing the following vectors In this case I equals the number of harmonics in V , since Q α,α+β = 0 for β larger than that number. Both recurrences (in p or in x) can be employed at convenience as Risken already did in the classical case [14].

example: the cosine potential
Let us consider the stationary solutions W (t → ∞; x, p) for a particle in a cosine potential, tilted by a constant force (see fig. 7.1): To this end we setĊ n = 0 (orĊ α = 0) and solve the vector recurrences. With the c's so obtained we reconstruct W (x, p) ∝ c α n u α (x)ψ n (p) [cf. Eq. (6.41)], which is plotted in figure 7.2 (still classical results).
Let us follow Risken in the interpretation of these graphs. In the absence of applied forces (F = 0) the stationary solution is the equilibrium Boltzmann distribution ∝ e −p 2 /2−V (x) . In the opposite case, when the force is very large, the substrate potential matters only a little, and we can think of a damped particle dragged by a constant force. In this limit p = F/γ. Therefore, from small to large forces the distribution shifts in p, from the Boltzmann case with p = 0 towards the stationary solution ∝ e − 1 2 (p−F/γ) 2 . In other words, a transition from "locked" to "running" solutions as a function of the force (with bistability in between). In the Hamiltonian case (γ = 0) the particle can be "locked", oscillating in one well, or "running" with p = 0. These two cases are depicted by solid trajectories. When coupling the particle to a bath it experiences fluctuations and dissipation. These can lead to the eventual trapping of a running solution into one well, or to the thermal launching of a locked particle over the potential barrier (dashed line) possibly into a running solution. Risken's graphs: classical distributions for a particle in the tilted periodic potential (7.9). The distribution is periodic in x so we have just plotted two periods for illustration. The damping is γ = 0.05 and the forces F = 0, 0.075, 0.15, and 0.2 (left to right, top to bottom). Note the p = F/γ ∼ 4 for the largest force, and the bistability between running and locked solutions for intermediate forces.

the quantum regime (K < ∞)
When K is finite, we have to take into account the rest of the terms in the Moyal series in equation (7.1). Here is when we start to appreciate the advantages of the phase-space formalism: one part of the equation is already solved -Risken did it for us! Now we are left "only" with handling the last sum, i.e., ∞ p). To this end, we make use of the result where recall that a = ∂ p + p/2 and a + = ∂ p − p/2. Then, using the Hermite p-basis each term in the sum will couple n with n ± (2s + 1). That is [77], the Moyal term will in principle give coupling to all indexes in the quantum generalization of the Brickmann hierarchy (6.46)Ċ n = m Q n,n+m C n+m .
The explicit form of the matrices Q can be found in PAPER I, sections 4 & 5. It is important to remark that for polynomial potentials V (2s+1) = 0, after some s. Then the recurrence would be finite, providing the finite I we required in the previous chapter to solve by continued fractions. On the other hand, for non-polynomial potentials there exist non-vanishing derivatives at all orders, like in a cosine potential. We are back in the generic case of the quantum Brickmann hierarchy (7.11) with infinite coupling range, making continued fractions of no use. But we still have the x-recurrence, let us see if we are more lucky with it . . .

the case of periodic potentials
In this case we proceed expanding the x dependence in plane waves (7.6). The "new part" (the Moyal sum in powers of ) brings in the derivatives V (s) = q V (s) q e iqx . To get their contribution to the matrix elements A αβ = dx u α A u β , we use the following result in the plane-wave basis: Therefore we have Q α,α+β n,n+m = 0 whenever β is larger than the number of harmonics in the potential. Then, using the vector recurrence (7.8) in the α index, we have a finite coupling range, just as in the classical case. The matrices Q α,α+β will contain more stuff, -dependent terms, but what is important is that the index I is finite and the same as in the classical problem: equal to the number of harmonics. In this case we could not chose between x or p recurrence (of no use), but luckily one of them works.

example: cosine potential
As we did in the classical case, let us compute the stationary solutions of the master equation in a cosine potential. This potential has one harmonic (I = 1) yielding a 3-term recurrence relation. SettingĊ α = 0 and solving, we get the distributions of Fig. 7.3. Here we have set a "very" quantum regime, with K = 1, 2 (recall K ∝ S 0 / ). In order not to leave the validity range γ/k B T ≪ 1 of the Caldeira-Leggett master equation (Sect. 3.5.3), we have used dimless γ = 10 −6 and T = 1 (ensuring γ/T K ≪ 1). The first striking feature (already announced in chapter 5)) is the distribution becoming negative in some regions (the white "islands" in figure 7.3). Besides, compared with the classical graphs 7.2, the distribution is quite delocalized in x. Specially for K = 1; for K = 2 we start to recognize the first symptoms of localization.
To interpret with the hands the distributions of figure 7.3 we recall first that in the range γ/T = 10 −6 ≪ 1 the stationary solutions should be well approximated by canonical distributions (chapter 4). Recalling also how we build the Wigner function from the density matrix [Eq. (5.4)], the equilibrium W eq would be: dy e iyp/ ̺ eq (x − 1 2 y, x + 1 2 y) ; ̺ eq = dq e −βǫq |q q| (7.13) This will be approached by our t → ∞ stationary solutions. The Bloch waves |q , for small K, can be constructed perturbatively ( [78], section "potential treated as a perturbation") (7.14) Thus, when the characteristic action is "small" compared with (small K/2π), the Bloch particle behaves almost as a free particle (see the energy bands in figure 7.4). Inserting eventually the above |q in the equilibrium (7.13) we obtain: This remarkable form captures qualitatively the features of the numerically obtained distributions. The first term (a Gaussian in momentum) "is" the free particle e iqx → exp(−p 2 /2|η| 2 ), while the second comes from shifted plane waves e i(q±1)x . They are responsible for the delocalization along the x axis. The last two terms in (7.15) follow from the cross products between e i(q+1)x & e i(q−1)x and e iqx & e i(q±1)x , respectively. Specifically, the third term incorporates a modulation due to the potential profile (as it is of order |η| 2 is less noticeable for K = 1). The fourth term, finally, produces negative regions. The above terms, and hence the negative islands, can be seen as an expression of the wave character of quantum mechanics. By this we mean that W generates the averages and hence it is made of squared wave functions [ψ * (x − u/2)ψ(x + u/2), section 5.2.1]. Seen in this way, W is like an intensity (its Fourier transform) and the last terms in (7.15) depict a typical interference pattern between sinusoidal waves.

computation of observables: transport problems
We did not mention it in this chapter, but we recall that once the distribution W (x, p) is reconstructed from the expansion coefficients c's, we can compute any observable following essentially the classical prescription [Eq. (5.
The "classical observable" A(x, p) is actually the W A transform of the operator A in Hilbert space, Eq. (5.2).
In some cases, it is not required to construct W explicitly, as some observables of interest can be expressed directly in terms of expansion coefficients. For example, in transport problems the current p can be written as The coefficients I α and K (1) 2n+1 depend on the prefactor w in Eq. (6.41), and on properties of the integrals of Hermite polynomials and plane waves [given explicitly in Eqs. (29) and (30) of PAPER I].

example I: quantum transport in a cosine potential
Let us focus now on transport problems, and the effects of the dissipative bath. We consider first that the substrate potential is the simplest possible, the cosine, which is however tilted by a constant force −x · F (see figure 7.1). For this problem we solve the Caldeira-Leggett equation and compute the current in stationary conditions p (t → ∞). Figure 7.4 shows the results for several values of the damping. The phenomenology is that studied by Chen and Lebowitz, which did a difficult path-integral calculation treating the potential perturbatively [79,80]. The curves can be understood as follows. When the force is small the particle moves on the band bottom, which for K = 1 has a near free particle form ( figure 7.4, right). Then p = F/γ as in a damped and dragged classical free particle. Increasing the force, we approach the band limit. There Bragg scattering reduces the group velocity ∝ ∂E/∂q, and hence p decreases (see any solid state book, e.g. [81]). But increasing further the force, we go into the next band, progressively recovering the free particle behavior p = F/γ. Figure 7.4 also shows how increasing the damping the curves approach this free particle current. This is understood in terms of the folk argument of bath as blurring the energy levels (recall figure 4.2), which in our case amounts to bridge the band gap, responsible for the decrease.

example II: directional motion in a ratchet potential
Up to this point we have treated the cosine potential. In this section we consider the minimal extension of adding a second harmonic with r the ratchet parameter (r = 0 brings us back to the one-harmonic case). This potential is asymmetric, as seen in Figs. 2.3 and 7.5. Looking at the potential profile we realize that the slope to the right is less steep (easy side) than to the left (hard side). An important property of ratchet potentials is that they can accommodate directional motion. By this we mean that the system can have a net current even if the average of the external forces is zero (on average, the force does not do work). One also speaks of rectifying zero-mean perturbations (generating a net response). How is this possible? First, if the system is out of equilibrium, the second law of thermodynamics does not hold. Besides, we need the potential being asymmetric, as otherwise there would not  exist a privileged direction [82,83]. Then a system modelized by an asymmetric potential like (7.17), and coupled to a dissipative bath, can show such a phenomenology if a driving force x · f (t) keeps it out of equilibrium.
For the numerical solution we are going to assume, for simplicity, that f (t) has a square wave profile, which changes between positive and negative values ±F with a period 2π/ω (see Fig. 7.6). Besides, we use the adiabatic trick: the period is long enough (compared with the rest of time scales of the problem), so that most of the time inside each step +F or −F the current is the stationary one at the corresponding force. Then the net current (averaged over time) would be: with p ±F the corresponding stationary solutions. 2 Working in the stationary regime, we setĊ α = 0 in the recurrence (6.47) and we are back in the case discussed in section 6.5.1. As for the continued-fraction solution, note that having two harmonics in the potential, the associated recurrences now have I = 2 as coupling range (section 6.4.3).  [14] in the overdamped regime γ → ∞. The vertical lines mark the two critical forces (7.19).

directed motion in the classical limit
Let us try to understand first the classical problem in order to address later the effects of adding quantum effects. The classical rectified current in shown in figure 7.7. The main feature is that the net current is positive, meaning that the particle moves on average to the right, the easy way. A way of understanding this is introducing the critical forces for barrier disappearance, to the left and to the right: Thus, the barrier felt by the particle when tilting the potential with +F is lower than that when we tilt with −F , so it is easier to move to the right. To understand the other features, we will think first on the T = 0 limit; that is, only dissipation, no fluctuations. The associated deterministic system has two more characteristic forces F 1 and F 2 , obeying F 1 < F 2 < F c . When F < F 1 the dynamically stable solution (the attractor) is the locked solution, while for F > F 2 the running solution is the stable one. Well, if the system is locked close to a minimum (F < F 1 ), or if it is in a wild running mode, in both cases the potential asymmetry (and other details) become less relevant. Therefore, we understand that the rectification will be more noticeable in some intermediate force range, some "window", as seen in figure 7.7. Besides, the characteristic forces F 1 and F 2 are damping dependent, and happen to decrease with γ [14]. Thus the window of maximum rectification shifts down along the force axis when we decrease γ [84]. Finally, to restore T in the picture, we can think of its effect as allowing jumping between attractors, smoothing the sharp ranges, and in particular allowing rectification below F + c (PAPER I and PAPER II).

quantum corrections
Our goal here is to study how the rectified currents of the previous section are modified when the quantum properties of the system are taken into account. To this end we focus on one of the curves (we fix one γ) and progressively decrease the quantumness parameter K ∝ S 0 / . This makes quantum effects more and more important. In figure 7.8 we see that for large amplitudes of the force wave |F | the quantum corrections increase the rectified current p +F + p −F , whereas the rectification is decreased at low forces with respect to the classical case.
To understand the reduction/increase, in PAPER I and PAPER II we resorted to the semiclassical calculation of the transmission through a saw-tooth potential. This is an asymmetric potential, much as the ratchet potential, with easy and hard sides (figure 7.9), but where analytical calculations are possible. We studied the transmission coefficient when deforming the potential with +F and −F (appendix F, PAPER I). We found that for energies over the barrier the transmission is higher to the easy side, because the wave-reflection is less intense than to the hard side (schematically indicated by the broken arrows in figure 7.9). On the other hand, when the particle has an energy lower than the barrier, the transmission is higher to the hard side, as tunneling is more probable towards that direction (solid arrows in figure 7.9). 3 Then, for large enough forces, particles are more energetic and there is more probability of being launched over the barrier by thermal fluctuations; then the wave-reflection dominates and we see an amplification of the rectification by quantum effects. On the Figure 7.9: Saw-tooth potential in the absence of force, and as deformed by ±F . The barrier is lower to the easy side than to the hard side. The solid arrows indicate the transmission for energies smaller than the barrier (tunnelling) and the broken arrows the transmission for energies larger than the barrier (affected by overbarrier wave reflection). The length of the arrows indicates qualitatively the probability of those processes. other hand, for small forces/energies, it is more probable to find particles with energies lower than the barrier. Here tunnel events win, and the rectified current is reduced.
In PAPER I we presented results for fixed amplitude F , studying the dependence of the current on the temperature.

example III: AC response in periodic potentials (LRT)
In the previous examples we have addressed stationary transport in a cosine potential, subjected to a constant force, and then rectified currents in a ratchet potential, with a driving force in the adiabatic regime. Thus, we could reduce the studies to the computation of stationary solutions, and combinations of them. Now we are going to address a genuine AC problem, with the driving out of the adiabatic regime.
Let us return to the cosine potential (7.9), and let f (t) = δ · e iωt . If the amplitude is small δ ≪ 1, we can use linear response as in section 6.1. The observable of interest now is the AC current [cf. Eq. (6.10)] (7.20) To compute the amplitude χ(ω) we can use the continued-fraction method as explained in section 6.5.2. Note that, being in linear response, p is proportional to the force, so that f (t) = 0 gives a zero time average of the current, irrespective the potential. Then, the asymmetry of the potential does not bring much new, and for simplicity we return to the cosine potential. Figure 7.10 shows details of the response in the low frequency range (left panel) and at high frequencies (right). The low frequency spectrum shows a Debye profile, while at high frequencies we see a pattern of absorption lines (we just plotted the imaginary parts). These figures can be understood qualitatively recalling the discussion of section 6.3. There The reason of the Debye part is the following. Along with the erasing of the nondiagonal terms in the density matrix (decoherence), the system has to readjust the populations of the wells (the cosine potential is now tilted by the AC force). This leads to the relaxation time (slow dynamics) we find in the low frequency range.
As for the high frequency curves, the widths of the absorption peaks can be calculated semiclassically (PAPER I, appendix F). We performed a WKB calculation of the widths W b of the lower bands for the example displayed (K = 200), finding W b = (4/K) exp[−K(1 − ǫ)/2] with ǫ = E/V 0 . Therefore the widths are exponentially small, and for K = 200 they can be considered as discrete levels (recall Fig. 2.3). Besides, their position can be estimated from the energy levels of the anharmonic oscillator obtained expanding the cosine to fourth order cos(x) ∼ = ax 2 + bx 4 . Indeed, plotting the associated energy differences, we find that the absorption peaks are centered around them, as expected. Now, if we increase the damping, we see how the absorption curves broaden, with a width proportional to γ ∝ 1/T 2 . Eventually, as the damping increasingly blurs the energy levels, the curves approach a smooth peakless "classical" profile.

summary
The example of the particle in a periodic potential has allowed us to review many of the concepts and tools introduced in the preceding chapters: Fokker-Planck equations, quantum master equations, phase-space methods, relaxation and decoherence times, . . .
The continued-fraction technique, plus phase space, has made possible to solve the problem as an extension of Risken approach to the classical problem. In this way we have bypassed the problem of the underlying spectrum (demanding, as is continuous and non-bounded). During the whole chapter we had the classical limit as reference, in order to understand our results as corrections (recall that the Caldeira-Leggett equation is semiclassical).
We finally remark the differences between the classical and quantum distributions (recall the wave-mechanical interpretation). We have discussed how the interplay of overbarrier wave reflection and tunneling can increase or decrease the rectified current in an asymmetric potential, depending on the force/energy range. We have also seen how the damping can make the curves to approach the classical counterparts. We closed the chapter with the study of the dynamical susceptibility spectra and characteristic times. 88 Chapter 8 applications II: statics & dynamics of superparamagnets Quantum nanomagnets are an example of mesoscopic objects described effectively by "one" degree of freedom (their net spin). In this chapter we are going to study the statics and dynamics of these systems, where the interaction with the environment is important. We start reintroducing the model and investigating its equilibrium properties. Next we specify the master equation to be solved to study non-equilibrium behavior. As examples we will address the relaxation mechanisms (analogous to the calculation of T 1 with Bloch's equations), by analyzing the longitudinal susceptibilities (linear and non-linear). We will conclude with the calculation of the transverse susceptibility (study of T 2 ).

spin Hamiltonian
Superparamagnets are solids or clusters of nanometric size where each molecule can be described by a spin Hamiltonian [21]: Here D is the anisotropy parameter, originating from the spin-orbit coupling, while −B z S z is the Zeeman term accounting for the coupling to external fields. As a first approximation the magnetic molecules can be considered as independent of one another, so the total Hamiltonian is a sum of contributions like (8.1). Its eigenvalues are ǫ n = −Dm 2 − B z m. 1 For D ≥ 0 the anisotropy term in (8.1) yields a bistable structure with minima at m = ±S, while the Zeeman term tilts this "double well" (figure 2.4, chapter 2). In equilibrium the spin orientations are like those of a paramagnet (the magnetization and susceptibility are given by Brillouin and Curie type laws), but due to the large S they are called superparamagnets. Physical realizations of this model are the Mn 12 or Fe 8 clusters, among others [22], which have a moderately large spin S ∼ 10, and where quantum effects of the spectrum's discreteness can be studied [4].
It is convenient to introduce scaled variables (β = 1/T ): The parameter h is called the reduced field. When there is no external field obviously h = 0, whereas h = 1 just when the barrier disappears. In these variables, the thermal Hamiltonian reads βH s = −σ(m/S) 2 − ξ(m/S) and we will take the classical limit S → ∞ keeping σ and ξ fixed (note |m/S| ≤ 1). In this way the energies are kept constant, while introducing more and more levels towards a continuum. Thus the spin number S plays the same role as K in the particle problems (characteristic action over ). In that case when K → ∞ the number of bands below the barrier was increased, but the potential height was kept fixed. From a formal point of view, varying parameters allows the study of problems with equispaced spectrum (D = 0), or problems with more featured spectra (D = 0). In the case D > 0, we have mentioned that the spectrum has a double-well structure, while D < 0 has a single well (figure 2.4). The level spacings are so that, when B z = 0, the levels m and −m become degenerate.

Hamiltonian coupling to the environment
As these spins "live" in a crystal, they interact with the lattice phonons, providing an example where the coupling to a bath of oscillators is well justified. Thus we write [cf. (3.26)]: where H s is the Hamiltonian accounting for renormalizations (recall sect. 3.4). As we see, these systems can be studied within the open-system formalism of the previous chapters. The coupling between the spin and the phonons includes where the complex conjugate constants η − = η * + ensure hermiticity, the anticommutator is represented as [ , ] + , and v(S z ) is any function of the operator S z . Note that the coupling is realized through the ladder operators S ± = S x ±iS y , which obey S ± |m = ℓ m,m±1 |m±1 with the custom factors ℓ m,m±1 = S(S + 1) − m(m ± 1).
In the case v(S z ) = cte, the coupling F ∝ η + S − +η − S + is linear in the system variables. This is analogous to the coupling used in the particle chapter, and here will serve us as 8.3. EQUILIBRIUM PROPERTIES control model. On the other hand, v(S z ) ∝ S z for spin-phonon (magneto-elastic) coupling [85,86], and we have to use the more realistic F ∼ [S z , S ± ] + . In the former, linear coupling case we use an Ohmic distribution of bath oscillators, J(ω) = λ ω, as we did previously. But when studying the magneto-elastic model we will use a super-Ohmic bath with a spectral distribution J(ω) = λ ω 3 appropriate for phonons (or photons) in three dimensions [1].
Both in the equilibrium properties and in the dynamics we will need the explicit form for the matrix elements F nm ≡ n|F |m : Note that when v(m) = cte (spin-phonon) the L m,m ′ add an extra dependence on m in the coupling/damping, compared with the linear case; if you wish, the spin analogue of position-dependent damping in particle problems [14].

equilibrium properties
In chapter 4 we discussed the equilibrium properties of open systems. In the previous chapter we already showed stationary Wigner distributions for a particle in a periodic cosine potential. There, as we worked with γ/T ≪ 1, the stationary solution was essentially the canonical distribution at zero damping (figure 7.3). Now we want to study explicitly corrections to the thermodynamical properties when λ = 0 in spin problems. Let us recall that to obtain the damping-dependent corrections we had to compute certain Z (2) . Then Z tot /Z b ∼ = Z s + Z (2) with Z tot , Z b and Z s the partition functions of the whole system, the bath and the closed system, respectively, Eq. (4.9). From Z we can obtain the free energy, and hence all other thermodynamical quantities, as in section 4.3.

example 0: equilibrium of the isotropic spin
If we consider D = 0 (isotropic spin) the spectrum becomes equispaced and the arguments of the Φ in (8.7) are simply ±B z (the constant level difference). As these do not depend on the index m, they can be taken out of the sum, as in the example of the harmonic oscillator in section 4.4. Then Z (2) can be explicitly calculated [cf. Eq. (4.20)]:

dynamics
It will not surprise the reader if to study the non-equilibrium properties we use the quantum master equations of chapter 3. Specifically, we will use the second order master equation (3.41), in what we called the improved rotating-wave approximation (RWA) d̺ nm /dt = −i∆ nm ̺ nm + R nn+1,mm+1 ̺ n+1m+1 + R nn,mm ̺ nm + R nn−1,mm−1 ̺ n−1m−1 .
(8.11) At the end of this chapter we will briefly touch the comparison between crude RWA, and its improved version.
At some point we would like to add a small transverse field, e.g., B x , to study the transverse response χ x (ω). Then the system Hamiltonian H s becomes with B ± = B x ± iB y . But for arbitrary B ± we do not know the eigenvalues and eigenvectors of H s analytically, so we cannot calculate the unperturbed evolution of the coupling term F nm (τ ) = F nm e −i∆nmτ (which was required to calculate the relaxation term R in section 3.5.1). However, restricting ourselves to B ± ≪ 1 we can approximate the required evolution by the diagonal part. Naturally, we can and will fully keep the transverse field in the unitary evolution. All things considered, we will have to solve where ∆ nm are the energy differences for the longitudinal part (8.3). For the explicit form of the R's, see PAPER IV. A final comment on the resolution method. Note that the master equation (8.13), can be formally written as̺ nm = α,β Q mm+β nn+α ̺ n+α,m+β . Then, introducing the matrices [Q n,n+α ] m,m+β = Q m,m+β n,n+α and the vectors [C n ] m = ̺ nm , we havė C n = Q n,n−1 C n−1 + Q n,n C n + Q n,n+1 C n+1 . (8.14) That is, we have rewritten (8.13) in a form suitable for the use of continued fractions (chapter 6, PAPER IV).

exampleI: longitudinal relaxation & linear susceptibility
Our first example of dynamics is a study of the relaxation mechanisms by means of the linear susceptibility. We will consider that no transverse fields are applied (B ± = 0), and analyze the time evolution of M z = S z : In this problem we just need the diagonal density-matrix elements to construct the response M z . Besides, for B ± = 0 the diagonal and off-diagonal elements in (8.13) get decoupled. Then it suffices to solve the system of 2S + 1 coupled differential equations for the populations N m (balance equations): where the transition probability is P m|m ′ = R mm ′ mm ′ . 3 We will calculate the linear response of M z to an ac perturbation δB z e −iωt . To this end, recall that χ z (ω) could be constructed as in Eq. (6.31), using the eigenvalues and eigenvectors of the balance equation without perturbation; and that equation (6.11) relates the response in the time and frequency domains.
We can compute χ z (ω) directly with continued fractions too, as in section 6.5.2. However, the eigenroute not only gives χ z (ω), but allows to analyze and characterize its contributions (while continued fractions are a kind of a black box). The remainder eigenvalues are positive (ensuring convergence to the stationary solution as t → ∞) and real (the density matrix is Hermitian). Figure 8.2 also shows that the eigenvalues are organized in two "groups". Far away from the others we find Λ 1 , which generates the slow long-time dynamics (τ 1 = Λ −1 1 ). Indeed the structure of the eigenmode |Λ 1 (figure 8.2, 3rd panel), with population decrease in one side and increase in the other, corresponds to transfer from one well to the other. As tunnel is not present (no transverse field), this population transfer should be channelled over the barrier top, and hence we call it the overbarrier process. On the other hand, the rest of eigenvalues give population readjustments inside each well (fast dynamics); we call these the intrawell processes. The former, Λ 1 , depends exponentially on the barrier height, whereas the intrawell modes depend only polynomially. Whence the large separation of the scales of these two processes seen in figure 8.2 (note the logarithmic scale). Well, having the eigenvalues and eigenvectors, 4 we can construct now the dynamical susceptibility. The eigenvalues, being real, would produce a susceptibility in the form of a sum of Debye profiles (section 6.3)

eigenvalues and eigenvectors of the relaxation matrix
We have plotted a few susceptibility curves in figure 8.3. The response is given, at most, by two main contributions. We associate the low frequency peak with the slow overbarrier process. While the peak at high frequency, when present, is to be attributed to the fast intrawell dynamics.

bimodal approximation (analytical)
To analyze the susceptibility curves, given the mode's structure above, we tried a relaxation profile consisting of two effective modes (as Kalmykov et al. did in the classical case [89]) Here χ eq is the equilibrium susceptibility and τ 1 = Λ −1 1 , which follows an Arrhenius-Kramers law: This exponential dependence on the barrier height ∆U is typical of problems of escape out of a metastable minimum [90,91]. The specific calculation for our problem is given in PAPER III.
The susceptibility above has two more parameters, a 1 and τ w , which can be determined following the idea proposed by Kalmykov et al. in the classical problem [89]. They argued that we have two parameters to determine, while χ z (ω) can be written in closed form in both the limits of low and high frequencies; in terms of τ int for ω → 0 [τ int is the area below the relaxation curve, Eq. (6.32)] and in terms of τ ef at high ω [the initial slope of the decay (6.33)]. The advantage is that both τ int and τ ef can be computed analytically. Then, these two limits would provide the two unknown parameters a 1 and τ w , so yielding a closed formula without fitting parameters (PAPER III). The results of the analytical ansatz are the lines in figure 8.3, which account quite well for the symbols, obtained constructing χ z (ω) from the numerical eigenstuff.
Let us comment on the curves. The relative importance of the two effective modes is governed by the parameter a 1 : when a 1 = 1 there is only overbarrier response, while switching to a 1 = 0, the whole weight is shifted to the intrawell response. Figure 8.4 shows how this is realized by changing the field B z , with a 1 = 1 at zero field and a 1 decreasing towards zero at large fields. Indeed, at large h there is no barrier, so the overbarrier contribution is naturally expected to vanish. However, figure 8.4 shows that a 1 ∼ = 0 much before the barrier has disappeared. We close with the isotropic spin. If we set D = 0 , we always have a single peak (see [92] and PAPER IV). In a sense, this is clear now for us; in the isotropic spin there is no barrier imparting a large time scale separation, and this kind of relaxation measurement does not resolve the differences between the fast modes.

example II: non-linear susceptibility
Let us illustrate, with the example of the longitudinal relaxation, how to proceed if we want to go beyond the linear response regime. In this and the next section we will see how going to higher orders in the response, one can obtain more information on the system [93,94]. In the nonlinear case, we do not have anymore an equation like (6.11) relating the behaviour in the time and frequency domains. Here we are going to study directly the response as a function of the frequency of the driving field.
To get the dynamical susceptibilities, we excite our system with δB z e iωt and proceed as we did in section 6.5.2. That is, we start from the master equation written formally as: and we Fourier expand the evolution operator L = L 0 + δ · e iωt L 1 + δ 2 · e i2ωt L 2 + · · · as well as the solution ̺ = ̺ 0 + δ · e iωt ̺ 1 + δ 2 · e i2ωt ̺ 2 + . . . Here we have explicitly extracted the temporal dependencies e kiωt from ̺ and L. Therefore, ̺ k and L k (k = 1, 2, . . . ) are time independent (of course, ̺ 0 and L 0 too). Now, equating order by order we have (we have written to second order only, but the perturbative structure is clear). Once again, this chain of equations can be solved by continued fractions, as in 6.5.2. First we solve for ̺ 0 ; this is then inserted in the second equation, and we solve for ̺ 1 ; knowing ̺ 0 and ̺ 1 we can solve for ̺ 2 in the third one.
In some cases this set can be solved analytically. For example, for spin 1 2 we had the equation of motion for the magnetization given by Bloch equations (section 6.2) Then L = Γ r (a number) and operating like in the chain (8.21)-(8.23) we find simple expressions. To first order we recover the Debye (6.24), while for the first nonlinear susceptibility we have Here we have written τ 1 for the relaxation time, while χ eq k = ∂ k M z /∂B k z , k = 1, 2 are the equilibrium linear and nonlinear susceptibilities. As we see, the nonlinear response gives information, not only of τ 1 , but of its B z -derivative as well, τ ′ 1 . As mentioned above χ 2 (ω) seems more sensitive to some parameters of the problem.
Well, this was for S = 1 2 , but we want to get the response for any S. In general, the response will include relaxation "blocks" of the type (8.25), and maybe cross-terms. As we did for the linear susceptibility, we will try to find an approximate analytical modelization.

generalization of the bimodal approximation to nonlinear response
The reader could rightly guess that we are going to try a similar reduction of the 2S modes to 2 effective ones (intrawell and overbarrier) to account for the nonlinear response. Thus, we construct an effective model with dynamics governed by the two relevant relaxation times of our problem: τ 1 and τ w . Although now we need to use 2 × 2 matrices in (8.21)-(8.23), one readily gets: where This is a generalized bimodal formula for the nonlinear response, without any free parameter, as before. The new coefficient, a 2 , tells when the intrawell modes enter the scene (decreasing from 1), and was already plotted in figure 8.4. We see that for a given field a 2 > a 1 , so that the fast modes cannot be detected at lower fields using χ 2 (ω). Had we obtained a 2 < a 1 , we would have used χ 2 (ω) to detect the intrawell modes without needing as large fields as when using χ 1 .
In figure 8.5 we compare the new formulas with the numerical results obtained by solving the chain (8.21)-(8.23) with continued fractions. We see how the bimodal ansatz works very well in the nonlinear case too. We also see how the intrawell modes as less noticeable in the nonlinear susceptibility χ 2 (ω), as we anticipated from a 2 (B z ).
8.8 example III: τ 1 and tunneling -experiments in Mn 12 In this section we will use the nonlinear formalism for the analysis of experiments in the superparamagnetic molecular cluster Mn 12 . In particular, we will see how to get more information on the relaxation time τ 1 from the χ 2 (ω) vs. ω spectra.
Well, in actual superparamagnets the Hamiltonian is a bit more complicated than (8.1) [4]: The term −E(S 2 x − S 2 y ) is a consequence of the distortion of the anisotropy axis, while B x accounts for possible transverse fields. But when these terms are small they can be treated perturbatively.
Let us first discuss (with the hands) how these extra terms, which do not commute with S z , affect the relaxation time τ 1 (8.19). In the absence of transverse terms, the states |m and |−m are degenerated at B z = 0. Now, the terms E and B x mix |m and |−m with ǫ |a − ǫ |s = Ω the tunnel splitting, which can be calculated perturbatively [86,95]. These states are delocalized (in m), and the spin can "tunnel".
If we now apply a field B z the degeneration is lifted and we block the tunnel channel. This can be seen from the perturbed eigenstates: The transverse terms are assumed small (Ω ≪ 1), so we have α + ∼ = 1 and α − ∼ = 0 at not very large fields B z . Therefore, out of resonance we have |s ∼ = |m and |a ∼ = | − m , and we are back in the behavior with no transverse terms. However, increasing the field up to B z = D, the levels |−m and |m − 1 enter in resonance, making possible tunneling through these states. This is repeated at all crossing fields, B z = kD, where the states |−m and |m − k become degenerate ( figure 2.4). Therefore, by changing the field B z , we "open" and "close" tunnel channels.
The tunnel splitting decreases exponentially from the states near the barrier top to the lower laying states near the bottom of the wells. This can be understood because the terms S 2 ± and S x couple states with m ± 2 and m ± 1 in each perturbative order: then, to connect S all the way to −S we need to go to higher orders in the perturbation. 5 On the other hand, in the temperature ranges we are going to address (T > 2 K) tunnel proceeds incoherently. That is, the decoherence time is so short, that coherent oscillations, left-right, are never observed [93]. In this incoherent regime, the relaxation time τ 1 can be written as [4]: where the sum extends to the resonant pairs m and m ′ . As we said, for the lowest laying pairs τ m,m ′ ∼ = 0. As a matter of fact, the barrier that the spin has to overcome decreases in the crossing fields (as the tunnel channels get open) and we have τ 1 ∝ e β∆U ef with ∆U ef < ∆U . In figure 8.7 we have plotted the formula (8.31) for a simple model with transverse field and displayed experimental data for Mn 12 . At the crossing fields, where the levels enter into resonance, the relaxation time shows dips, as a consequence of the decrease of the effective barrier [95]. How could we see this behavior studying χ 2 (ω) curves? Well, if we work with not too strong static fields, we have a 1 ∼ = a 2 ∼ = 1. Besides, if we probe the system with frequencies of the order of τ −1 1 , the fast intrawell contribution to (8.26) is negligible (ωτ w ≪ 1, since τ w ≪ τ 1 ). Then we can work with an approximate χ 2 (ω) dominated by the overbarrier contribution At fields B z close to the crossing fields, the derivative τ ′ 1 becomes large (we are inside the dips). Therefore, we expect that the term with τ ′ 1 would dominate in χ 2 (ω) in those ranges. On the other hand, the derivative changes sign, and so would do χ 2 (ω) at the crossing fields.  We have plotted χ 2 (ω) as a function of B z in figure 8.8. The experimental data shows the same qualitative features as those just discussed on the basis of the formula (8.32). We thus see that χ 2 (ω) is a good tool to detect tunnel in these systems, as it provides resonant features (peaks and sign changes) not easily mistaken. And all this because χ 2 (ω) is more sensitive to the variations in τ 1 that the linear susceptibility (PAPER V).

example IV: transverse response
The linear susceptibility is still useful, as it highlights some other aspects of the problem. In this last example we are going to compute the transverse susceptibility, as we already did for the spin 1 2 case. To this end, we excite with a transverse perturbation δ · B x e −iωt and we monitor the response of M x = S x M x = 1 2 S + + S − = 1 2 m (ℓ m−1,m ̺ m−1,m + ℓ m+1,m ̺ m+1,m ) (8.32) In this problem, the unitary part in the master equation (8.13) mixes off-diagonal and diagonal elements, incorporating the ladder factors ℓ m,m±1 = S(S + 1) − m(m ± 1). Therefore, we have to consider all the density-matrix elements. This can be done solving with the continued fraction method (section 6.5.2). The response obtained is displayed in figure 8.9.
Based on the results of the spin-1/2 in 6.2.2, one could expect some absorption spectrum. The imaginary part is indeed made of Lorentzians centered at the level differences. But due to the anisotropy, the H s spectrum is not equispaced now, and we can find up to 2S + 1 peaks (cf. the cosine problem in figure 7.10).
Following our custom, we will try to draw an approximate formula over the exact curves. Recall that under the crude rotating-wave approximation Eq. (3.40), the time evolution of the non-diagonal elements is simply given by ̺ m±1,m ∝ exp[−i(∆ m ± 1m + R m±1m±1,mm )t]. Then the response in the frequency domain follows from Eq. (6.11) as χ x (ω) = 1 2 χ eq x m a + m i∆ m−1m + R m−1,m i(ω + ∆ m−1m ) + R m−1,m + a − m i∆ m+1m + R m+1,m i(ω + ∆ m+1m ) + R m+1,m (8.33) Here χ eq x is the equilibrium susceptibility, while R m±1,m ≡ R m±1m±1,mm are relaxation coefficients and the amplitudes a ± m depend essentially on the population difference between the states m and m ∓ 1.
This constitutes a generalization of equation (6.28) for S > 1 2 . This formula is plotted as the continuous line in figure 8.9, where we see that it fits quite well the numerics (which correspond to the improved RWA, Eq. (8.13), symbols).
If we now increase S or the coupling λ (that is, decreasing the quotient ∆ mm±1 /λ), we have observed that the crude RWA (with which we got the formula) starts to fail (see also [37] for more RWA failures). We could still try to get the classical limit of the susceptibility [96]. But there the formula would be valid only in the limit λ = 0, anyhow recovering Gekht's formula for the absorption line of a classical superparamagnet [97,98,96]. In any case, for the couplings and spin values in ordinary quantum superparamagnets, equation (8.33) could be a valuable tool to fit experimental absorption curves.

summary
In this last chapter, we have sought to apply all the formalism discussed along the thesis to the example of quantum superparamagnets. What we learnt in chapter 4 has been used to calculate dissipative corrections to the celebrated Brillouin law of quantum paramagnets. With the master equation of chapter 3 together with the methods presented in chapter 6, we have studied relaxation mechanisms in superparamagnets. We could contribute to the study of the phenomenon of tunnel of the magnetic moment, studying the nonlinear susceptibility. Finally, we derived approximate expressions for several quantities; the formula for the dynamical transverse susceptibility can be used to estimate the decoherence times in superparamagnets.

summary & conclusions
As it is the custom, let us conclude remarking the most important results presented in this thesis. • The dissipation theory is now formulated with a Hamiltonian: we can quantize it. We discussed the concept of quantum decoherence and derived quantum master equations for the reduced density matrix. The master equations used here are valid in the regimen of weak system-bath coupling, compared with the self-correlation times of the bath.
• We showed that the master equations are thermodynamically consistent, that is, that they have as stationary solution the reduced distribution obtained with quantum statistical mechanics. We used the opportunity to discuss the differences between equilibrium in the classical and quantum cases. In particular, in contrast to the classical distribution, the quantum reduced distribution depends on the damping. We concluded the equilibrium properties discussing how to properly define and compute thermodynamical quantities.
• With the phase-space formalism we closed the circle. We transformed the equations for the density matrix ̺, into equations for the Wigner distribution. Working in phase space and taking the limit → 0, the quantum equations reduce to their classical counterparts, explicitly relating the classical and quantum theories of dissipation.
• We discussed the continued-fraction technique to solve master equations (obtainment of stationary solutions and the linear and nonlinear responses). In spite of the demanding pre-treatment and specificity, when it does work it happens to be a quite useful method. In particular, it very suited to solve the Caldeira-Leggett master equation in phase space for nonlinear potentials. It also allows to study spin problems with S ≫ 1, a limit very demanding with other techniques.
• We applied the formalisms and the methods of the first part of the thesis to the problem of a damped particle in a periodic potential. The main results are: -We obtain the Wigner stationary solutions, and we see how the wave nature of quantum mechanics takes reflection in the negative zones in the "distribution", persisting at long times.
-Computing the currents in a cosine potential, we see how increasing the damping we get closer and closer to the classical limit.
-We tackled the problem of directional motion in potentials without spacereflection symmetry. We computed the first quantum corrections to the rectified current. The deviations from the classical limit can be understood appealing to two complementary quantum phenomena: tunnel and wave reflection.
-We obtained a estimation of the decoherence times and the inter-well jumps from the calculation of the linear susceptibility.
• Finally, we addressed the equilibrium and dynamical properties of quantum superparamagnets. The main results are: -We obtained dissipative corrections to the Brillouin magnetization law. In particular, we could address the entanglement between system and bath at very low temperatures from the equilibrium magnetization.
-We fully characterized the relaxation mechanisms in these systems: overbarrier dynamics, and intrawell processes in the potential wells. This allowed us to obtain a simple formula for the longitudinal linear response and closed form expressions for the relaxation times for any value of the spin S.
-We generalized our formalism to the nonlinear response. We showed how to detect tunnel by means of the nonlinear susceptibility.
-We produced a simple expression for the transverse susceptibility, which accounts for the absorption spectra in the ranges of parameters relevant for ordinary experiments.