Heat fluctuations in quantum spin chains

Gabriel T. Landi

Instituto de Física da 
Universidade de São Paulo

Encontro de física - Natal, RN
07 de setembro de 2016

Introduction

In collaboration with Dragi Karevski, Université de Lorraine.

Idea: to study the heat flux exchanged between two XX quantum spin chains.

Prepare the system in two halves at different temperatures (quantum quench).

r=0.5;
xmax = 7;
graph=Graphics[{

Text["(Subscript[T, 1])",{-1.,0}],
Text["(Subscript[T, 2])",{xmax+1,0.}],

Opacity[0.1],Red,EdgeForm[Directive[Black,Dashed]],
Table[Rectangle[{x-r,-r},{x+r+3,r}],{x,{0}}],

Opacity[0.1],Blue,EdgeForm[Directive[Black,Dashed]],
Table[Rectangle[{x-r-3,-r},{x+r,r}],{x,{xmax}}],
Opacity[1],
Black,
PointSize[0.02],
Table[Point[{x,0}],{x,0,xmax}]
},BaseStyle->{FontFamily->"Times",30},
ImageSize->600]

At t = 0 we turn on the two halves and monitor the flow of heat between the two bodies.

Unitary dynamics.

According to classical thermodynamics:

Heat should flow from the hot to the cold body.

The flow should continue until both bodies acquire the same temperature.

The total heat flow after a long time will be

Q = ((U((Subscript[T, 2]))-U((Subscript[T, 1])))/2)

The Jarzynski-Wójcik fluctuation theorem

In mesoscopic and microscopic systems, fluctuations become relevant.

We must now speak of Q as a random variable with a probability density (Subscript[P, t])(Q)

C. Jarzynski and D. K. Wójcik, PRL. 92, 230602 (2004)

Exchange fluctuation theorem (valid for weak coupling)

(((Subscript[P, t])(Q))/((Subscript[P, t])(-Q)))= (e^(Δβ Q)),    Δβ = (1/(Subscript[T, 1]))-(1/(Subscript[T, 2]))

“It is exponentially more likely to observe the flow of heat in the ‘right’ direction.”

Consequences:

⟨(e^(-Δβ Q))⟩ = 1           =>            ⟨Q⟩ ≥ 0            ((Subscript[T, 2])>(Subscript[T, 1]))

“On average, heat always flows in the ‘right’ direction.”

In practice, it is easier to study the characteristic function:

(Subscript[F, t])(r) = ⟨(e^( i r Q))⟩ = tr{(e^( i r (Subscript[H, 1])(t)))(e^(-i r (Subscript[H, 1])(0)))(Subscript[ρ, th])}

Then the distribution of heat is found from the inverse Fourier transform

(Subscript[P, t])(Q) = (1/(2π))∫ (e^(-i r Q))(Subscript[F, t])(r) dr

Goal of this work

To study (Subscript[P, t])(Q) for a many-body model.

Observe the transition micro -> meso -> macro.

Study the influence of a quantum phase transition.

GTL, D. Karevski, PRE 93 032122 (2016)

The XX chain

We consider two XX chains weakly coupled to each other

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

Each chain can be individually diagonalized with a Jordan-Wigner + sine transform

(Subscript[H, 1]) = (∑_( k)) (Subscript[λ, k]) (Subscript[a, k]^†)(Subscript[a, k]) 
(Subscript[H, 2]) = (∑_( k)) (Subscript[λ, k]) (Subscript[b, k]^†)(Subscript[b, k]) 
(Subscript[λ, k ])= h - 2  cos(k)

Here (Subscript[a, k]) and (Subscript[b, k]) are Fermionic operators.

The chains undergo a quantum phase transition at (Subscript[h, c]) = 2.

A gap opens in the excitation spectrum.

Plot[Evaluate@Table[h-2 Cos[k],{h,{0,1,2,3}}],{k,0,π},
BaseStyle->{FontFamily->"Times",27},
GridLines->{None,{0}},
FrameTicks->{{Automatic,Automatic},{{0,π/2,π},Automatic}},
FrameLabel->{"k","(Subscript[λ, k])"},
ImageSize->500,
PlotRangePadding->None,
GridLinesStyle->Directive[Gray, Dashed],
Epilog->{
Text["h = 0",Scaled[{0.9,0.44}]],
Text["h = 1",Scaled[{0.9,0.61}]],
Text["h = 2",Scaled[{0.9,0.75}]],
Text["h = 3",Scaled[{0.9,0.89}]],
Arrowheads[{-0.025,0.025}],
Arrow[{{0.2,0},{0.2,0.95}}],
Text["gap",{0.4,0.6}]
}
]

This model allows numerically exact solutions

The interaction potential becomes

V = (∑_( k, q))(Subscript[G, k,q])((Subscript[a, k]^†)(Subscript[b, q]) + (Subscript[b, q]^†)(Subscript[a, k])),       
(Subscript[G, k,q]) = g sin(L k) sin(q)                 g = ((2(Subscript[g, 0]))/(L+1))

The evolution after t = 0 will be unitary with

H = (∑_( k)) (Subscript[λ, k])((Subscript[a, k]^†)(Subscript[a, k]) + (Subscript[b, k]^†)(Subscript[b, k])) + (∑_( k, q))(Subscript[G, k,q])((Subscript[a, k]^†)(Subscript[b, q]) + (Subscript[b, q]^†)(Subscript[a, k]))

The Hamiltonian is quadratic.

Time evolution is entirely determined by the covariance matrix (in momentum space):

(Subscript[θ, k,q]) = ⟨(Subscript[c, k]^†)(Subscript[c, q])⟩    where c ∈ {a,b}

Or in position space

(Subscript[θ, i,j]) = ⟨(Subscript[c, i]^†)(Subscript[c, j])⟩

They satisfy

((dθ)/dt) = i [W,θ],      W = some matrix

Example: one fermion starting on the left.

Animate[MatrixPlot[Abs@θsol7[t],ImageSize->300],{t,0,200},DefaultDuration->34]

Analytical calculations

The exact Hamiltonian is

H = (∑_( k)) (Subscript[λ, k])((Subscript[a, k]^†)(Subscript[a, k]) + (Subscript[b, k]^†)(Subscript[b, k])) + (∑_( k, q))(Subscript[G, k,q])((Subscript[a, k]^†)(Subscript[b, q]) + (Subscript[b, q]^†)(Subscript[a, k]))       
(Subscript[G, k,q]) = g sin(L k) sin(q)                   g = ((2(Subscript[g, 0]))/(L+1))

The 2nd term is naturally weak since it depends ~1/L.

V is a “surface” term

(Subscript[H, 1]) and (Subscript[H, 2]) are “volume” terms.

We may treat V as a perturbation. We find:

H = (∑_( k)) { (Subscript[λ, k])((Subscript[a, k]^†)(Subscript[a, k]) + (Subscript[b, k]^†)(Subscript[b, k]))+ (Subscript[G, k ])((Subscript[a, k]^†)(Subscript[b, k]) + (Subscript[b, k]^†)(Subscript[a, k])) }
(Subscript[G, k]) = g sin(L k) sin(k) = (±1) g (sin^2)k

The entire problem factors into L independent subspaces, one for each pair of normal modes ((Subscript[a, k]),(Subscript[b, k])).

Only modes with the same momentum can exchange energy.

Characteristic function

The heat exchanged by each mode are independent random variables.

Q = (∑_( k )) (Subscript[Q, k])

(like a random walk)

Charac. Func. (Subscript[F, t])(r) will factor into the product of L independent functions

(Subscript[F, t])(r) = (∏_( k )) (Subscript[F, t])(r,k)
(Subscript[F, t])(r,k) = (Subscript[p, k]^0) + (e^( i r (Subscript[λ, k])))(Subscript[p, k]^(+)) + (e^(-i r (Subscript[λ, k])))(Subscript[p, k]^(-))

During the evolution, there are only 3 possible events for (Subscript[Q, k]):

(Subscript[p, k]^(+))  = (Subscript[n, k]^2)(1-(Subscript[n, k]^1))(sin^2)((Subscript[G, k])t)       
                (Subscript[p, k]^(-))  = (Subscript[n, k]^1)(1-(Subscript[n, k]^2))(sin^2)((Subscript[G, k])t)   
  (Subscript[p, k]^0) = 1 - (Subscript[p, k]^(+))-(Subscript[p, k]^(-))
  (Subscript[n, k]^(1,2)) = ((((e^((Subscript[λ, k])/(Subscript[T, 1,2])))+1))^(-1))

These probabilities individually satisfy the fluctuation theorem

((Subscript[p, k]^(+))/(Subscript[p, k]^(-))) = (e^(Δβ (Subscript[λ, k])))

Average heat

The average heat is

(Subscript[⟨Q⟩, t]) = (∑_( k)) (Subscript[λ, k]) ((Subscript[p, k]^(+)) - (Subscript[p, k]^(-))) = (∑_( k)) (Subscript[λ, k]) ((Subscript[n, k]^2)-(Subscript[n, k]^1)) (sin^2)((Subscript[G, k])t)

If we set (Subscript[T, 2])=T and (Subscript[T, 1]) = T + ΔT and convert the sum to an integral (thermodynamic limit), we get

(Subscript[⟨Q⟩, t]) = ((ΔT L)/π)(∫_0^( π))(Subscript[c, k])(T) (sin^2)(g t (sin^2)k) dk
(Subscript[c, k ])(T) = (Subscript[λ, k])((∂(Subscript[(OverBar[n]), k]))/(∂T)) = specific heat of each mode

At high temperatures we obtain a somewhat uninteresting behavior:

img=load["aveQ1.pdf",0.6];
lab[img,"High temp: T = 10"]

At low temperatures we see a strong influence of the quantum phase transition.

(*img1 = load["aveQ2a.png",0.45];
img2 = load["aveQ2b.png",0.45];
img3 = load["aveQ2c.png",0.45];
img4 = load["aveQ2d.png",0.45];
lab[Grid[{{lab[img1,"h = 0.5"],lab[img2,"h = 1.0"]},{lab[img3,"h = 1.5"],lab[img4,"h = 1.95"]}}],"T = 0.01"]
*)
Grid@Partition[Table[
Plot[(10^(3 )) AveQ[∞,h,0.01,t],{t,0,300},
PlotRange->{0,If[h==2.1,0.11,All]},
Frame->True,
PlotRangePadding->0,
(*GridLines->{None,{Evaluate[(10^(3 )) AveQ[∞,h,0.01,∞]]}},*)
GridLinesStyle->Directive[Black,Dashed],
BaseStyle->{FontFamily->"Times",18},
FrameLabel->{"gt","⟨Q⟩/LΔT (× (10^(-3)))"},
ImageSize->335,
AspectRatio->1,
ImagePadding->{{60,15},{53,10}},
PlotStyle->Black,
PlotPoints->200,
Epilog->{
Text[Style["h = "<>ToString@h,22],Scaled[{0.77,0.08}]]
}
],{h,{0.5,1.,1.5,1.95,2,2.1}}],2,2,1," "]

The average heat tends to the thermodynamic value

(lim_( t -> ∞ )) (Subscript[⟨Q⟩, t]) = (Subscript[Q, thermo]) = ((U((Subscript[T, 2]))-U((Subscript[T, 1])))/2)

But this does not mean that the system’s density matrix will reach equilibrium.

It is widely expected for integrable systems that the system will tend to a generalized Gibbs state.

see:  M. Collura & D. Karevski, PRB, 89(21) (2014)

Finite size effects

When the size of the system is fininte, there will be intereference from wavepackets that reflect at the boundaries of the sample.

img1=lab[load["aveQ3a.pdf",0.65],"T = 2, h = 0"];
img2=lab[load["aveQ3b.pdf",0.65],"T = 0.2, h = 0"];
img3=lab[load["aveQ3c.pdf",0.65],"T = 0.2, h = 1.8"];
Column[{img1,img2,img3}]

The distribution of heat as a function of time

We can find (Subscript[P, t])(Q) from an inverse Fourier transform of the characteristic functio

(Subscript[F, t])(r) = (∏_( k )){ (Subscript[p, k]^0) + (e^( i r (Subscript[λ, k])))(Subscript[p, k]^(+)) + (e^(-i r (Subscript[λ, k])))(Subscript[p, k]^(-)) }

This can be done using the FFT algorithm.

img1=lab[load["PQ_a.pdf",0.45],"L = 13"];
img2=lab[load["PQ_b.pdf",0.45],"L = 30"];
img3=lab[load["PQ_c.pdf",0.45],"L = 60"];
img4=lab[load["PQ_d.pdf",0.45],"L = 100"];
img5=lab[load["PQ_e.pdf",0.45],"L = 200"];
img6=lab[load["PQ_f.pdf",0.45],"L = 400"];
Grid[{{img1,img2},{img3,img4},{img5,img6}}]

Influence of finite size effects on the time evolution

Animate[Show[frames1[[i]],ImageSize->400],{i,1,Length@frames1,1},
AnimationRate->15]

Instabilities due to reflections at the boundaries

Animate[Show[frames2[[i]],ImageSize->400],{i,1,Length@frames2,1},
AnimationRate->22]

Next step

Couple the two halves to Lindblad baths

r=0.5;
xmax = 7;
graph=Graphics[{

Text["(Subscript[T, 1]), (Subscript[μ, 1])((Subscript[h, 1]))",{-2.4,0}],
Text["(Subscript[T, 2]), (Subscript[μ, 2])((Subscript[h, 2]))",{xmax+2.4,0.}],

Opacity[0.1],Red,EdgeForm[Directive[Black,Dashed]],
Table[Rectangle[{x-r,-r},{x+r+3,r}],{x,{0}}],

Opacity[0.1],Blue,EdgeForm[Directive[Black,Dashed]],
Table[Rectangle[{x-r-3,-r},{x+r,r}],{x,{xmax}}],
Opacity[1],
Black,
PointSize[0.02],
Table[Point[{x,0}],{x,0,xmax}]
},BaseStyle->{FontFamily->"Times",25},
ImageSize->600]

In collaboration with Pedro Guimarães and Mario de Oliveira (USP).

Was discussed in poster session, 05/09.

P. H. Guimarães, GTL and M. J. de Oliveira, submitted for publication (2016)

Thank you all for your attention.


and thanks FAPESP for the financial support.

GTL, D. Karevski, PRE 93 032122 (2016)

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[n, k]^1)-(Subscript[n, k]^2)) dk
(Subscript[J, E]) = ((4 (Subscript[g, 0]^2))/π)(∫_0^( π)) (Subscript[λ, k]) (sin^3)k ((Subscript[n, k]^1) - (Subscript[n, k]^2)) dk

Now the currents 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[n, k])(sin^3)k  dk
    (Subscript[J, E]) = δμ ((∂G)/(∂μ)) +δT ((∂G)/(∂T)),            G = ((4 (Subscript[g, 0]^2))/π)∫(Subscript[λ, k])(Subscript[n, 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

We attribute this to our choice of quantum master equation, which satisfies detailed balance.

P. H. Guimarães, GTL, M. J. de Oliveira, submitted for publication (2016).

Funções auxiliares

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

SetDirectory[NotebookDirectory[]]
<<"LinLib`";
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->{FontFamily->"Times",20},
ImageSize->400];

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";

fnames1=FileNames["*tiff","Frames_1"];
frames1=Import/@fnames1;

fnames2=FileNames["*tiff","Frames_2"];
frames2=Import/@fnames2;

DumpSave["θsol.mx",θsol7];

<<"θsol.mx";

Clear[AveQ]; 
AveQ[L_,h_,T_,t_]:=Module[{λ,k,s}, 

λ = h-2 Cos[k];
s = If[t===∞,(1/2),((Sin[t ((Sin[k])^2)])^2)];
If[L===∞,
(1/π)NIntegrate[((((λ/T)))^2)((E^(λ/T))/((((E^(λ/T))+1))^2))s,{k,0,π},Method->{Automatic,"SymbolicProcessing"->False}],
(1/L) Total@Table[((((λ/T)))^2)((E^(λ/T))/((((E^(λ/T))+1))^2))s,{k,(π/(L+1)),((L π)/(L+1)),(π/(L+1))}] 
]
];