Open quantum systems applied to condensed matter physics

Gabriel T. Landi
Institute of Physics, University of São Paulo

XXXIX Congresso Paulo Leal Ferreira
October 26th, 2016

Motivation: from micro to macro

The main goal of condesed matter physics is to understand the macroscopic properties of matter. e.g.,

The conductivity of a metal.

The specific heat of a gas.

The viscosity of a fluid.

The magnetization of a piece of iron.

etc.

Long ago we learned that, in order to do that, we must look at the microscopic world.

We now understand these properties as being emergent properties of the underlying microscopic interactions.

Thus, we must learn how to obtain macroscopic results from a microscopic theory.

In general that is a hard question:

The micro-world has ~(10^23) particles, all interacting like crazy (Think, for a second, about how complicated a glass of water is).

However, the great hero Josiah Willard Gibbs, continuing the work started by Maxwell and Boltzmann, found an answer for one specific class of systems: those in thermal equilibrium.

He found that the probability of finding the system at a quantum state |n⟩ is

(Subscript[P, n]) = ((e^(-β (Subscript[E, n])))/Z),    β = (1/((Subscript[k, B])T))

This formula represents a bridge between the two worlds:

It tells you how to compute macroscopic properties from knowledge of the microscopic energies.

According to Robert Millikan,

" Gibbs did to statistical mechanics and thermodynamics what Laplace did for celestial mechanics and Maxwell for eletromagnetism."

Please note: the Gibbs formula is atemporal: it does not matter how we arrived at equilibrium, all that matters is that we are there.

Example

load["Grafeno_screen.png", 0.8]

Prediction of the theory (this is a calculation any undergraduate student can do)

m = tanh(((μH)/((Subscript[k, B])T)))

Perfect agreement with experiment (F/C := μ):

load["Grafeno_screen2.png",0.6]

The Gibbs formula does not solve all problems, but... it tells you where to start.

All the information about the liquid-gas transition is contained within this integral.

Z = ∫exp[-(1/((Subscript[k, B])T))(∑_( pairs)) V((Subscript[r, i,j]))] (d^3)(Subscript[r, 1])(d^3)(Subscript[r, 2])...(d^3)(Subscript[r, N])

We have no idea how to compute it, but that is not surprising:

The integral is difficult because nature is difficult.

But at least we know what we should do.

And now we can look for alternative ways to tackle the problem:

Monte Carlo, cluster expansions, mean-field theories, etc.

Non-equilibrium systems

Now consider a physical system that is maintained in contact with two reservoirs at different temperatures.

Graphics[{
EdgeForm[Black],
Lighter[Gray,0.4],
Rectangle[{1,.5},{6,2.5}],

Lighter@Red,
Rectangle[{0,0},{1,3}],
Lighter@Blue,
Rectangle[{6,0},{7,3}],
Black,
Text["(Subscript[T, A])",{0.5,1.5}],
Text["(Subscript[T, B])",{6.5,1.5}],

Arrow[{{2.,1.5},{5.,1.5}}],
Text["J",{3.5,1.8}]
},
ImageSize->{400,Automatic},
BaseStyle->20]

This system will never reach equilibrium, because the warmer reservoir will be constantly giving out heat to the colder reservoir.

NESS: "Non-equilibrium steady-state".

Its a "steady-state" because things are no longer changing in time.

But its not thermal equilibrium.

What are the rules governing non-equilibrium systems?

What is the analog of the Gibbs formula for non-equilibrium systems?

What are the rules governing non-equilibrium systems?

What is the analog of the Gibbs formula for non-equilibrium systems?

We have no idea!

We do not have a general framework for handling non-equilibrium systems.

If we want to study the NESS, we must go back to the fundamental laws (Newton, Schrödinger).

Even though we are interested in something atemporal, we need to include the full dynamics!

There is no unified framework for dealing with non-equilibrium physics.

Classical systems:

Stochastic processes: Langevin/Fokker-Planck equation

Molecular dynamics.

Quantum systems:

Green-Kubo (linear response)

Landauer-Büttiker

Boltzmann equation.

My research interests: Lindblad dynamics

My main research interest is in building a consistent framework to study (mostly) the NESS in condensed matter systems.

And I do that by applying techniques initially developed in atomic physics and quantum optics.

In particular, today I want to focus on Lindblad dynamics

This has applications in

Electron transport

Heat transport

Spin transport (magnon currents).

Motivation: the quantum harmonic oscillator

(Subscript[E, n]) = ℏω(n+1/2)

xr =linspace[5,100,10]//Sqrt;
Plot[(x^2),{x,-10.5,10.5},
Prolog->{Blue,Table[Line[{{-x,(x^2)},{x,(x^2)}}],{x,xr}]},
AspectRatio->1,
Frame->False,
AxesOrigin->{0,0},
Axes->True,
AxesStyle->Arrowheads[{0.05}],
ImageSize->300,
BaseStyle->{FontFamily->"Times",20},
Ticks->None,
AxesLabel->{"x","(1/2)(ω^2)(x^2)"},
ImagePadding->{{20,20},{50,50}}
]

The contact with a heat bath will cause transitions between the different energy levels.

Imagine that QHO = diatomic molecule and Heat Bath = Electromagnetic field.

The molecule can absorb and emit photons.

The bath determines the average energy level (OverBar[N]).

nbrange ={0,0.5,1,2,4};
n0 = 4;
steps = 6000;
Δt = 0.002;
Clear[nv,nb];
Do[
nv = ConstantArray[0,steps];
nv[[1]] = n0;
Do[
nv[[i+1]]=RandomChoice[
{
Δt  nb (nv[[i]]+1),
Δt (nb+1) nv[[i]],
1-Δt(nb (nv[[i]]+1)+  (nb+1) nv[[i]])}->{nv[[i]]+1,nv[[i]]-1,nv[[i]]}],{i,1,steps-1}];
NV[nb]=nv;
,{nb,nbrange}];

Animate[
ListLinePlot[NV[nb],
(*GridLines->{None,{3}},*)
AspectRatio->1/5,
ImageSize->700,
BaseStyle->{FontFamily->"Times",20},
Frame->True,
Axes->None,
FrameLabel->{"tempo","n"},
PlotRange->{{Max[0,t-1000],Max[2000,t+1000]},{-0.5,Max[NV[nb]]}},
PlotRangePadding->None,
PlotStyle->Black,
Epilog->{Red,
Text[Style["●",16],{t,NV[nb][[t]]}]}
]

,{t,1,Length@nv,1},{{nb,4,Style["(OverBar[N])",22]},nbrange,ControlType->Setter},DefaultDuration->40,AnimationRunning->False]

Lindblad dynamics

Idea: to describe the contact with a bath in a fully quantum mechanical way.

Extend von Neumann’s equation with an additional term:

((dρ)/dt)= -i[H,ρ] + 𝒟(ρ)

𝒟(ρ) is called the dissipator and describes the contact with the heat bath.

Must be completely positive trace preserving.

tr(𝒟(ρ))=0

Main work of Lindblad

𝒟(ρ) = (∑_( α)) ((Subscript[L, α])ρ (Subscript[L, α]^†) - (1/2){(Subscript[L, α]^†)(Subscript[L, α]),ρ})

Constructing the dissipator may be a difficult task.

We will discuss here two approaches:

Phenomenological dissipators.

Microscopic derivations.

The most famous example: the quantum harmonic oscillator

The dissipator for a quantum harmonic oscillator with H = ω a†a is

𝒟(ρ) = γ(1+(OverBar[n]))[a ρ a† -(1/2){a†a,ρ}] + γ (OverBar[n])[a†ρ a -(1/2){a a†,ρ}]

where γ = coupling strength to the bath and

(OverBar[n]) =(1/((e^(βω))-1))

This dissipator correctly thermalizes the system: it has only one fixed point:

𝒟((e^(-β H))) = 0

It also satisfies detailed balance

((γ(1+(OverBar[n])))/(γ (OverBar[n]))) = (e^(β ω))

This is a physical property, related to the fact that the underlying dynamics (Schrödinger’s equation for the system+bath) is time-reversal invariant.

This dissipator produces

((d⟨a†a⟩)/dt) = γ((OverBar[n])-⟨a†a⟩)

The average occupation of the oscillator relaxes toward its equilibrium value (Subscript[⟨a†a⟩, ∞])=(OverBar[n]).

Two coupled oscillators

Now consider two coupled oscillators, with bosonic operators a and b. The Hamiltonian is taken to be

H = ω(a†a+b†b) -g(a†b+b†a)

Each system is now assumed to be connected to a dissipator, just like before:

(Subscript[𝒟, a])(ρ) = γ(1+(Subscript[(OverBar[n]), a]))[a ρ a† -(1/2){a†a,ρ}] + γ (Subscript[(OverBar[n]), a])[a†ρ a -(1/2){a a†,ρ}]
(Subscript[𝒟, b])(ρ) = γ(1+(Subscript[(OverBar[n]), b]))[b ρ b† -(1/2){b†b,ρ}] + γ (Subscript[(OverBar[n]), b])[b†ρ b -(1/2){b b†,ρ}]

The complete dynamics of the system will then be governed by

((dρ)/dt)= -i[H,ρ] + (Subscript[𝒟, a])(ρ)+(Subscript[𝒟, b])(ρ)

First problem: what should we take for (Subscript[(OverBar[n]), a,b])? Let’s leave them generic for now.

Suppose first that the two baths are at the same temperature ((Subscript[(OverBar[n]), a])=(Subscript[(OverBar[n]), b])). Even in this case, these dissipators no longer take the system correctly to equilibrium!

Thermal equilibrium is ρ = (e^(-β H)). But now the two oscillators interact so we cannot factor (e^(-β H)) as a product.

Moral of the story

If you want the dissipators to model Heat baths, then their shape will depend on H

Change any tiny detail on H and you should change your 𝒟s

The example above is what I call “phenomenological dissipators” or “local dissipators”

They should give reasonable answers when the two oscillators are weakly coupled.

Note: we no longer can define (Subscript[T, a]) and (Subscript[T, b]). But we can still play with (Subscript[(OverBar[n]), a]) and (Subscript[(OverBar[n]), b]).

Let’s examine the NESS of this problem. We find

⟨a†a⟩= (((Subscript[(OverBar[n]), a])+(Subscript[(OverBar[n]), b]))/2)+((γ^2)/((γ^2)+(g^2)))((((Subscript[(OverBar[n]), a])-(Subscript[(OverBar[n]), b])))/2)
⟨b†b⟩= (((Subscript[(OverBar[n]), a])+(Subscript[(OverBar[n]), b]))/2)-((γ^2)/((γ^2)+(g^2)))((((Subscript[(OverBar[n]), a])-(Subscript[(OverBar[n]), b])))/2)
⟨a†b⟩ =((2i γ g)/((γ^2)+(g^2))) ((((Subscript[(OverBar[n]), a])-(Subscript[(OverBar[n]), b])))/2)

The two oscillators are correlated: ⟨a†b⟩ ≠ 0.

This produces a flow of energy

J ∝⟨a†b⟩ = ((ω γ g)/((γ^2)+(g^2))) ((((Subscript[(OverBar[n]), a])-(Subscript[(OverBar[n]), b])))/2)

There will continue to be a flux as long as the bath occupations are different.

We will come back to this problem later and compare with a microscopic model.

Fun fact: this result for the current is actually the same for a 1D chain of N oscillators.

This defines a ballistic current (independent of size).

It is related to the fact that in these models the excitations propagate like waves.

Spin 1/2 chains

Now I want to talk about spin chains.

Consider a single spin 1/2 particle with H = - μ B (Subscript[σ, z])

A local dissipator correctly thermalizing this particle is

𝒟(ρ) = (γ/2)(1+f)[(SubPlus[σ])ρ (SubMinus[σ])-(1/2){(SubMinus[σ])(SubPlus[σ]),ρ}]+(γ/2)(1-f)[(SubMinus[σ])ρ (SubPlus[σ])-(1/2){(SubPlus[σ])(SubMinus[σ]),ρ}]

where

f= tanh(((μ B)/((Subscript[k, B])T)))

This dissipator takes the system to the correct thermal state: 𝒟((e^(β μ B (Subscript[σ, z]))))=0.

Moreover, it again satisfies detailed balance

((γ(1+f))/(γ(1-f))) = (e^(2β μ B))

Quantum spin chains

Now we apply this to a 1D quantum XXZ spin chain.

H = (1/2)(∑_( i = 1)^( N-1))((Subscript[σ, i]^x)(Subscript[σ, i+1]^x)+(Subscript[σ, i]^y)(Subscript[σ, i+1]^y) + Δ (Subscript[σ, i]^z)(Subscript[σ, i+1]^z))

When Δ = 0 (XX chain) the model is exactly solvable.

It maps into a model of free fermions.

The Δ-term can be seen as an interaction between quasi-particles

We use the local dissipators

(Subscript[𝒟, i])(ρ) = (γ/2)(1+(Subscript[f, i]))[(Subscript[σ, i]^(+))ρ (Subscript[σ, i]^(-))-(1/2){(Subscript[σ, i]^(-))(Subscript[σ, i]^(+)),ρ}]+(γ/2)(1-(Subscript[f, i]))[(Subscript[σ, i]^(-))ρ (Subscript[σ, i]^(+))-(1/2){(Subscript[σ, i]^(+))(Subscript[σ, i]^(-)),ρ}]

with i = 1 or N.

load["diagram.png",0.5]

We are interested here in the spin current, defined from the continuity equation

((d⟨(Subscript[σ, i]^z)⟩)/dt)=(Subscript[J, i-1])-(Subscript[J, i])

In this model it becomes

(Subscript[J, i])=⟨(Subscript[σ, i]^x)(Subscript[σ, i+1]^y)-(Subscript[σ, i]^y)(Subscript[σ, i+1]^x)⟩

In the NESS ((d⟨(Subscript[σ, i]^z)⟩)/dt)=0  ⟶ (Subscript[J, i]) = J  for all sites (the current is uniform through the chain).

XX chain (Δ =0)

In the case of an XX chain the current becomes simply

J = (γ/(1+(γ^2)))((((Subscript[f, 1])-(Subscript[f, N])))/2)

It is almost identical to the case of the quantum harmonic oscillator.

D. Karevski and P. Platini, Phys. Rev. Lett. 102, 207207 (2009)
   GTL and D. Karevski, Phys. Rev. E. 93, 032122 (2016)

Interpretation: the XX chain behaves like free fermions and the spin current is ballistic.

You inject particles on the left and collect on the right.

The particles just travel freely through the chain, so the total current does not depend on the size.

Δ +≠ 0: spin rectification

We have carried out numerical studies of the current under non-uniform fields

H = (1/2)(∑_( i = 1)^( N-1))((Subscript[σ, i]^x)(Subscript[σ, i+1]^x)+(Subscript[σ, i]^y)(Subscript[σ, i+1]^y) + Δ (Subscript[σ, i]^z)(Subscript[σ, i+1]^z)) + (∑_( i = 1)^( N)) (Subscript[h, i])(Subscript[σ, i]^z)

We choose the field as a linear profile, varying linearly between -h and +h.

This produces an asymmetry in the problem.

load["rectification_flux.pdf",0.5]

The spin current is now rectified.

If you reverse the baths, you reverse the sign of the current but the magnitude does not necessarily remains the same.

We may define the rectification coefficient

R = ((J(f) + J(-f))/(J(f) - J(-f)))

It measures the asymmetry of the current.

load["rectification_2.png",0.6]

The Rectification is zero when Δ = 0, even if the system is assymetric.

GTL, et. al., Phys. Rev. E. 90, 042142 (2014)
   L. Schuab, et. al. Phys. Rev. E. 94 042122 (2016)

Heisenberg chain (Δ = 1)

We have been able to obtain an exact solution, in the presence also of boundary fields.

H = (1/2)(∑_( i = 1)^( N-1))((Subscript[σ, i]^x)(Subscript[σ, i+1]^x)+(Subscript[σ, i]^y)(Subscript[σ, i+1]^y) + (Subscript[σ, i]^z)(Subscript[σ, i+1]^z))+h((Subscript[σ, 1]^z)-(Subscript[σ, N]^z))

Column[{load["drawing.pdf",0.4],load["drawing2.pdf",0.4]}]

The solution is based on the idea of matrix product states.

ρ = S S†
S = ⟨0|(Ω^(⊗N))|0⟩
Ω = (Subscript[S, z])(σ^z)+(SubPlus[S])(SuperPlus[σ])+(SubMinus[S])(SuperMinus[σ])

Each spin site is augmented with an auxiliary space, modeled by operators (Subscript[S, i]) satisfying the SU(2) algebra.

This allowed us to find an exact method to compute the current for arbitrary sizes.

D. Karevski, et. al., Phys. Rev. Lett. 110, 047201 (2013)
   GTL and D. Karevski, Phys. Rev. B. 91, 174422 (2015)

Here are the results for h = 0

load["Heisen1.png",0.7]

Here are the results as a function of the boundary field h

load["Heisen2.png",0.7]

“Semi-microscopic theories”

Now I want to turn to microscopic derivations.

Instead of using a phenomenological dissipator, we start with a microscopic theory of the system-bath interaction.

To warm-up we start with a “semi-microscopic theory”.

Consider again the XX chain: it can be mapped into a fermionic chain

H   = (1/2)(∑_( i = 1)^( N-1))((Subscript[σ, i]^x)(Subscript[σ, i+1]^x)+(Subscript[σ, i]^y)(Subscript[σ, i+1]^y) )
= (∑_( i = 1)^( N-1))((Subscript[η, i]^†)(Subscript[η, i+1]) + (Subscript[η, i+1]^†)(Subscript[η, i]))

We ask what is the dissipator that would

Take this chain to the correct thermal equilibrium state (e^(-β H)) and

Do so while satisfying detailed balance

We can construct this by hand. First we diagonalize the chain in Fourier space

(Subscript[c, k]) = (∑_( i = 1)^( N)) (Subscript[S, i,k]) (Subscript[η, i]),          (Subscript[S, i,k]) = (Sqrt[(2/(N+1))]) sin(((π i k)/(N+1)))

The Hamiltonian then becomes

H = (∑_( k)) (Subscript[ϵ, k ])(Subscript[c, k]^†)(Subscript[c, k]),        (Subscript[ϵ, k]) = cos(k)

This is now a problem of uncoupled Fermions. We know the dissipator which does the job for a single fermion. It reads

(Subscript[𝒟, k])(ρ) = γ (1-(Subscript[f, k])) [ (Subscript[c, k])ρ (Subscript[c, k]^†) - (1/2){(Subscript[c, k]^†)(Subscript[c, k]),ρ}] + γ (Subscript[f, k]) [(Subscript[c, k]^†)ρ (Subscript[c, k]) -(1/2){(Subscript[c, k])(Subscript[c, k]^†),ρ}]
(Subscript[f, k ])= (1/((e^(((Subscript[ϵ, k])-μ)/T))+1))

So a way to thermalize the total chain would be through the master equation

((dρ)/dt)=-i [H,ρ] + (∑_( k ))(Subscript[𝒟, k])(ρ)

Two coupled chains

But we want to study transport so we divide a chain in two or three parts and couple the different parts to a bath in this way.

load["multi_site_1.png",0.5]

This problem can be solved exactly if we assume that the different chains are weakly coupled.

The profile of the inner chain correctly reproduces those of the outside chains when γ is small

load["multi_site_2.png",0.4]

Particle and energy current

The particle and energy currents in the steady-state and in the thermodynamic limit will be

(Subscript[J, N]) = ((4 (Subscript[g, 0]^2))/π)(∫_0^( π))(sin^3)k ((Subscript[f, k]^1)-(Subscript[f, k]^2)) dk
(Subscript[J, E]) = ((4 (Subscript[g, 0]^2))/π)(∫_0^( π)) (Subscript[ϵ, k]) (sin^3)k ((Subscript[f, k]^1) - (Subscript[f, k]^2)) dk

They have the structure of the Landauer-Büttiker formula.

Both depend on temperatures (Subscript[T, 1,2]) and chemical potentials (Subscript[μ, 1,2])

Expanding for infinitesimal unbalances, we obtain

(Subscript[J, N]) = δμ ((∂F)/(∂μ)) +δT ((∂F)/(∂T)),          F = ((4 (Subscript[g, 0]^2))/π)∫(Subscript[f, k])(sin^3)k  dk
    (Subscript[J, E]) = δμ ((∂G)/(∂μ)) +δT ((∂G)/(∂T)),            G = ((4 (Subscript[g, 0]^2))/π)∫(Subscript[ϵ, k])(Subscript[f, k])(sin^3)k  dk

Manipulate[
lab[Grid[{{plμ[LL],plT[LL]},{pEμ[LL],pET[LL]}}],"L = "<>ToString@LL],{{LL,1,"L"},LLrange},ControlType->Slider]

Onsager reciprocal relations

The heat flux is defined as

(Subscript[J, Q]) = (Subscript[J, E]) - μ (Subscript[J, N])

We found that (Subscript[J, N]) and (Subscript[J, Q]) can be put in Onsager’s canonical form

(({{(Subscript[J, N])}, {(Subscript[J, Q])}})) = (({{(Subscript[L, 11]), (Subscript[L, 12])}, {(Subscript[L, 21]), (Subscript[L, 22])}}))(({{(δμ)/T}, {-δ(1/T)}}))

The Onsanger coefficients are all obtained analytically:

(Subscript[L, 11]) = T ((∂F)/(∂μ))               (Subscript[L, 12]) = (T^2)((∂F)/(∂T))
(Subscript[L, 21]) = T (((∂G)/(∂μ))- μ ((∂F)/(∂μ)))         (Subscript[L, 22]) = (T^2)(((∂G)/(∂T))-μ ((∂F)/(∂T)))

Moreover, they satisfy the reciprocal relations

(Subscript[L, 21]) = (Subscript[L, 12])
(Subscript[L, 11])(Subscript[L, 22]) - (Subscript[L, 21])(Subscript[L, 12 ])≥ 0

Reciprocal relations are related to detailed balance.

The original dissipators satisfy detailed balance.

But when we couple the two chains, DB is violated.

It only continues to hold approximately when γ is small.

P. H. Guimarães, GTL and M. J. de Oliveira, Phys. Rev. E. 94 032139 (2016)

Full microscopic theory for a many-body theory

We have just been able to construct a full microscopic derivation for a many-body theory.

We consider the 1D bosonic tight-binding chain

H=ϵ (∑_( i =1)^( N))(Subscript[a, i]^†)(Subscript[a, i ]) -g(∑_( i = 1)^( N-1))((Subscript[a, i]^†)(Subscript[a, i+1]) + (Subscript[a, i+1]^†)(Subscript[a, i]))

The first and last sites are coupled to heat baths (i=1 or N)

(Subscript[H, B,i]) = (∑_( ℓ ))(Subscript[Ω, ℓ]) (Subscript[b, ℓ,i]^†)(Subscript[b, ℓ,i])
(Subscript[H, int,i]) = (∑_( ℓ ))(Subscript[c, ℓ]) ((Subscript[a, i])+(Subscript[a, i]^†))((Subscript[b, ℓ,i])+(Subscript[b, ℓ,i]^†))

Even though the coupling is local (only at the first and last sites), we find that the Lindblad dissipator will be intrinsically non-local, affecting all normal modes of the system:

(Subscript[𝒟, i])(ρ) = (∑_( k ))(Subscript[S, i,k]^2)(Subscript[γ, i,k]) (1 + (Subscript[(OverBar[n]), i,k]))[(Subscript[a, k]) ρ (Subscript[a, k]^†) -(1/2){(Subscript[a, k]^†)(Subscript[a, k]),ρ}] + (Subscript[S, i,k]^2)(Subscript[γ, i,k]) (Subscript[(OverBar[n]), i,k])[(Subscript[a, k]^†)ρ (Subscript[a, k]) -(1/2){(Subscript[a, k]) (Subscript[a, k]^†),ρ}]

where (Subscript[γ, k]) is the bath spectral density

(Subscript[γ, k]) =2 π (∑_( ℓ)) (Subscript[c, ℓ]^2) δ((Subscript[ϵ, k]) - (Subscript[Ω,  ℓ]))

This dissipator correctly thermalizes the entire system while satisfying detailed balance.

Moreover, it has information about the position of which site is coupled to the bath.

Therefore, it allows us to study the currents through the system.

They again have the form of Landauer’s formula

(Subscript[J, N]) = (1/2)(∑_( k ))(Subscript[γ, k]) ( (Subscript[(OverBar[n]), 1,k])- (Subscript[(OverBar[n]), N,k]))
(Subscript[J, E]) = (1/2)(∑_( k ))(Subscript[ϵ, k])(Subscript[γ, k]) ( (Subscript[(OverBar[n]), 1,k])- (Subscript[(OverBar[n]), N,k]))

The transfer integral is now (Subscript[γ, k]), the system-bath spectral density

J. P. Santos and GTL, arXiv:1610.05126 (2016)

Conclusions

Main goal of condensed matter is to go from micro to macro.

For systems in thermal equilibrium, Gibbs told us the answer:

(Subscript[P, n]) = ((e^(-(Subscript[βE, n])))/Z)

This formula is amazing: working with non-equilibrium physics has taught me to value the equilibrium formalism.

NESS:

Much more complicated: must start with stochastic dynamics.

Enormous number of applications in condensed matter.

Lindblad dynamics:

Offer a rich platform for fully quantum mechanical studies.

Phenomenological theories give reasonable results, but are of limited extent.

Microscopic derivations are more difficult, but now I think we got it!

Next steps:

Full microscopic theory for the XX spin chain and the Ising model in a transverse field.

Study effects of quantum phase transition.

Magnons (modeled as bosons): study the magnon current.

Quantum field theory description and path integrals.

Microscopic theory for electrons: much harder because baths are usually bosonic.

Thank you all for your attention.

Funções auxiliares

(*Show[Import["","Image",ImageSize->600]]*)

SetDirectory[NotebookDirectory[]];
<<"LinLib.wl";
<<"CustomTicks`";
load[filename_,size_]:=Show[Import[filename],ImageSize->Scaled[size]];
lab[img_,lbl_]:=Labeled[img, lbl,Top,LabelStyle->{FontFamily->"Times",26,Blue}];


SetOptions[Plot,Frame->True,Axes->False,BaseStyle->20,ImageSize->400,PlotStyle->{Black}];

SetOptions[InputNotebook[], 
DefaultNewCellStyle -> "Item",
ShowCellLabel->"False",
CellGrouping->Manual,
FontFamily->"Times",
DefaultNewCellStyle->{"Text",FontFamily->"Times"},
BaseStyle->{FontFamily->"Times"},
MultiLetterItalics->False,
SingleLetterItalics->Automatic
]

Clear[ℐ];
ℐ[La_,γ_,k_,t_:1]:=Module[{a = (γ^2)+ (t^2)((Sin[k])^2),b = 2 γ t Cos[k],qr},
qr = Range[(π/(La+1)),((π La)/(La+1)),(π/(La+1))];
(((Sin[k])^2)/(La+1))Sum[(((Sin[q])^2)/((γ^2)+(t^2)(((Cos[k]-Cos[q]))^2))),{q,qr}]
]


Clear[ℐtl];
ℐtl[γ_,k_,t_:1]:=Module[{a = (γ^2)+ (t^2)((Sin[k])^2),b = 2 γ t Cos[k]},
(1/(t^2))( (1/(γ(Sqrt[2])))(Sqrt[a+(Sqrt[(a^2)+(b^2)])])-1)((Sin[k])^2) 
]


Clear[𝕁];
𝕁[La_,{μa_,μc_},{Ta_,Tc_},γ_:1]:=Module[{mr,𝕖,𝕟},

mr = Range[(π/(La+1)),((π La)/(La+1)),(π/(La+1))]; 
𝕖[k_]:=-2 Cos[k];
𝕟[μ_,T_,k_]:=If[T==0,HeavisideTheta[μ-𝕖[k]],(1/(Exp[(((𝕖[k]-μ))/T)]+1))];
 ((4γ (*(g^2)*))/(((La+1))^2))Total@Flatten@Table[((((Sin[k])^2)((Sin[q])^2)(𝕟[μa,Ta,k]- 𝕟[μc,Tc,k]))/((γ^2)+(((Cos[k]-Cos[q]))^2)))    ,{k,mr},{q,mr}]
];

Clear[𝕁tl];
𝕁tl[{μa_,μc_},{Ta_,Tc_},γ_:1]:=Module[{mr,𝕖,𝕟,a,b,int},

𝕖[k_]:=-2 Cos[k];
𝕟[μ_,T_,k_]:=If[T==0,HeavisideTheta[μ-𝕖[k]],(1/(Exp[(((𝕖[k]-μ))/T)]+1))];

a = (γ^2)+ ((Sin[k])^2);
b = 2 γ Cos[k];
(*int =( -1  + (1/(2γ))((Sqrt[a+I b])+(Sqrt[a-I b])))((Sin[k])^2);*)
int =( (1/(γ(Sqrt[2])))(Sqrt[a+(Sqrt[(a^2)+(b^2)])])-1)((Sin[k])^2);

 ((4γ (*(g^2)*))/π) NIntegrate[(𝕟[μa,Ta,k]- 𝕟[μc,Tc,k]) (int),{k,0,π}]
];

Clear[𝔽];
𝔽[La_,μ_,T_,γ_:1]:=Module[{mr,𝕖,𝕟},

mr = Range[(π/(La+1)),((π La)/(La+1)),(π/(La+1))]; 
𝕖[k_]:=-2 Cos[k]; 
𝕟[k_]:=If[T==0,HeavisideTheta[μ-𝕖[k]],(1/(Exp[(((𝕖[k]-μ))/T)]+1))] ;
 ((4γ (*(g^2)*))/(((La+1))^2))Total@Flatten@Table[((((Sin[k])^2)((Sin[q])^2)(𝕟[k]))/((γ^2)+(((Cos[k]-Cos[q]))^2)))    ,{k,mr},{q,mr}]
];

Clear[𝔽tl];
𝔽tl[μ_,T_,γ_:1]:=Module[{mr,𝕖,𝕟,a,b,int},

𝕖[k_]:=-2 Cos[k];
𝕟[k_]:=If[T==0,HeavisideTheta[μ-𝕖[k]],(1/(Exp[(((𝕖[k]-μ))/T)]+1))];

a = (γ^2)+ ((Sin[k])^2);
b = 2 γ Cos[k];
int =( (1/(γ(Sqrt[2])))(Sqrt[a+(Sqrt[(a^2)+(b^2)])])-1)((Sin[k])^2);

 ((4γ (*(g^2)*))/π) NIntegrate[(𝕟[k]) (int),{k,0,π}]//Chop
];

ℐEtl[γ_,k_,t_:1]:=Module[{F,a,b},

a = (γ^2)+(t^2)((Sin[k])^2);
b = 2 γ t Cos[k];

F=(4/t)Cos[k] - ((Sqrt[2])/(γ (t^2)))(t Cos[k](Sqrt[(Sqrt[(a^2)+(b^2)])+a]) + γ  ((Cos[k])/(Sqrt[((Cos[k])^2)]))(Sqrt[(Sqrt[(a^2)+(b^2)])-a]));

(*{-2t  Cos[k] ℐtl[γ,k,t] , ((Sin[k])^2)F,-2t  Cos[k] ℐtl[γ,k,t] + ((Sin[k])^2)F}*)
-2t  Cos[k] ℐtl[γ,k,t] + ((Sin[k])^2)F
];
Clear[𝕁E];
𝕁E[La_,{μa_,μc_},{Ta_,Tc_},γ_:1]:=Module[{mr,𝕖,𝕟},

mr = Range[(π/(La+1)),((π La)/(La+1)),(π/(La+1))]; 
𝕖[k_]:=-2 Cos[k];
𝕟[μ_,T_,k_]:=If[T==0,HeavisideTheta[μ-𝕖[k]],(1/(Exp[(((𝕖[k]-μ))/T)]+1))];
 ((2γ (*(g^2)*))/(((La+1))^2)) Total@Flatten@Table[(((Sin[k])^2)((Sin[q])^2)(𝕟[μa,Ta,k]- 𝕟[μc,Tc,k])(𝕖[k]+𝕖[q]))/((γ^2)+(((Cos[k]-Cos[q]))^2))    ,{k,mr},{q,mr}]
];

Clear[𝕁Etl];
𝕁Etl[{μa_,μc_},{Ta_,Tc_},γ_:1]:=Module[{mr,𝕖,𝕟,ie},

𝕖[k_]:=-2 Cos[k];
𝕟[μ_,T_,k_]:=If[T==0,HeavisideTheta[μ-𝕖[k]],(1/(Exp[(((𝕖[k]-μ))/T)]+1))];
ie = ℐEtl[γ,k];

 ((2γ (*(g^2)*))/π) NIntegrate[ie(𝕟[μa,Ta,k]- 𝕟[μc,Tc,k])  ,{k,0,π}] 
];

LLrange= {1,4,10,20,40,60,80,100,120,160};
δμ = 0.001;
TT = 0.02;
Jtl = Chop@Quiet@Table[{μ,𝕁tl[{μ+δμ,μ},{TT,TT}]/δμ},{μ,linspace[-3,3,100]}];

Do[
plμ[LL] = Show[
Plot[Evaluate[{𝕁[LL,{μ + ((δμ)/2),μ- ((δμ)/2)},{TT,TT}]/δμ}],{μ,-3,3},
PlotRange->{0,All},
AspectRatio->1,
ImagePadding->{{50,10},{60,10}},
FrameLabel->{"μ",None},
PlotLabel->"(Subscript[J, N]),μ",
BaseStyle->{FontFamily->"Times",20},
ImageSize->260,
PlotStyle->Black,
Frame->True,
PlotStyle->Directive[Black]
],
ListLinePlot[Jtl,PlotStyle->Directive[Red,Dashed]]
]; ,{LL,LLrange}]

δT = 0.001;
TT = 0.02;
JtlT = Chop@Quiet@Table[{μ,𝕁tl[{μ,μ},{TT+δT,TT}]/δT},{μ,linspace[-3,3,100]}];

Do[
plT[LL] = Show[
Plot[Evaluate[{𝕁[LL,{μ,μ},{TT+((δT)/2),TT-((δT)/2)}]/δT}],{μ,-3,3},
PlotRange->All,
AspectRatio->1,
ImagePadding->{{50,10},{60,10}},
FrameLabel->{"μ",None},
PlotLabel->"(Subscript[J, N]),T",
BaseStyle->{FontFamily->"Times",20},
ImageSize->260,
PlotStyle->Black,
Frame->True
],
ListLinePlot[JtlT,PlotStyle->Directive[Red,Dashed]]
]; ,{LL,LLrange}]


δμ = 0.001;
TT = 0.02;
JEtl = Chop@Quiet@Table[{μ,𝕁Etl[{μ+δμ,μ},{TT,TT}]/δμ},{μ,linspace[-3,3,100]}];

Do[
pEμ[LL] = Show[
Plot[Evaluate[{𝕁E[LL,{μ + ((δμ)/2),μ- ((δμ)/2)},{TT,TT}]/δμ}],{μ,-3,3},
PlotRange->All,
AspectRatio->1,
ImagePadding->{{50,10},{60,10}},
FrameLabel->{"μ",None},
PlotLabel->"(Subscript[J, E]),μ",
BaseStyle->{FontFamily->"Times",20},
ImageSize->260,
PlotStyle->Black,
Frame->True,
PlotStyle->Directive[Black]
],
ListLinePlot[JEtl,PlotStyle->Directive[Red,Dashed]]
]; ,{LL,LLrange}]


δT = 0.001;
TT = 0.02;
JEtlT = Chop@Quiet@Table[{μ,𝕁Etl[{μ,μ},{TT+δT,TT}]/δT},{μ,linspace[-3,3,100]}];

Do[
pET[LL] = Show[
Plot[Evaluate[{𝕁E[LL,{μ,μ},{TT+((δT)/2),TT-((δT)/2)}]/δT}],{μ,-3,3},
PlotRange->All,
AspectRatio->1,
ImagePadding->{{50,10},{60,10}},
FrameLabel->{"μ",None},
PlotLabel->"(Subscript[J, E]),T",
BaseStyle->{FontFamily->"Times",20},
ImageSize->260,
PlotStyle->Black,
Frame->True,
PlotStyle->Directive[Black]
],
ListLinePlot[JEtlT,PlotStyle->Directive[Red,Dashed]]
]; ,{LL,LLrange}]
DumpSave["plots.mx",{plμ,plT,pEμ,pET}];

LLrange= {1,4,10,20,40,60,80,100,120,160};
<<"plots.mx";