Advanced Applications of Bosonic Second Quantization in Wolfram Language

Bruno Tenorio, WRSA

#WolframTechConf

ABSTRACT

In this talk, we discuss the formalism of bosonic second quantization. Building on recent developments in the Wolfram Quantum Framework and Wolfram Language, we explore advanced applications in quantum optics. We show how interference phenomena of light can be described in terms of “bright” and “dark” states, offering a fully quantum perspective. Symbolic methods are presented for handling the non-commuting algebra of bosonic operators, enabling alternative approaches to calculations requiring operator simplification. We also examine multi-mode phase-space representations, focusing on the joint Wigner function for two-mode systems. The emphasis is on making advanced quantum ideas more approachable through computation, bridging abstract theory and practical exploration in both learning and research.

What is this talk about?

Bright and dark states of light

Non commutative algebra in quantum optics

Two mode phase space representation: joint Wigner function

Load Wolfram Quantum Framework

Install the paclet:

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

Load the Quantum Framework:

Needs["Wolfram`QuantumFramework`"]

Load bosonic second quantization tools:

Needs["Wolfram`QuantumFramework`SecondQuantization`"]

?Wolfram`QuantumFramework`SecondQuantization`*

Bright and dark states of light

Light-Matter interaction

Energy exchange Hamiltonian between field and atom sensor, where the field operators (E^((+))), (E^((-))) are related to positive and negative frequency components

Probability of a photon from a single mode detected by a sensor atom:

Set truncation size and define a single mode operator:

SetFockSpaceSize[5];
a := AnnihilationOperator[];

One can show that for a generic state (Ket[ψ])=(∑_n)(Subscript[c, n])(Ket[n]), the mean value of number operator is zero, ⟨ψ|(SuperDagger[a])a|ψ⟩=0, only if ∀n>=1, (Subscript[c, n])=0

Reduce[With[{ψ=a@QuantumState[Array[c,$FockSize,0],$FockSize]},((SuperDagger[ψ])@ψ)["Scalar"]]==0]

Evidently the only state that is undetectable is the vacuum:

(SuperDagger[FockState[0]])@((SuperDagger[a])@a)@FockState[0]

Now let’s consider two modes, we will get more interesting results than the single mode case. Each mode has an annihilation and creation operator (a, b) , and including a relative phase we can write

Second mode operator:

b:= AnnihilationOperator[{2}]

Find two instances of undetectable states:

FindInstance[With[{ψ=(a+b)@QuantumState[Flatten@Array[c,{$FockSize,$FockSize},{0,0}]]},((SuperDagger[ψ])@ψ)["Scalar"]]==0,Flatten@Array[c,{$FockSize,$FockSize},{0,0}],2]//Transpose//TableForm

In fact, there are infinite states that are undetectable such that H |g⟩ |Ψ⟩ = 0. Really different from the single mode case!

Collective operators and states

Consider a larger truncation size of Fock space:

SetFockSpaceSize[10];

Symmetric and antisymmetric collective operators:

c :=(a+b (E^(I θ)))/(Sqrt[2]);
d := (-a (E^(-I θ))+b)/(Sqrt[2]);

The following commutators hold:

Commutator[a,b]

Commutator[c,d]

Photon preserving property,

randomState = QuantumState["RandomPure",{$FockSize,$FockSize}];

Chop[((SuperDagger[randomState])@((SuperDagger[a])@a+(SuperDagger[b])@b)@randomState)["Scalar"]]

Chop[Simplify[((SuperDagger[randomState])@((SuperDagger[c])@c+(SuperDagger[d])@d)@randomState)["Scalar"]]]

Collective states defined for 0<= n <= N, fixing some number of maximum photons N:

One can show that an equivalent definition is:

We will use the new definition of  since in the previous one, the matrix exponentials/powers have high computational cost:

ClearAll[ψNn]
ψNn[𝒩_,n_ ]/; 0<=n<=𝒩<$FockSize:=
QuantumState[Sum[(((E^(I θ (k-n)))(Sqrt[k!(𝒩-k)!]))/((Sqrt[n!(𝒩-n)!(2^𝒩)])))Sum[Binomial[n,j]Binomial[𝒩-n,k-j](((-1))^(k-j)),{j,Max[0,k-(𝒩-n)],Min[n,k]}]FockState[{k,𝒩-k}],{k,0,𝒩}],"Parameters"->θ]

Try some example, 𝒩 = 8, n = 4:

ψNn[8,4]

Properties of collective states

For a given total number of photons N, the states (Ket[((Subscript[ψ, n])^N)]) form a complete basis

Table[Simplify@Sum[ψNn[n,i]["Projector"],{i,0,n}],{n,$FockSize-1}]//Column//TraditionalForm

Connection of collective state (N,n) and (N-1, n-1) by the operator (OverHat[c])

Table[(c[ψNn[Np,j]]-Sqrt[j] ψNn[Np-1,j-1])["Norm"],{Np,1,$FockSize-1},{j,1,Np}]

Mean value of number operator for mode a  :

And@@@Table[Simplify[With[{ψ=a[ψNn[Np,j]]},ψ["Norm"]^2==Np/2]],{Np,1,$FockSize-1},{j,0,Np}]

This also holds for mode b

Defining the bright and dark states

A  perfectly dark state (PDS), or simply dark state, are states generated with n=0 such  that they cannot excite the sensor atom

Checking the amplitudes of the analytical expression hold (“AmplitudeList” returns non zero amplitudes):

Table[ψNn[Np,0]["AmplitudeList"]==Table[(Sqrt[((Np!)/(2^Np))])(((((-1))^m)(E^(I m θ)))/(Sqrt[m!(Np-m)!])),{m,0,Np}],{Np,$FockSize-1}]

(E^( (+)))(r,t) is also proportional to the collective operator (OverHat[c]) , dark states must be eigenstates of (OverHat[c]) with eigenvalue 0:

Table[c[ψNn[Np,0]]//Simplify//TraditionalForm,{Np,$FockSize-1}]

This means they are not detectable

On the other hand, the states generated when n=N were coined as super-radiant states in the context of atom decay. They are referred as maximally super-radiant states (MSS) or just bright states

Comparing amplitudes:

Table[ψNn[Np,Np]["AmplitudeList"]==Table[(E^(-I Np θ))(Sqrt[((Np!)/(2^Np))])((E^(I m θ))/(Sqrt[m!(Np-m)!])),{m,0,Np}],{Np,$FockSize-1}]

The bright states (MSS)  have a transition rate  (Sqrt[2N]),  that is higher than the single-mode case

Atomic part definitions:

atomBasis = QuantumBasis[{"e","g"}];
σplus = QuantumOperator["+",atomBasis];
σminus = QuantumOperator["-",atomBasis];
{excited,ground}= atomBasis["BasisStates"];

Hamiltonian for two modes involving the interaction part:

H =  QuantumTensorProduct[σplus,a + b (E^(I θ))]+QuantumTensorProduct[σminus,(SuperDagger[a])+(SuperDagger[b])(E^(- I θ))];

Indeed , stronger coupling compared to the single mode:

Table[Normal[H[ψNn[Np,Np][ground]]["StateVector"]==
(Sqrt[2Np])ψNn[Np-1,Np-1][excited]["StateVector"]],{Np,1,$FockSize-1}]

State “ladder” of collective states (figure taken from the reference paper):

Explaining double-slit: Interference from coherence, not from a “wavefront”

What we see and measure can  be derived purely from quantum correlation functions and the process of photon detection

Double slit diagram:

The probability for a sensor atom ( to a “click” detector) to absorb a photon at a position  is proportional to the first-order correlation

with  the state prepared at the double-slit.

Single photon case

For the simple case of one-photon, one gets:

Assuming  |(Subscript[k, 1])|=|(Subscript[k, 2])|the time dependent part can be factored and discarded (global phase). If we define ϕ = (Subscript[k, 2]).(Subscript[d, 2])-(Subscript[k, 1]).(Subscript[d, 1]) we get:

Also θ = ((Subscript[k, 2])-(Subscript[k, 1])). r, so we have

Therefore, the first-order correlation reads as:

Collective basis using :

qb = QuantumBasis[Association[Table[Superscript[Subscript["ψ",j],i]->If[j<=i,ψNn[i,j]["StateVector"],
SparseArray[{},{($FockSize^2)}]],{i,0,$FockSize-1},{j,0,$FockSize-1}]//Flatten]];

Define the state |S⟩:

Sstate = (1/(Sqrt[2]))(FockState[{1,0}]+(E^(-I ϕ))FockState[{0,1}])

Change of basis:

QuantumState[Sstate,qb]

With an extra global it can take a nicer form. Depending on the value of δ = θ - ϕ, we can get only a bright state, a dark state or a superposition:

FullSimplify[(E^(-I ((θ-ϕ)/2)))%]

Only the bright state contributes to the first-order correlation

Coherent state case

SetFockSpaceSize[10];

Collective basis of  states, using θ = 0:

qb:=QuantumBasis[Association[Flatten[Table[(Superscript[(Subscript["ψ", j]), i])->If[j≤i,ψNn[i,j][0]["StateVector"],SparseArray[{},{($FockSize^2)}]],{i,0,$FockSize-1},{j,0,$FockSize-1}]]]];

It is possible to show in-phase coherent states decompose only on the MSS subspace (considering θ = 0)

Two mode not normalized coherent state with arbitrary amplitudes:

twoModecoherent := QuantumState[QuantumTensorProduct[CoherentState[][α1],CoherentState[][α2]],"Parameters"->{α1,α2}];

Decomposing the state using a change of basis to collective states, we get only bright states:

MSSdecomposition = QuantumState[twoModecoherent[α,α],qb]//Simplify

, this implies enhanced absorption by a factor of 2 compared to the single mode case. 
Exhibits constructive interference

SetFockSpaceSize[20];

H = QuantumOperator[ QuantumTensorProduct[σplus,a + b (E^(I θ))]+ QuantumTensorProduct[σminus,(SuperDagger[a])+(SuperDagger[b])(E^(- I θ))],"Parameters"->θ];

Numerically the error is  really small up to α = 2.5, for this truncation size of 20:

With[{sim=(H[0][twoModecoherent[x,x]["Normalize"]@ground]-2x twoModecoherent[x,x]["Normalize"]@excited)["Norm"]},
Plot[sim,{x,0,2.5},]]

And two coherent states with opposite phases decompose only in the PDS subspace

SetFockSpaceSize[10];

Decomposing in the collective states basis, we get only dark states:

PDSdecomposition = QuantumState[twoModecoherent[α,-α],qb]//Simplify

This can be interpreted as destructive interference for two fields with same amplitude but opposite phases

Naturally this superposition of dark states  has to be undetectable

In the context of the two slit experiment used with coherent states, can be described by a two mode coherent state with the same |α| but a relative phase dependent on the geometry of the slits and detection point (see diagram in single photon case)

In general in the collective basis this state has all components:

QuantumState[twoModecoherent[α,α (E^(I β))],qb]["Amplitude"]//Keys

Only when β = 2 l π or β = (2l+1) π, with l an integer we get decomposition only of MSS or PDS respectively

Non commutative algebra: bosonic algebra

Heisenberg picture: commutators and evolution of operators

Differential equation for the operator in Heisenberg picture:

Single mode annihilation and creation operators

Define scalar variables:

alg = NonCommutativeAlgebra[<|"ScalarVariables"->{h ,ω,m}|>];

Hamiltonian for a single mode of radiation  ℏω ((SuperDagger[a])a+1/2):

Hsinglemode = h ω((SuperDagger[\[FormalA]])**\[FormalA]+1/2);

Annihilation and creation operators hold the commutation relation [a, (SuperDagger[a])] = 1 :

commRel = {Commutator[\[FormalA],(SuperDagger[\[FormalA]])]-1};

Calculate (i/h) [H,a] :

(I/h) NonCommutativePolynomialReduce[Commutator[Hsinglemode,\[FormalA]],commRel,{\[FormalA],(SuperDagger[\[FormalA]])},alg][[-1]]

Calculate (i/h) [H,(SuperDagger[a])] :

(I/h) NonCommutativePolynomialReduce[Commutator[Hsinglemode,(SuperDagger[\[FormalA]])],commRel,{\[FormalA],(SuperDagger[\[FormalA]])},alg][[-1]]

Both cases have exponential solutions, operators at time 0 coincide with the corresponding Schrodinger picture ones

And trivially the number operator commutes with the Hamiltonian:

(I/h) NonCommutativePolynomialReduce[Commutator[Hsinglemode,(SuperDagger[\[FormalA]])**\[FormalA]],commRel,{\[FormalA],(SuperDagger[\[FormalA]])},alg][[-1]]

Position and momentum in quantum harmonic oscillator

Define the algebra and QHO Hamiltonian H  = ((p^2)/(2m))+(1/2)m (ω^2)(x^2) :

alg = NonCommutativeAlgebra[<|"ScalarVariables"->{m,ω, h}|>];
Hharmonic = ((p**p)/(2 m))+((m (ω^2) x**x)/2);

Canonical commutator [x, p] = i ℏ :

canonicalCommutator = {Commutator[x,p]- I h};

Calculate i /ℏ  [H , x] :

dtx = (I/h)NonCommutativePolynomialReduce[Commutator[Hharmonic,x],canonicalCommutator,{x,p},alg][[-1]]

Calculate i / ℏ  [H , p] :

dtp = (I/h)NonCommutativePolynomialReduce[Commutator[Hharmonic,p],canonicalCommutator,{x,p},alg][[-1]]

So now we can solve the associated differential equations:

DSolve[
{D[x[t],t]==dtx/.p->p[t],
D[p[t],t]==dtp/.x->x[t],
x[0]==(Subscript[\[FormalX], 0]),
p[0]==(Subscript[\[FormalP], 0])},
{x[t],p[t]},t][[1]]//Simplify//Column

Normal ordering

Define a power of (a + (SuperDagger[a])):

op =GeneralizedPower[NonCommutativeMultiply,(\[FormalA]+(SuperDagger[\[FormalA]])),3]

Expansion is not normal ordered:

NonCommutativeExpand[op,alg]

Perform normal ordering:

NonCommutativePolynomialReduce[op,commRel,{\[FormalA],(SuperDagger[\[FormalA]])}][[-1]]

Non commutative power (tweaking for automatic simplification for zeroth power):

genPower[x_,k_]:=If[k==0,1,GeneralizedPower[NonCommutativeMultiply,x,k]]

Define a function to calculate [((SuperDagger[a])^(Subscript[n, 1]))(a^(Subscript[m, 1])),((SuperDagger[a])^(Subscript[n, 2]))(a^(Subscript[m, 2]))] in the normal order:

ClearAll[com];
com[{n1_,m1_},{n2_,m2_}]:=NonCommutativePolynomialReduce[Commutator@@NonCommutativeMultiply@@@Partition[genPower@@@{{(SuperDagger[\[FormalA]]),n1},{\[FormalA],m1},{(SuperDagger[\[FormalA]]),n2},{\[FormalA],m2}},2],commRel,{\[FormalA],(SuperDagger[\[FormalA]])}]

Manually the calculation has this form

Define the formula:

ClearAll[comManual];
comManual[{n1_,m1_},{n2_,m2_}]:=
Sum[
k!(Binomial[m1,k]Binomial[n2,k]-Binomial[m2,k]Binomial[n1,k])NonCommutativeMultiply@@genPower@@@{{(SuperDagger[\[FormalA]]),n1+n2-k},{\[FormalA],m1+m2-k}},{k,Max[Min[n1,m2],Min[n2,m1]]}]//NonCommutativeExpand

Test for random integers:

With[{c=RandomInteger[{0,10},{2,2}]},Print["c = " <>ToString[c]];{comManual@@c == Last[com@@c],comManual@@c}]

Displacement Operator

Resource functions for calculating the nth degree term of Baker Campbell Hausdorff and Zassenhaus expansions:

BCH = ResourceFunction["BakerCampbellHausdorffTerms"]

Zassenhaus = ResourceFunction["ZassenhausTerms"]

Define an algebra with scalar variables will be used:

alg = NonCommutativeAlgebra[<|"ScalarVariables"->{α,β,ϵ}|>];

Now we want to manipulate the displacement operator 
we write relations in a redundant manner on purpose:

relations = {Commutator[-α\[Conjugate] \[FormalA],Commutator[-α\[Conjugate]\[FormalA],α (SuperDagger[\[FormalA]])]],Commutator[α (SuperDagger[\[FormalA]]),Commutator[-α\[Conjugate]\[FormalA],α (SuperDagger[\[FormalA]])]],Commutator[-α\[Conjugate]\[FormalA],α (SuperDagger[\[FormalA]])]+((Abs[α])^2)};

Groebner basis indicates the only necessary relation is fundamental commutator, which is expected:

gb = NonCommutativeGroebnerBasis[relations,{\[FormalA],(SuperDagger[\[FormalA]])},alg]//FullSimplify

Zassenhaus simplification for displacement operator, normal ordered form has this form :

Sum[NonCommutativePolynomialReduce[Zassenhaus[{α (SuperDagger[\[FormalA]]),-α\[Conjugate]\[FormalA]},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg],{n,2,10}][[-1]]

Anti-normal order  , we write the summands in the opposite order when calling Zassenhaus:

Sum[NonCommutativePolynomialReduce[Zassenhaus[{-α\[Conjugate]\[FormalA],α (SuperDagger[\[FormalA]])},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg],{n,2,10}][[-1]]

Or we can go the other way around using the BCH formula:

(1/2)((Abs[α])^2)+Sum[NonCommutativePolynomialReduce[BCH[{-α\[Conjugate]\[FormalA],α (SuperDagger[\[FormalA]])},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg],{n,1,10}][[-1]]//FullSimplify

-(1/2)((Abs[α])^2)+Sum[NonCommutativePolynomialReduce[BCH[{α (SuperDagger[\[FormalA]]),-α\[Conjugate]\[FormalA]},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg],{n,1,10}][[-1]]//FullSimplify

Making the central operator also an exponential

D[Exp[Sum[NonCommutativePolynomialReduce[BCH[{\[FormalA] Conjugate[α]-α (SuperDagger[\[FormalA]]), ϵ \[FormalA],α (SuperDagger[\[FormalA]])-α\[Conjugate]\[FormalA]},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg][[-1]],{n,1,6}]],ϵ]/.ϵ->0

Similarly:

D[Exp[Sum[NonCommutativePolynomialReduce[BCH[{\[FormalA] Conjugate[α]-α (SuperDagger[\[FormalA]]), ϵ (SuperDagger[\[FormalA]]),α (SuperDagger[\[FormalA]])-α\[Conjugate]\[FormalA]},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg][[-1]],{n,1,6}]],ϵ]/.ϵ->0

We get a phase factor and the displacement operator term z (SuperDagger[a])-(SuperStar[z])a with z = α + β:

Sum[NonCommutativePolynomialReduce[BCH[{α (SuperDagger[\[FormalA]])-α\[Conjugate]\[FormalA], β (SuperDagger[\[FormalA]])-β\[Conjugate]\[FormalA]},n,NonCommutativeMultiply],gb,{(SuperDagger[\[FormalA]]),\[FormalA]},alg][[-1]],{n,1,7}]

Joint Wigner Function

Mathematical Definition

One way to define the joint Wigner function for two modes 𝑎, 𝑏 is through a trace expression that depends on the displacement and phase-shift (parity) operators.  Explicitly, it can be written as:

For convenience, this expression can be rewritten in terms of a kernel 𝒦 defined as:
 with

And the associated matrix elements in the Fock space are:

Implementation

𝒦 kernel implementation:

ClearAll[𝒦mn];
𝒦mn[γ_,{m_,n_}]:=With[{ζ=2γ},(((-1))^n)Exp[-(((Abs[ζ])^2)/2)]If[m>=n,(Sqrt[((n!)/(m!))])If[m==n,1,(ζ^(m-n))]LaguerreL[n,m-n,((Abs[ζ])^2)],(Sqrt[((m!)/(n!))])If[m==n,1,((Conjugate[-ζ])^(n-m))]LaguerreL[m,n-m,((Abs[ζ])^2)]]]

Numeric compiled version of (Subscript[𝒦, mn]):

Kmncomp=FunctionCompile[Function[{Typed[γ,"ComplexReal64"],Typed[m,"Integer64"],Typed[n,"Integer64"]},Module[{ζ,res},ζ=2.0*γ;
res=If[m>=n,Sqrt[Gamma[n+1.0]/Gamma[m+1.0]]*If[m==n,1.+0.I,ζ^(m-n)]*LaguerreL[n,m-n,Abs[ζ]^2],Sqrt[Gamma[m+1.0]/Gamma[n+1.0]]*If[m==n,1.+0.I,Conjugate[-ζ]^(n-m)]*LaguerreL[m,n-m,Abs[ζ]^2]];
(-1)^n*Exp[-Abs[ζ]^2/2.0]*res]]];

Definitions without explicit trace (speed will depend on how sparse the density matrix is):

ClearAll[NWigner2mod]
NWigner2mod[rho_SparseArray,α_,β_]:=
With[
{d=√First[rho["Dimensions"]],vals=rho["ExplicitValues"],pos=rho["ExplicitPositions"]},
(4./(π^2)) Dot[vals,Kmncomp[α,Quotient[#2-1,d],Quotient[#1-1,d]]
			Kmncomp[β,Mod[#2-1,d],Mod[#1-1,d]]
				&@@@pos]]

ClearAll[Wigner2mod]
Wigner2mod[rho_SparseArray,α_,β_]:=
With[
{d=√First[rho["Dimensions"]],vals=rho["ExplicitValues"],pos=rho["ExplicitPositions"]},
(4/(π^2)) Dot[vals,𝒦mn[α,{Quotient[#2-1,d],Quotient[#1-1,d]}]
			𝒦mn[β,{Mod[#2-1,d],Mod[#1-1,d]}]
				&@@@pos]]

Properties

Global inversion symmetry (always true):

Partial transpose of the bipartite system does the transformation to the Wigner function:

Parity symmetry (not general)

For a separable state ρ = (Subscript[ρ, A]) ⊗ (Subscript[ρ, B]), the joint Wigner function factorizes:

Fock states |m,n⟩

Density matrices for states |m,0⟩:

SetFockSpaceSize[8];
ρs=Array[FockState[{#,0}]["DensityMatrix"]&,6];

Spatial range:

range = N@Subdivide[-2,2,20];

Sample points for 3D slices((Subscript[p, 2])=0):

sampsFock=Table[Flatten[Table[{x1,p1,x2,NWigner2mod[ρs[[i]],((x1+I p1)/(Sqrt[2])),(x2/(Sqrt[2]))]},{x1,range},{p1,range},{x2,range}],2],{i,6}];

Show the plots:

Grid@Partition[Table[ListDensityPlot3D[sampsFock[[n]],],{n,6}],2]

Fock state |6,2⟩:

Check symmetry properties for state |2,3⟩:

ψ = FockState[{2,3}];

Wigner function:

wig = Function[{x1,p1,x2,p2},Evaluate@Simplify[Wigner2mod[ψ["DensityMatrix"],((x1+I p1)/(Sqrt[2])),((x2+I p2)/(Sqrt[2]))],Element[{x1,p1,x2,p2},Reals]]];

Global symmetric:

Simplify[wig[x1,p1,x2,p2]-wig[-x1,-p1,-x2,-p2]]

Parity symmetric:

Simplify[wig[x1,p1,x2,p2]-wig[-x1,-p1,x2,p2]]

PT transformation:

Wigner2mod[QuantumPartialTranspose[ψ,{2}]["DensityMatrix"],α,Conjugate[β]]-Wigner2mod[ψ["DensityMatrix"],α,β]//Simplify

Separability:

Wigner2mod[ψ["DensityMatrix"],α,β]-wignerTrace[FockState[2]["DensityMatrix"],α]wignerTrace[FockState[3]["DensityMatrix"],β]

NOON State

SetFockSpaceSize[10];

NOON state with N = 4:

ψ =(1/(Sqrt[2]))( FockState[{0,4}]+FockState[{4,0}]);

Joint Wigner function:

Wfunc = Function[{x1,p1,x2,p2},Simplify[Wigner2mod[ψ["DensityMatrix"],((x1+I p1)/(Sqrt[2])),((x2+I p2)/(Sqrt[2]))],Element[{x1,p1,x2,p2},Reals]]];

Check global symmetry:

Wfunc[x1,p1,x2,p2]-Wfunc[-x1,-p1,-x2,-p2]

Parity symmetry for NOON states only work when N is even:

Simplify[Wfunc[x1,p1,x2,p2]-Wfunc[-x1,-p1,x2,p2]]

3D slices ( (Subscript[p, 2])=0):

func = Wigner2mod[ψ["DensityMatrix"],α,β]/.{α->((x1+I p1)/(Sqrt[2])),β->(x2/(Sqrt[2]))};

DensityPlot3D[Re@func,{x1,-3,3},{p1,-3,3},{x2,-3,3},]

GIFS:

SystemOpen["C:\\Users\\Lenovo\\Documents\\NOON_WignerContour.GIF"]

SystemOpen["C:\\Users\\Lenovo\\Documents\\NOON_Wigner400.GIF"]

Joint Wigner representation of bright states

Pick bright state with 𝒩 = n = 4:

ρ = ψNn[4,4][0]["DensityMatrix"];

Wigner function in terms of symbolic variables:

func = Wigner2mod[ρ,α,β]/.{α->((x1+I p1)/(Sqrt[2])),β->(x2/(Sqrt[2]))};

Visualize the 3D slices:

DensityPlot3D[Re@func,{x1,-3,3},{p1,-3,3},{x2,-3,3},]

(Legended[, (BarLegend)])

X slices visualization:

SliceDensityPlot3D[Re@func,{"XStackedPlanes",5},{x1,-3,3},{p1,-3,3},{x2,-3,3},]

References

Villas-Boas, C. J., Máximo, C. E., Paulino, P. J., Bachelard, R. P., & Rempe, G. (2025). Bright and Dark States of Light: The Quantum Origin of Classical Interference. Physical Review Letters, 134(13), 133603. https://doi.org/10.1103/PhysRevLett.134.133603

Bright and dark states of light: reinterpreting interference https://community.wolfram.com/groups/-/m/t/3461741

Exploring non-commutative algebra: commutators, BCH expansion and calculations of bosonic algebra https://community.wolfram.com/groups/-/m/t/3533713

Two-mode phase space representation: joint Wigner function https://community.wolfram.com/groups/-/m/t/3542358

Thanks !

brunot@wolfram.com

$Assumptions = θ ∈ Reals;

wignerTrace[ρ_SparseArray, α_]:= With[{pos = ρ["ExplicitPositions"], vals= ρ["ExplicitValues"]}, (2/π)Dot[vals, 𝒦mn[α, {#2-1,#1-1}]& @@@ pos]]