Second Quantization and Its Applications in Quantum Optics

Bruno Tenorio
Wolfram Research South America

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

What is this talk about?

Bosonic second quantization

Atom-field interactions

Jaynes-Cummings

Rabi oscillations

Decoherence and dissipation interactions

Evolution of states in dissipative environments

Master equations in quantum optics

Load Quantum Framework

PacletInstall["https://www.wolfr.am/DevWQCF",ForceVersionInstall->True];

<<Wolfram`QuantumFramework`
<<Wolfram`QuantumFramework`SecondQuantization`

?Wolfram`QuantumFramework`SecondQuantization`*

Truncate Fock space size

We truncate the space to favor the discrete computation capabilities of the Quantum Framework.

a := AnnihilationOperator[]

a

It’s possible to change the default size so the definitions can be called without specifying size:

SetFockSpaceSize[20];

a

Or specify the desired size  for a particular calculation:

DisplacementOperator[β,21]

Simple case examples

The annihilation operator applied to the vacuum returns the scalar 0

a @ FockState[0]

Photon distribution of a coherent state

CoherentState[][<|\[FormalAlpha]->1+2I|>]["Normalize"]["ProbabilityPlot"]

Two mode Fock Hilbert space

a["Bend"]@QuantumTensorProduct[FockState[2],FockState[5]]//TraditionalForm

Ways to calculate the mean number of photons of a state

ψ = QuantumState["RandomPure",$FockSize];

(QuantumMeasurementOperator[N[(SuperDagger[a])@a]]@ψ)["Mean"]

((SuperDagger[ψ])@ (SuperDagger[a])@a @ ψ)["Scalar"]//Chop

Dot[Range[$FockSize]-1,Diagonal[ψ["DensityMatrix"]]] //Chop

Quasi-probability representations

For Fock number state n = 5

xvec = pvec = Range[-5,5,0.05];

wigner = WignerRepresentation[FockState[5],xvec,pvec];
husimi = HusimiQRepresentation[FockState[5],xvec,pvec];

Row[{ListPlot3D[wigner,ColorFunction->"ThermometerColors",PlotLegends->Automatic,PlotRange->All,DataRange->MinMax/@{xvec,pvec},AxesLabel->{"x","p","W(x,p)"},LabelStyle->14,PlotLabel->"Wigner"],ListPlot3D[husimi,
ColorFunction->"ThermometerColors",PlotLegends->Automatic,PlotRange->All,DataRange->MinMax/@{xvec,pvec},AxesLabel->{"x","p","Q(α)"},LabelStyle->14,PlotLabel->"Husimi"]}]

Integrate over all phase space:

wigInterp= interpolateData[wigner];

Integrate[wigInterp[x,p],{x,-5,5},{p,-5,5}]

Calculate  ⟨((OverHat[x])^2)⟩ from the Wigner representation:

NIntegrate[(x^2)wigInterp[x,p],{x,-5,5},{p,-5,5}]

Atom-Field Interactions

Semiclassical approach

For two levels with energies (Subscript[E, a]) and (Subscript[E, b]), with (Subscript[E, a])=0 and letting the Bohr frequency be  (Subscript[E, b])-(Subscript[E, a]) = ℏ (Subscript[ω, 0]), the full Hamiltonian takes the form:

atomBasis = QuantumBasis[{"a","b"}];

atomicHamiltonian = (QuantumOperator["X",atomBasis](Subscript[Ω, 1]) Cos[ω t +ϕ]+ QuantumOperator[DiagonalMatrix[{0,(Subscript[ω, 0])}],atomBasis])

Considering an initial state : |ψ(t)⟩ = (Subscript[c, a])(t) |a⟩+(Subscript[c, b])(t) (e^(-i (Subscript[ω, 0 ])t))|b⟩

ψi =QuantumState[{(Subscript[c, \[FormalA]])[t],(Subscript[c, \[FormalB]])[t](E^(-I (Subscript[ω, 0]) t))},atomBasis];

Using RWA approximation we can solve Schrodinger equation

schrodingerEq = (I D[ψi,t] - (atomicHamiltonian)@ψi);
coefs =SolveValues[ Thread[schrodingerEq["AmplitudeList"]==0],
{(Subscript[c, \[FormalA]])'[t],(Subscript[c, \[FormalB]])'[t]}][[1]] //TrigToExp;
eqs = Thread[{(Subscript[c, \[FormalA]])'[t],(Subscript[c, \[FormalB]])'[t]} == coefs /.RWAδrules]

Obtaining the solution

coefSols = FullSimplify[DSolve[Join[eqs,{(Subscript[c, \[FormalA]])[0]== 1, (Subscript[c, \[FormalB]])[0] == 0}],{(Subscript[c, \[FormalA]])[t],(Subscript[c, \[FormalB]])[t]},t]] /. rabiRules

solvedState = QuantumState[Values @ coefSols[[1]],atomBasis,"Parameters"->{t,δ,Ω,(Subscript[Ω, 1]),ϕ}];

Rabi oscillations

probFuncs =solvedState[t,#,(Sqrt[1+((#)^2)]),1,0]["Probability"][(Ket["b"])]& /@ Range[0,3] ;

Grid[{{Plot[probFuncs,{t,0,6 π},],}}]

({{(Legended[, (LineLegend["δ = 0", "δ = \!\(\*SubscriptBox[\(Ω\), \(1\)]\)", "δ = 2\!\(\*SubscriptBox[\(Ω\), \(1\)]\)", "δ = 3\!\(\*SubscriptBox[\(Ω\), \(1\)]\)"])]), }})

Jaynes–Cummings Hamiltonian

SetFockSpaceSize[32];
atomBasis = QuantumBasis[{"e","g"}];
σplus =  QuantumOperator["JZ-",atomBasis];
σminus = QuantumOperator["JZ+",atomBasis];

(SuperMinus[σ])= |g⟩ ⟨e|,  (SuperPlus[σ])= | e⟩ ⟨ g|

Free Hamiltonian

H0 = QuantumOperator[
h ν QuantumOperator[(SuperDagger[a])@a,{2} ]+ 
1/2 h ω QuantumOperator["Z",atomBasis],
"Parameters"->{ν,ω}];

Atom-field interaction term

H1 = h g (QuantumTensorProduct[σplus,a]+
QuantumTensorProduct[σminus,(SuperDagger[a])]);

Interaction picture Hamiltonian (at resonance)

HInteract = Simplify[Exp[I H0/h t]@ H1 @ Exp[-I H0/h t]];
U = Exp[-I HInteract[ν,ν]/h t];

Case with initial field in Fock state |1⟩

Rabi oscillations when the field is initially in a Fock state

ψ0 = QuantumTensorProduct[QuantumState["1",atomBasis],FockState[1]]
ψt = U @ ψ0

proj =QuantumOperator[{{1,0},{0,0}},atomBasis];
probE =(QuantumMeasurementOperator[proj]@ ψt)["Mean"];

Plot[probE /. g->1 ,{t,0,3π},]

Case with initial field in a coherent state

Initial state   with |α| = 4

α = 4Exp[I RandomReal[{0,2π}]];

ψ0=QuantumTensorProduct[QuantumState["1",atomBasis],CoherentState[][α]["Normalize"]];

ψt = QuantumState[U @ ψ0,"Parameters"->{g,t}];

{cet[t_],cgt[t_]} =Values[(#["Dagger"]@ψt[1,t])["Amplitudes"]]& /@ atomBasis["BasisStates"];

Collapse and revivals

Row[{Plot[Total @(Abs[cet[t]]^2-Abs[cgt[t]]^2),
{t,0,40},],}]

(RowDefault)

Evolution of the photon number distribution

Animate[ListLinePlot[((Abs[cet[t]])^2)+((Abs[cgt[t]])^2),],{t,0,30,0.1`},]

Case with the initial field in a thermal state
(density matrix approach)

Atom initially in  the superposition (1/(Sqrt[2]))(|e⟩+|g⟩) and field in a thermal state with (OverBar[n])=5

ρ0 = QuantumTensorProduct[QuantumState["+",atomBasis]["MatrixState"],
ThermalState[5.]];

Evolution of the density matrix  ρ(t) =(OverHat[U])(t) ρ(0) (SuperDagger[(OverHat[U])])(t)

ρt  = QuantumState[(SuperDagger[U])[ρ0],"Parameters"->{g,t}];

atomMatrix = Normal[QuantumPartialTrace[ρt[1,t],{2}]["DensityMatrix"]]//ComplexExpand//Simplify;

Plot[Evaluate[atomMatrix[[1,{1,2}]]],{t,0,200},]

Atom gets “thermalized”, time average of the coherences → 0

timeAverages = Integrate[(atomMatrix[[1,2]]),{t,0,tf}]/tf;

Plot[timeAverages,{tf,5,200},]

Dissipation , decoherence and master equations

Decay of a coherent state

SetFockSpaceSize[15];
α = 2 Exp[I RandomReal[{0,2Pi}]];
ρ0 = coherent[α];
γ = RandomReal[{1,2}];

Considering standard evolution with dissipation

ℋ = QuantumOperator["Hamiltonian"[(SuperDagger[a])@a,{a},{γ}]];

ρt =QuantumState[Exp[-I ℋ t][ρ0],"Parameters"->{t}];

time = Range[0,1,0.01];
fock =Range[$FockSize]-1;
qs = QuantumSimilarity[coherent[α (E^(-I #)) (E^(-(1/2) # γ))],ρt[#]]& /@ time;
meanN = Chop[fock.Diagonal[ρt[#]["DensityMatrix"]]]& /@ time ;

With photon loss mechanism the state remains coherent with amplitude α(t) = α(0) (e^(- γ t /2)) (e^(- i ω t)), so the mean number of photons decay

Labeled[Grid[{{ListLinePlot[Thread[{time,meanN}],],ListLinePlot[Thread[{time,qs}],]}}],]

Visualizing the Wigner representation (spiral motion)

Generation of coherent states from decoherence: optical balance

G = RandomReal[{1,2}];

ℋ= QuantumOperator["Hamiltonian"[G(a+a["Dagger"]),{a},{γ}]];

ρ0 = QuantumState["RandomPure",$FockSize];

Steady state solution is a coherent state |-2 i (G/γ)⟩⟨-2 i (G/γ)|

expℋ =QuantumEvolve[ ℋ,None,{t,0,3}];

ρt := QuantumState[expℋ[t][ρ0],"Parameters"->t];

fidelities= Table[{t,QuantumSimilarity[ρt[t],coherent[-I 2G/γ],"Fidelity"]},{t,0,3,0.1}];

ListLinePlot[fidelities,]

For initial (OverHat[ρ])(0) = |0⟩ ⟨0| evolves to coherent state | -i G t ⟩⟨- i G t | for short times, where the complex amplitude is independent of γ

ρ0 = QuantumState["0",$FockSize];

time = Range[0,0.8,0.02];
qs = QuantumSimilarity[coherent[-I G #],ρt[#]]& /@ time;
meanN = Chop[fock.Diagonal[ρt[#]["DensityMatrix"]]]& /@ time ;

ListLinePlot[{Thread[{time,qs}],Thread[{time,meanN}]},]

Kerr Interaction

We explore  two photon driven KNR Hamiltonian in a frame rotating at resonator frequency to prepare and stabilize states in a manifold of coherent states

where (Subscript[ℒ, loss])(ρ) is the jump operator corresponding (OverHat[a])

Steady state solution driven KNR

SetFockSpaceSize[12];

(* Based on Puri et.al.(2017) *)
Ka = -1.;     (* Kerr non-linear strength*)
(Subscript[e, p]) = -4.;     (* Drive strength*)
κ = 8. ;        (* Photon loss*)

ρ0 = FockState[0];

H0 = -Ka ((SuperDagger[a])@(SuperDagger[a])-(Subscript[e, p])\[Conjugate]/Ka )@(a@a-(Subscript[e, p])/Ka )+((Abs[(Subscript[e, p])])^2)/Ka ;

ℋ = QuantumOperator["Hamiltonian"[H0,{a},{κ}]];
expℋ= QuantumEvolve[ℋ,None,{t,0,2}];

Steady state (1/2)(|(Subscript[α, 0])⟩ ⟨(Subscript[α, 0])|+|-(Subscript[α, 0])⟩ ⟨-(Subscript[α, 0])|)

r0=((((((4 ((Subscript[e, p ])^2)-(κ^2)/4))/(4 (Ka^2)))))^(1/4)); θ0=ArcTan[(κ/(Sqrt[16 ((Subscript[e, p])^2)-(κ^2)]))]/2;

ρs=1/2coherent[#]["MatrixState"]&/@{r0 (E^(I θ0)),-r0 (E^(I θ0))} //Total;

QuantumSimilarity[expℋ[2.][ρ0],ρs,"Fidelity"]

Transforming the ground state to the coherent manifold

Adiabatic initialization of cat state

We use a time dependent two photon drive: (Subscript[ϵ, p])(t)=(Subscript[ϵ, 0])(1-exp((t^4)/(τ^4)))

(*Based on Puri et.al.(2017)*)
Ka = 4.;                                                        (* Kerr non linear strength*)
ep0 = 4 Ka;                                                  (* Drive strength *)
(Subscript[e, p]) = ep0(1-E^(-t^4/tau^4));     (* Time dependent drive*)
κ = 0.;                                                          (* Photon loss*)
tau = 5/Ka;                                                (* Adiabatic parameter*)

H0 =-Ka ((SuperDagger[a])@(SuperDagger[a])-(Subscript[e, p])\[Conjugate]/Ka )@(a@a-(Subscript[e, p])/Ka )+ ((Abs[(Subscript[e, p])])^2)/Ka ;

|0⟩ evolves to  even cat state and |1⟩ to odd cat state, the two degenerate eigenstates of KNR without driving:

ℋ = QuantumOperator["Hamiltonian"[H0]];

expℋ= QuantumEvolve[ℋ, None,{t,0,2}];

αt[tx_]:=(((((4 ((((Subscript[e, p]) /. t->tx))^2)-(κ^2)/4)/(4 (Ka^2)))))^(1/4))Exp[I ArcTan[(κ/(Sqrt[16 ((((Subscript[e, p]) /. t->tx))^2)-(8^2)]))]/2]

Computing similarities respect to even and odd cat state:

QuantumSimilarity[CatState[][αt[2.],0 ],expℋ[2.][FockState[0]]]

QuantumSimilarity[CatState[][αt[2.],Pi],expℋ[2.][FockState[1]]]

Wigner representation adiabatic evolution

frames =Table[Legended[Show[ListPlot3D[
WignerRepresentation[expℒ[t][ρ0],xvec,pvec],
],],BarLegend[]],{t,0,2,0.1}];

Jaynes- Cummings with dissipation

(Subscript[ℒ, decay, ])(Subscript[ℒ, loss]) are the jump operators corresponding to (SuperMinus[(OverHat[σ])])and (OverHat[a])

Considering the  initial state (Ket[e])(Ket[n]),  with n=2

SetFockSpaceSize[8];

ν = 10.;   (* field frequency *)
ω = 10.;   (* atomic frequency *)
gc = 15.; (* interaction coupling strength *)

α = 2 Exp[I RandomReal[{0,2π}]];
field0 = FockState[2];
atom0 = QuantumState["0",atomBasis];

H =  ν QuantumOperator[(SuperDagger[a])@a,{2} ]+ 
1/2 ω QuantumOperator["Z",atomBasis]+
gc (QuantumTensorProduct[σplus,a]
+QuantumTensorProduct[σminus,(SuperDagger[a])]);

ψ0 = QuantumTensorProduct[atom0,field0];

decayAtom = QuantumTensorProduct[σminus,QuantumOperator["I"[$FockSize]]];
lossField = QuantumTensorProduct[QuantumOperator["I",atomBasis],a];

Damped Rabi oscillations

atomProb[evolved_]:=Table[{t,QuantumPartialTrace[evolved[t],{2}]["DensityMatrix"][[1,1]]//Chop},{t,0,1,0.01}]

probsE =Table[ℋ = QuantumOperator["Hamiltonian"[H,{decayAtom,lossField},{x,x}]];
expℋ = QuantumEvolve[ℋ,None,{t,0,1}];evolved = QuantumState[expℋ[t][ψ0],"Parameters"->t];
atomProb[evolved],
{x,{0.8,0.4,0.2}}];

ListLinePlot[probsE,]

(Legended[, (LineLegend["γ,κ = 0.8", "γ,κ = 0.4", "γ,κ = 0.2"])])

Grid[{{%134,%138},{%142,%145}}]

Labeled[%149,Style["Evolutions for different states for γ=κ=0.8",FontFamily->"Arial"],Top]

Summary

An implementation of second quantization was integrated in the Quantum Framework, showcasing versatility of the quantum phenomena that can be simulated.

The tool can be useful to develop educational content for courses in introductory quantum optics and quantum mechanics.

Research-oriented simulations can be realized in problems that can be modelled with master equations.

Thanks!

Bruno Tenorio
brunot@wolfram.com

References

Second Quantization Tech Note https://resources.wolframcloud.com/PacletRepository/resources/Wolfram/QuantumFramework/tutorial/SecondQuantization.html

Gerry, C. C., & Knight, P. L. (2005). Introductory quantum optics. Cambridge university press.

Grynberg, G., Aspect, A., & Fabre, C. (2010). Introduction to quantum optics: from the semi-classical approach to quantized light. Cambridge university press.

Scully, M. O., & Zubairy, M. S. (1997). Quantum optics. Cambridge university press.

Puri, S., Boutin, S., & Blais, A. (2017). Engineering the quantum states of light in a Kerr-nonlinear resonator by two-photon driving. npj Quantum Information, 3(1), 18.
https://www.nature.com/articles/s41534-017-0019-1

PacletInstall["https://www.wolfr.am/DevWQCF",ForceVersionInstall->True];

interpolateData[data_]:= Interpolation[Flatten[MapThread[,Join[,{data}],2],1]];

$Assumptions = {t>0, δ>0 , (Subscript[Ω, 1])>0,ω ∈ Reals, h ∈ Reals,ν ∈ Reals};

RWAδrules ={(E^(I ϕ + I t ω + I t (Subscript[ω, 0])))->0,(E^(-I ϕ - I t ω - I t (Subscript[ω, 0])))->0,(I t ω-I t (Subscript[ω, 0]))-> I t δ,(I t (Subscript[ω, 0])-I t ω)-> -I t δ};

rabiRules = {(Sqrt[_])->Ω,1/(Sqrt[_])->1/ Ω};

coherent := CoherentState[]["Normalize"];