Negative photon number with a coherent driving

[eric940109]

Dear QuantumCumulants developer,

I attempted to add a coherent driving term \eta*(a'-a) in the Hamiltonian, and the resultant time evolution of the photon number <a'a> contains huge negative values when \eta was simply chosen to be 1 (ode solver failed when \eta becomes larger, e.g. 10^5). When \eta is 0, the time evolution of the photon number looks correct. Please see the attached figure.

May I ask whether QuantumCumulants is suitable for simulating the coupling between many atoms and a single cavity mode under a coherent driving of the cavity mode?

Thank you!

Best wishes,
Eric

[ChristophHotter]

Hi @eric940109

Yes QuantumCumulants is usually well suited for this problem (depending on the parameter regime).
I think in your case the mistake is that the Hamiltonian is not hermitian (which is unphysical).
For a real eta (eta = 1) the drive term needs to be eta (a' + a) or i eta (a' - a).

[eric940109]

Hi @ChristophHotter

Thanks a lot for your comment! After correcting the drive term, the negative photon numbers seem to still remain. Does this mean my setting of parameters are out of the reasonable parameter regime?

The code and typical results are attached.

M = 2 #order
@cnumbers N Δc Δs g κ γ12 η γ22
hf = FockSpace(:cavity)
ha1 = NLevelSpace(:atom, 2)
ha = ClusterSpace(ha1, N, M)
h = tensor(hf, ha);

@qnumbers a::Destroy(h)
σ(i, j, k) = Transition(h, :σ, i, j)[k]

#Hamiltonian
H = Δc*a'*a + Δs*sum(σ(2,1,i)*σ(1,2,i) for i = 1:M) + g*sum(a'*σ(1,2,i) + a*σ(2,1,i) for i=1:M) + η*(a'+a)

#Collapse operators
J = [a;[σ(1,2,i) for i = 1:M];[σ(2,2,i) for i = 1:M]]
rates = [κ;[γ12 for i = 1:M];[γ22 for i = 1:M]]

#Derive equation
eqs = meanfield([σ(2,2,1)],H,J;rates = rates, order = 2)
eqs_c = complete(eqs, order = 2)

using OrdinaryDiffEq, ModelingToolkit
@named sys = ODESystem(eqs_c)
u0 = zeros(ComplexF64, length(eqs_c))
u0[1] = 1.0
ps = [N, Δc, Δs, g, κ, γ12, η, γ22] #N Δc Δs g κ γ12 η γ22
p0 = [1e14, 0.0, 0.0, 2*π*0.022, 2*π*0.18e6, 1e4, 1e5, 2*π*0.42e6]
prob0 = ODEProblem(sys, u0, (0.0, 8e-5), ps.=>p0)
alg = RK4()
sol0 = solve(prob0, alg)
t0 = sol0.t
q = real.(sol0[a'*a])
using Plots
plot(t0,q)```

Best wishes,
Eric

[ChristophHotter]

The program itself seems fine.
I think the problem is your parameter set, you have very small and very large numbers at the same time, which is usually hard to solve numerically.

[ChristophHotter]

I just saw this again. You could also try to reduce the tolerance of the numerics with e.g. sol0 = solve(prob0, Tsit5(); abstol=1e-10, reltol=1e-10.

[eric940109]

@ChristophHotter Thanks a lot for your reply! I will have it a go!