Sistemas quânticos abertos fora do equilíbrio

Gabriel T. Landi

Universidade Estadual de Londrina
20 de maio de 2016

Sistemas abertos: atrito/viscosidade

Perguntas simples as vezes possuem respostas complicadas.

Estamos acostumados a lidar com o atrito/viscosidade na mecânica clássica.

Por exemplo, a 2a lei de Newton para um paraquedista pode ser escrita como

m(dv/dt)=mg -γv

Esta é uma fórmula empírica.

Inclusive, sabe-se que em altas velocidades ela deixa de ser válida.

Do ponto de vista da engenharia, não há muito problema.

Ao desenhar um avião, usamos um programa que tem tabelado valores da viscosidade do ar para diferentes materiais e geometrias.

Mas do ponto de vista da física, como podemos entender este efeito?

Este é o objetivo da física, assim como todas as outras ciências naturais: entender.

Do micro ao macro

Para entender a viscosidade, precisamos olhar para o mundo microscópico.

Sabemos hoje que a viscosidade tem origem na interação (eletrostática) entre o nosso material e as moléculas do ar.

É como um urso sendo picado por abelhas.

load["bear2.jpg",0.3]

A picada de uma única abelha pode ser irrelevante, mas a picada de diversas abelhas produz um efeito macroscópico.

Fenômeno emergente.

Dissipação e flutuação

A viscosidade (dissipação) é na verdade um aspecto de um fenômeno mais geral.

Quando um sistema é suficientemente pequeno, os "chutes" de cada molécula podem ser suficientes para produzir deslocamentos observáveis.

Movimento browniano (video)

Descoberto por um botanista, Robert Brown, em 1827.

Oscilações incessantes de grãos de pólen em líquidos.

Ocorre com todos os materiais;

Mas em materiais macroscópicos as flutuações se tornam muito pequenas e observamos apenas a viscosidade.

O movimento Browniano é o exemplo mais famoso de um processo estocástico.

Evolução temporal com fatores aleatórios.

Sistemas quânticos abertos

Quem já estudou mecânica quântica certamente percebeu que em nenhum momento há atrito no sistema.

Descrever o atrito na mecânica quântica é muito mais difícil.

Exemplo mais simples: oscilador harmônico.

Níveis discretos de energia: (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","V(x)"}
]

A dinâmica deste sistema, quando sujeito a um atrito, ocorre através de transições entre os níveis eletrônicos.

Processo estocástico quântico.

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]

Emissão e absorção

Seja (OverBar[N]) o número médio de quanta no oscilador harmônico [lembre-se: (Subscript[E, n]) = ℏω(n+1/2)]

As taxas de transição são

Prob(n->n+1) = Δt (OverBar[N]) (n+1)
Prob(n->n-1) = Δt ((OverBar[N])+1) n

n->n+1 significa absorção.

O oscilador absorve um fóton e pula para um nível excitado.

n->n-1 significa emissão.

O sistema emite um fóton e cai para um nível inferior.

Quando n=1 -> emissão espontânea:

O sistema naturalmente decai do estado n=1 -> n=0.

Mas quando n>1 -> emissão estimulada:

Quanto maior n, maior a probabilidade de emissão.

Formalismo de Lindblad e segunda quantização

Existem diversas formas de descrever sistemas quânticos abertos.

A minha favorita é a forma de Lindblad, que é baseada na discussão anterior sobre o oscilador harmônico.

O número n que aparecia no oscilador harmônico passa a representar o número de partículas em um certo estado.

Por exemplo, podemos falar em (Subscript[n, p]) que representa “o número de partículas com momento p”

A descrição do sistema é então feita descrevendo como estes números de ocupação evoluem no tempo e como os diferentes estados interagem uns com os outros.

Essa teoria se chama segunda quantização.

Ela é a base para a física da matéria condensada e a física de partículas.

Todos os estudantes deviam aprender segunda quantização no colegial!!

O formalismo de Lindblad é uma maneira de descrever a troca de partículas e energia entre o sistema e o ambiente.

É importante diferenciar entre Bósons e Férmions:

Férmions: possuem spin semi-inteiro. Satisfazem o princípio da exclusão de Pauli

(Subscript[n, i]) = 0 ou 1

Bósons: possuem spin inteiro. Não satisfazem o PEP.

(Subscript[n, i]) = 0, 1, 2, 3,...

Detalhe técnico (não saiam correndo!):

No formalismo de Lindblad modificamos a eq. de von Neumman, adicionando um termo que descreve o acoplamento com o ambiente.

((dρ)/dt)=-i[H,ρ] + D(ρ)
D(ρ) = γ (∑_( α))[ (1± (Subscript[(OverBar[N]), α]))((Subscript[a, α]) ρ (Subscript[a, α]^†) - (1/2){(Subscript[a, α]^†)(Subscript[a, α]),ρ})+(Subscript[(OverBar[N]), α])((Subscript[a, α]^†) ρ (Subscript[a, α]) - (1/2){(Subscript[a, α])(Subscript[a, α]^†),ρ})]

Aplicações para sistemas fora do equilíbrio

Fluxo através de uma cadeia.

Considere uma cadeia 1D através da qual partículas podem pular de um sítio para o outro.

As extremidades da cadeia estão conectadas a reservatórios mantidos usando temperaturas distintas.

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

Text["(Subscript[(OverBar[N]), 1])",{0,-2.r}],
Text["(Subscript[(OverBar[N]), 2])",{xmax,-2.r}],
(*Text["(a)",{-1,0}],
*)
Opacity[0.1],Red,EdgeForm[Directive[Black,Dashed]],
Table[Rectangle[{x-r,-r},{x+r,r}],{x,{0}}],

Opacity[0.1],Blue,EdgeForm[Directive[Black,Dashed]],
Table[Rectangle[{x-r,-r},{x+r,r}],{x,{xmax}}],
Opacity[1],
Black,
PointSize[0.02],
Table[Point[{x,0}],{x,0,xmax}],
Text["g",{2.5,0.5}]
(*Text["A",{0,.75}],
Text["C",{7,.75}],
Text["B",{3.5,.75}]
*)

},BaseStyle->{FontFamily->"Times",30},ImageSize->550]

Neste caso haverá uma corrente de partículas/energia através da cadeia.

Este problema é o exemplo típico de um NESS: non-equilibrium steady-state.

É um “steady-state” pois as coisas não estão mais mudando no tempo.

Mas não é um estado de equilíbrio devido à existência de uma corrente.

A corrente de partículas através deste sistema pode ser calculada analiticamente e vale:

J = ((g γ)/((g^2)+(γ^2)))((Subscript[(OverBar[N]), 1])-(Subscript[(OverBar[N]), 2]))
g = probabilidade de tunelamento
γ = acomplamento com o banho

Ou seja, o fluxo depende apenas da diferença entre as ocupações dos sítios das bordas.

Este resultado vale tanto para Férmions quanto para Bósons.

A. Asadian, et. al., PRE, 87, 012109 (2013)
G. T. Landi, et. al., PRE, 90, 042142 (2014)

Este é um exemplo de um fluxo balístico: não depende do tamanho do sistema.

A corrente elétrica em um material não é balística; ela é difusiva. Da lei de Ohm,

I = (V/R)= (V/((Lρ/A))) ∝ (1/L)

Se você duplica o tamanho do seu fio, a corrente cai pela metade.

No nosso caso J não depende de L.

Banhos de múltiplos sítios

Considere agora o caso onde temos uma cadeia dividida em duas partes.

Cada parte é mantida em contato com um reservatório diferente.

Temperatura T

Potencial químico μ

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

Text["(Subscript[T, 1]), (Subscript[μ, 1])",{1.5,-2.r}],
Text["(Subscript[T, 2]), (Subscript[μ, 2])",{5.5,-2.r}],

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",18}]

Os números de ocupação de cada cadeia são definidos pela distribuição de Fermi-Dirac

(Subscript[(OverBar[N]), k]) = (1/((e^(((Subscript[ϵ, k])-μ)/T))+1))

Neste caso a corrente através das cadeias assume a forma de Landauer-Büttiker para a condução balística de elétrons em estruturas mesoscópicas.

Quantização da condutância.

G = ((e^2)/h)

load["voltage.pdf",0.4]

P. H. Guimarães, G. T. Landi and M. J. de Oliveira, submetido ao PRE (2016)

A corrente em função do potencial químico passa a ser muito sensível ao tamanho da cadeia.

Manipulate[
Row[{plμ[LL],plT[LL]}],{{LL,1,"L"},LLrange},ControlType->Slider]

Cadeias de spin

Um material magnético pode ser visto como vários spins alinhados na mesma direção.

Os spins se alinham devido a uma interação entre eles (chamada interação de troca).

Graphics[
Table[Arrow[{{x,0},{x,1}}],{x,0,10,1}],
ImageSize->600]

Suponha agora que você da um chute no primeiro spin.

Graphics[
Table[
Rotate[Arrow[{{x,0},{x,1}}],If[x==0,-20Degree,0]],{x,0,10,1}],
ImageSize->600]

O segundo spin vai sentir este chute, devido à interação de troca.

Consequentemente, ele irá se mover de alguma maneira.

Esse movimento será sentido pelo terceiro spin, e assim por diante.

L = 10;
h0 = 0.5;
U = -Sum[mx[i][t] mx[i+1][t] + my[i][t] my[i+1][t] + mz[i][t] mz[i+1][t],{i,1,L-1}];
Do[
hx[i] = -D[U,mx[i][t]];
hy[i] = -D[U,my[i][t]];
hz[i] = -D[U,mz[i][t]]+ h0;
,{i,1,L}]

eqs = Flatten@Table[{
D[mx[i][t],t]==my[i][t] hz[i] - mz[i][t] hy[i],
D[my[i][t],t]==mz[i][t] hx[i] - mx[i][t] hz[i],
D[mz[i][t],t]==mx[i][t] hy[i] - my[i][t] hx[i]},{i,1,L}];
ic = Join[
{mx[1][0]==Cos[20°],my[1][0]==0,mz[1][0]==-Sin[20°]},
Flatten@Table[{mx[i][0]==1,my[i][0]==0,mz[i][0]==0},{i,2,L}]];
vars = Flatten@Table[{mx[i],my[i],mz[i]},{i,1,L}];

tf = 10;
sol = Partition[NDSolveValue[{eqs,ic},vars,{t,0,tf}],3];

func[x_,tt_]:=ArcTan[(sol[[x,1]])[tt],(sol[[x,3]])[tt]];
Animate[
Graphics[
Table[
Rotate[Arrow[{{x,0},{x,1}}],func[x, tt]],{x,1,10,1}],
ImageSize->600],{tt,0,tf},

AnimationRunning->False,DefaultDuration->10]

No final teremos uma onda se propagando pelo material.

Essa é uma onda de spin.

Note que não há nenhuma partícula se propagando. Todos os spins estão parados.

Mas, mesmo assim, há a propagação de uma onda.

É como uma ola no estádio de futebol!

Ao estudarmos este fenômeno usando a mecânica quântica, descobrimos que estas onda são quantizadas.

Um quanta de onda de spin se chama um magnon.

Para todos os fins e efeitos um magnon se comporta como uma partícula.

Essa é a grande beleza da matéria condensada:

Sistemas diferentes muitas vezes satisfazem leis parecidas.

“The same equations have the same solutions”, Feynman.

Circuitos magnônicos

É possível construir dispositivos que funcionam com magnons ao invés de elétrons.

Recentemente foi proposto um transistor magnônico.

load["YIG0.png",0.5]

As correntes de magnons são geradas através de micro-antenas que produzem um campo magnético AC.

O filme de YIG é modificado de tal maneira a permitir somente a propagação de alguns tipos de magnons.

load["YIG1.png",0.5]

Os magnons vermelhos ficam aprisionados e os magnons azuis se propagam livremente.

A densidade de magnons vermelhos controla as colisões entre magnons, controlando portanto a corrente de magnons.

“All magnon transistor”

load["YIG2.png",0.6]

Porta lógica magnônica

Recentemente fomos capazes de encontrar uma solução exata para o fluxo através da cadeia de Heisenberg.

Este é um modelo que inclui a interação entre os magnons.

G. T. Landi and D. Karevski, PRB 91, 174422 (2015)

Observamos que o fluxo agora depende da constante γ (a taxa com que injetamos magnons no sistema)

load["h0.pdf",0.35]

Em seguida estudamos o que ocorre quando aplicamos campos magnéticos apenas nas extremidades.

Dependendo do sinal do campo h, existem duas opçoes:

O campo na mesma direção do polarizador.

load["drawing1.pdf",0.4]

Ou o campo em uma direção oposta ao polarizador.

load["drawing2.pdf",0.4]

Observamos que o sistema apresenta uma transição descontínua dependendo do valor do campo.

Mas só quando γ é pequeno

Row[{load["h2.pdf",0.35],load["h1.pdf",0.35]}]

Explicação: h funciona como um aprisionador de magnons, aumentando portanto as interações magnon-magnon.

Plateau representa o fluxo balístico onde os magnons não interagem.

Parte inferior representa o fluxo sub-difusivo, marcado por fortes interações magnon-magnon.

O sistema funciona portanto como uma porta lógica para magnons.

Operando o dispositivo com h próximo da transição permite chavear entre dois estados com correntes muito distintas.

É assim que funcionam as portas lógicas e transistors.

Produção de entropia

Como caracterizar um sistema fora do equilíbrio?

Mudando de assunto...

Como caracterizar um sistema fora do equilíbrio?

A chave está na entropia.

Quando um sistema sofre um processo reversível a mudança na entropia do sistema é

ΔS = ∫ (dQ/T)  := ϕ    (Processos reversíveis)

Quando o processo é feito de forma irreversível a entropia será maior que este valor:

ΔS ≥ ϕ          (Processos irreversíveis)

Ou seja, entropia foi produzida devido à irreversibilidade.

O excesso é chamado de produção de entropia:

π = ΔS - ϕ    ≥ 0                     (Produção de entropia)

Essa fórmula foi escrita no livro original do Clausius (que foi quem introduziu o conceito de entropia) em 1867.

Ou seja, a variação total da entropia tem duas contribuições:

O fluxo de entropia do ambiente para o sistema.

A entropia produzida no sistema.

ΔS = π + ϕ

Produção de entropia no NESS

Considere um sistema que é mantido constantemente fora do equilíbrio.

Então haverá uma produção constante de entropia.

Podemos escrever

((ΔS)/(Δt))=(π/(Δt))+(ϕ/(Δt))

No limite Δt -> 0 definimos então

(dS/dt) = Π + Φ
	Π = taxa de produçao de entropia ("entropia produzida por segundo").

No NESS a entropia total não vai mudar no tempo:

(dS/dt) = 0    ⟹    Π = - Φ  ≥ 0

Ou seja: toda a entropia produzida é despejada no ambiente.

Note que isso não é equilíbrio! Isso é o NESS.

Equilíbrio ⟺ Π  = 0

Exemplo: no caso do fluxo de calor entre dois banhos

Π = J ((1/(Subscript[T, 1]))-(1/(Subscript[T, 2])))

Exemplo: resistor e bateria

load["Circuit1.png",0.35]

Esse simples sistema configura um NESS:

A bateria está constantemente passando corrente pelo circuito.

E o resistor está constantemente dissipando calor para o ambiente.

Neste caso a taxa de produção de entropia no NESS será

Π = ((V^2)/(R T))  = (P/T)=((Potência)/Temperatura)

Exemplo: dois circuitos RL acoplados

O campo magnético gerado no indutor 1 afeta o indutor 2 pela indutância mútua.

load["Circuit2.png",0.8]

Portanto, há um fluxo de entropia de um circuito para o outro.

A taxa de produção de entropia neste caso é

Π = ((Subscript[V, 1]^2)/((Subscript[R, 1])(Subscript[T, 1]))) + ((Subscript[V, 2]^2)/((Subscript[R, 2])(Subscript[T, 2]))) + (((m^2) (Subscript[R, 1])(Subscript[R, 2]))/(((Subscript[L, 1])(Subscript[L, 2])-(m^2))((Subscript[L, 2])(Subscript[R, 1])+(Subscript[L, 1])(Subscript[R, 2])))) (((((Subscript[T, 1])-(Subscript[T, 2])))^2)/((Subscript[T, 1])(Subscript[T, 2])))
m = indutância mútua

Vemos que há 3 contribuições:

Um termo para cada bateria.

E um termo devido ao fluxo de calor de um circuito para o outro

Note como esse termo depende de ((((Subscript[T, 1])-(Subscript[T, 2])))^2)

G. T. Landi, T. Tomé e M. J. de Oliveira, J. Phys. A.: Math. Theor. 46 (2013) 395001

Exemplo: oscilador opto-mecânico

Uma cavidade optomecânica (interferômetro de Fabry-Perot) consiste em um espelho fixo e outro ligado a um oscilador.

O sistema é irradiado com um laser, que portanto interage com o oscilador.

Objetivo: estudar a relação entre a radiação e o meio mecânico.

load["Cavity_Optomechanics.png",0.35]

A produção de entropia no NESS desse sistema se deve à presença do laser, que constantemente fornece energia.

Existem dois canais de dissipação de energia:

O próprio canal ótico por onde o laser incide.

O banho térmico do oscilador, que dissipa energia por atrito.

A taxa de produção de entropia no NESS lê-se (não se preocupe com o significado dos parâmetros):

Π = ((4(ℰ^2))/κ) + ((γ (G^2)(ℰ^2))/((ω^2)[2n(T)+1]((Δ^2)+(κ^2))))
	Azul               = dissipação do laser (intensidade ℰ) no próprio canal ótico.
Vermelho =  dissipação do laser no banho térmico do oscilador.

No entanto, há uma terceira contribuição à produção de entropia, que advém das flutuações quânticas.

(Subscript[Σ,  quant. ]) = (((G^2)κ ((Δ^2)+(κ^2)))/(2Δ [ω((Δ^2)+(κ^2))-((ΔG)^2)]))

Estas medidas experimentais foram realizadas pelo grupo do Prof. M. Aspelmayer da Universidade de Viena.

O modelo teórico concorda muito bem com os dados experimentais.

load["exp1.png",0.4]

M. Brunelli, et. al. Submetido ao Nature Communications (2016).

Condensado de Bose-Einstein em uma cavidade ótica

Fizemos também um cálculo análogo para um condensado de Bose-Einstein aprosionado em uma cavidade ótica.

Realização experimental do modelo de Dicke.

Este sistema apresenta uma transição de fase quântica em função do acoplamento entre o campo eletromagnético e os átomos de Rb aprisionados.

Experimentos realizados pelo grupo do Prof. T. Esslinger do ETH Zürich.

load["exp2.png",0.4]

A produção de entropia reflete a transição de fase quântica, apresentando uma divergência no ponto crítico.

Teoremas de flutuação

Considere agora o seguinte problema.

Inicialmente dois sistemas são preparados a diferentes temperaturas.

Em t = 0 acoplamos os dois sistemas e deixamos eles evoluírem em contato um com o outro.

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

Text["(Subscript[T, 1])",{1.5,-2.r}],
Text["(Subscript[T, 2])",{6.5,-2.r}],

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+1,-r},{x+r+1,r}],{x,{xmax}}],
Opacity[1],
Black,
PointSize[0.02],
Table[Point[{x,0}],{x,0,3}],
Table[Point[{x,0}],{x,5,8}]

},BaseStyle->{FontFamily->"Times",18}]

Podemos nos perguntar: qual o calor que flui de um corpo para o outro após um tempo t?

Este é o problema mais simples da termodinâmica.

Sabemos que o calor vai fluir do corpo mais quente para o corpo mais frio até que os dois entre em equilíbrio térmico.

Mas isso só é verdade para sistemas macroscópicos, que é onde a termodinâmica é definida.

Se o sistema for microscópico, podemos eventualmente observar situações onde o calor flui do corpo mais frio para o mais quente.

Ocorre devido ao movimento Browniano.

Ou seja, podem ocorrer violações locais da segunda lei.

Teorema de flutuação de Jarzynski-Wójcik

Quando o calor flui entre dois sistemas microscópicos, o calor Q se torna uma variável aleatória.

Podemos falar então na probabilidade (Subscript[P, t])(Q)

(Subscript[P, t])(Q) = Prob. de que um calor Q flua entre os dois sistemas no tempo t.

Em 2004 Jarzynski e Wójcik descobrirarm que esta grandeza satisfaz a seguinte relação:

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

Ou seja,

“É exponencialmente mais provável que o calor vá fluir na direção certa”.

Note que Δβ Q representa exatamente a produção de entropia π.

Portanto, podemos dizer que (Subscript[P, t])(-Q) é a probabilidade de observar uma produção de entropia negativa.

Fluxo de calor entre duas cadeias de spin

Recentemente estudamos o fluxo de calor entre duas cadeias de spin.

(vídeos)

Conclusões

Sistemas quânticos abertos.

Apesar de não ser trivial, existem maneiras de descrever o atrito na mecânica quântica.

A minha favorita é o formalismo de Lindblad.

Sistemas fora do equilíbrio.

Descrevem diversos problemas experimentais.

Caracterizados pela existência de uma corrente no NESS.

Produção de entropia.

Melhor maneira de caracterizar os estados fora do equilíbrio.

Possuem uma contribuição intrinsecamente quântica.

Funções auxiliares

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

SetDirectory[NotebookDirectory[]];
<<"LinLib`";
<<"CustomTicks`";
load[filename_,size_]:=Show[Import[filename],ImageSize->Scaled[size]];

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

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->{{80,10},{50,10}},
FrameLabel->{"μ/t","∂F/∂μ"},
BaseStyle->{FontFamily->"Times",20},
ImageSize->300,
PlotStyle->Black,
ImagePadding->{{73,10},{50,10}},
Frame->True,
PlotStyle->Directive[Black],
Epilog->{Text["(a) L = "<>ToString@LL,Scaled[{0.8,0.92}]]}
],
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->{{80,10},{50,10}},
FrameLabel->{"μ/t","∂F/∂T"},
BaseStyle->{FontFamily->"Times",20},
ImageSize->300,
PlotStyle->Black,
ImagePadding->{{73,10},{50,10}},
Frame->True,
PlotStyle->Directive[Black],
Epilog->{Text["(a) L = "<>ToString@LL,Scaled[{0.8,0.92}]]}
],
ListLinePlot[JtlT,PlotStyle->Directive[Red,Dashed]]
]; ,{LL,LLrange}]