How to plot Kerr resonator steady state photons

[albertomercurio]

Hello,

I have the following case of a single-photon driven Kerr resonator.

using QuantumCumulants
using HarmonicBalance
using CairoMakie

# %%

L = 10
κ = 0.1
U = 10κ / L
Δ = -10κ

h = FockSpace(:cavity)
@qnumbers a::Destroy(h)
@variables Δ₂::Real U₂::Real F₂::Real κ₂::Real
param = [Δ₂, U₂, F₂, κ₂]

H_qc = Δ₂* a' * a + U₂ * (a'^2 * a^2) / 2 + F₂ * (a + a')
ops_qc = [a, a']

# %%

eqs = meanfield(ops_qc, H_qc, [a]; rates=[κ₂], order=1)
eqs_completed = complete(eqs)

# %%

F_list = range(0, 10κ, 30)

fixed = (U₂ => U, κ₂ => κ, Δ₂ => Δ)
varied = (F₂ => F_list .* √(L) )
problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(
    eqs_completed, param, varied, fixed
)

# %%

result = get_steady_states(problem_c1, TotalDegree())
y = get_solutions(result, string(result.problem.variables[2] * result.problem.variables[1]); realify=false)

I don't know if I extracted the steady state photon number correctly, I see there are two variables and I guess they are $a^\dagger$ and $a$.

The system should show a dissipative phase transition around F = 2, with zero photons below a threshold and a large number of photons after this point. I'm not expecting this to be fully captured by the mean field equations, but at least some meaningful features. Here instead I just see some quadratic function.

Am I doing something wrong?

[oameye]

Hey Alberto,

Plotting the results with:

using QuantumCumulants
using HarmonicBalance
using CairoMakie

# %%

L = 10
κ = 0.1
U = 10κ / L
Δ = -10κ

h = FockSpace(:cavity)
@qnumbers a::Destroy(h)
@variables Δ₂::Real U₂::Real F₂::Real κ₂::Real
param = [Δ₂, U₂, F₂, κ₂]

H_qc = Δ₂* a' * a + U₂ * (a'^2 * a^2) / 2 + F₂ * (a + a')
ops_qc = [a, a']

# %%

eqs = meanfield(ops_qc, H_qc, [a]; rates=[κ₂], order=1)
eqs_completed = complete(eqs)

# %%

F_list = range(0, 10κ, 1000)

fixed = (U₂ => U, κ₂ => κ, Δ₂ => Δ)
varied = (F₂ => F_list .* √(L) )
problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(
    eqs_completed, param, varied, fixed
)

# %%

result = get_steady_states(problem_c1, TotalDegree())

branches = get_branches( result, "(aᵣ^2+aᵢ^2)/2"; realify=true, class=["stable"], not_class=[])
branches_unstable = get_branches( result, "(aᵣ^2+aᵢ^2)/2"; realify=true, class=["physical"], not_class=["stable"])

# %%

fig = Figure()
ax = Axis(fig[1, 1], xlabel="F", ylabel="n")

for i in eachindex(branches)
    all(isnan, branches[i]) && continue
    lines!(ax, F_list .* √(L), branches[i], linewidth=2, label="branch $i", color=Cycled(i))
end
for i in eachindex(branches_unstable)
    all(isnan, branches_unstable[i]) && continue
    lines!(ax, F_list .* √(L), branches_unstable[i], linewidth=2, label="branch $i", color=Cycled(i), linestyle=:dash)
end
Legend(fig[1,2], ax)
fig

I get

Which is the usual hysteresis of the Kerr resonator. You are far too much detuned to see the high-amplitude branch to become unstable. Indeed, if we plot the number of stable states in a $\Delta-F$ sweep, we get:

fixed = (U₂ => U, κ₂ => κ)
varied = (F₂ =>range(0, 10κ, 100) .* √(L),  Δ₂ => range(-15κ, 5κ, 100) .* √(L))
problem = HarmonicSteadyState.HomotopyContinuationProblem(
    eqs_completed, param, varied, fixed
)
result = get_steady_states(problem, WarmUp())

phase_m = phase_diagram(result, class=["stable"])

fig = Figure()
ax = Axis(fig[1, 1], xlabel="Δ", ylabel="F")
heatmap!(ax, range(-15κ, 5κ, 100) .* √(L), range(0, 10κ, 100) .* √(L), phase_m')
fig

I think you are referring too "quantum" dissipative phase transition described by Ciuti. I have looked at this many times before. This fluctuation induced phase transition will have no remnant indications in the mean-field limit. This is clear when plot the phase diagram next two the photon number computed with you package :). I copy here some plots from an old notebook:

[albertomercurio]

Thanks a lot! Indeed I'm aware that mean field doesn't reproduce it properly, but just need that curve.

What is the definition of $a_{r,i}$ in the code? Are they the real and imaginary parts of the coherence? If yes, shouldn't the number of photons be expressed as $a_r^2 + a_i^2$ instead (without the 1/2 factor)?

Now, I have the last question. Suppose that I want to do the same calculation but not in the drive frame, the Hamiltonian is time-dependent

$$ H(t) = \omega_0 a^\dagger a + \frac{U}{2} a^{\dagger 2} a^2 + F (a e^{i \omega_d t} + a^\dagger e^{-i \omega_d t}) $$

How can I simulate this? It should give the same results as it is linked by a unitary transformation.

I know that QuantumCumulants.jl doesn't support well time dependent terms. Should I directly define the equations in HarmonicBalance.jl?

[albertomercurio]

I tried to put $t$ as an independent variable, as suggested here

• timedependent parameters qojulia/QuantumCumulants.jl#93

So, like this

using Symbolics
using QuantumCumulants
using HarmonicBalance
using CairoMakie

# %%

L = 10
κ = 0.1
U = 10κ / L
Δ = -10κ

h = FockSpace(:cavity)
@syms t::Real
@qnumbers a::Destroy(h)
@variables ω0₂::Real ωd₂::Real Δ₂::Real U₂::Real F₂::Real κ₂::Real
param = [ω0₂, ωd₂, Δ₂, U₂, F₂, κ₂]

# H_qc = Δ₂* a' * a + U₂ * (a'^2 * a^2) / 2 + F₂ * (a + a')
H_qc = ω0₂ * a' * a + U₂ * (a'^2 * a^2) / 2 + F₂ * (a * exp(1im * ωd₂ * t) + a' * exp(-1im * ωd₂ * t))
ops_qc = [a, a']

# %%

eqs = meanfield(ops_qc, H_qc, [a]; rates=[κ₂], iv=t, order=1)
eqs_completed = complete(eqs)

# %%

F_list = range(0, 10κ, 1000)

fixed = (U₂ => U, κ₂ => κ, Δ₂ => Δ)
varied = (F₂ => F_list .* √(L) )
problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(
    eqs_completed, param, varied, fixed
)

But I get the error

ERROR: UndefVarError: `t` not defined in `HarmonicSteadyState.HC_wrapper`
Suggestion: check for spelling errors or missing imports.

Perhaps it is better to write the equations manually?

[oameye]

Ha yeah, the definition of $a_{r,i}$ is not that well documented. It is kind-off a niche feature. Every bosonic operator is expanded as $a= (a_r + i a_r)/\sqrt{2}$. I don't know why I normalized with $\sqrt{2}$ and probably should done with. I will document it now.

At the moment, we don't have an interface for time-dependent Hamiltonian system. Harmonicbalance.jl only works with algebraic equations of motion. I started working in this in FloquetExpansions.jl using SecondQuantizedAlgebra.jl. But the code of QuantumCumulant.jl was a bit spaghetti due to the indexing feature. Once I get QuantumCumulant.jl through for MTKv10. I hope to clean up SecondQuantizedAlgebra.jl and get FloquetExpansions.jl working for higher orders Floquet.

I would be cool to have a MTK for Quantum optic systems. I envision doing a perturbation theory on some second quantized liouvillian system, port it to HarmonicSteadyState.jl for some meanfield analyses and use QuantumToolbox.jl for a though Quantum analyses.

[oameye]

I tried to put t as an independent variable, as suggested here

• timedependent parameters qojulia/QuantumCumulants.jl#93

So, like this

using Symbolics
using QuantumCumulants
using HarmonicBalance
using CairoMakie

%%

L = 10
κ = 0.1
U = 10κ / L
Δ = -10κ

h = FockSpace(:cavity)
@syms t::Real
@qnumbers a::Destroy(h)
@variables ω0₂::Real ωd₂::Real Δ₂::Real U₂::Real F₂::Real κ₂::Real
param = [ω0₂, ωd₂, Δ₂, U₂, F₂, κ₂]

H_qc = Δ₂* a' * a + U₂ * (a'^2 * a^2) / 2 + F₂ * (a + a')

H_qc = ω0₂ * a' * a + U₂ * (a'^2 * a^2) / 2 + F₂ * (a * exp(1im * ωd₂ * t) + a' * exp(-1im * ωd₂ * t))
ops_qc = [a, a']

%%

eqs = meanfield(ops_qc, H_qc, [a]; rates=[κ₂], iv=t, order=1)
eqs_completed = complete(eqs)

%%

F_list = range(0, 10κ, 1000)

fixed = (U₂ => U, κ₂ => κ, Δ₂ => Δ)
varied = (F₂ => F_list .* √(L) )
problem_c1 = HarmonicSteadyState.HomotopyContinuationProblem(
eqs_completed, param, varied, fixed
)
But I get the error

ERROR: UndefVarError: t not defined in HarmonicSteadyState.HC_wrapper
Suggestion: check for spelling errors or missing imports.
Perhaps it is better to write the equations manually?

HarmonicSteadyState already assumes a time independent system. The goal is to find all steady-state and properties thereoff knowing the system is in a rotating frame.

[oameye]

So in short at the moment, we don't really support what you are looking for and it is something I would like to achieve with FloquetExpansions.jl. However, It is not their yet. And so it is easier at the moment to do you RWA yourself and just enter the time-dependent Hamiltonian in HarmonicSteadyState.jl.

[oameye]

Also if you think QuantumToolbox.jl could also benefit from SecondQuantizedAlgrabra.jl. I am happy to collab/discuss. Anyway I realised that my Phd is also not that long anymore so I will have to make priorities :p

[albertomercurio]

Thanks for the reply.

Also if you think QuantumToolbox.jl could also benefit from SecondQuantizedAlgrabra.jl. I am happy to collab/discuss.

That would be nice. I have a few cases where we could use symbolic expressions to be used in QuantumToolbox.jl, like the simulation of open systems using the Truncated Wigner Approximation. I would say that we could discuss once you are more free from the PhD then 😆.

So in short at the moment, we don't really support what you are looking for and it is something I would like to achieve with FloquetExpansions.jl. However, It is not their yet. And so it is easier at the moment to do you RWA yourself and just enter the time-dependent Hamiltonian in HarmonicSteadyState.jl.

My time-dependent example was just a simple case to compare the time-independent and time-dependent cases, as they are linked by a unitary transformation. However, the system I would like o study cannot remove time dependence due to the presence of counter rotating terms in the Hamiltonian.

I was wondering if I can simply follow the basic example here

• https://quantumengineeredsystems.github.io/HarmonicBalance.jl/v0.15/introduction/

where you simulate a time-dependent differential equation. My case should be very similar, where I should maybe split it into two real differential equations of $\alpha_r$ and $\alpha_i$. Isn't it possible? I mean it is just a system of differential equations.

[oameye]

I was wondering if I can simply follow the basic example here

• quantumengineeredsystems.github.io/HarmonicBalance.jl/v0.15/introduction

where you simulate a time-dependent differential equation. My case should be very similar, where I should maybe split it into two real differential equations of α r and α i . Isn't it possible? I mean it is just a system of differential equations.

Well we don't simulate the time-dependent differential equation. For that you can just use QC.jl with OrdinaryDiffEq.jl. We apply a symbolic perturbation theory, which makes the system time-independent and then solve for the steady-states using Homotopycontinuation.jl.

I would like to study cannot remove time dependence due to the presence of counter rotating terms in the Hamiltonian.

If you can map you system to a second differential equation of harmonic oscillator + correction. You can use the interface in HB.jl as you saw in the docs to do RWA. At the moment we don't have a proper interface for first order diff-eq to do symbolic expansion. Although this would actually not be that hard, but nobody has it implemented yet

So at this point it is just easier to do RWA yourself and use QC.jl symbolics to define the system in HarmonicSteadyState.jl. If you want to just directly want to give the realified equations you can use this non-public function.

If you don't want to get rid of the counter-rotating terms this is not what the package does.
Instead, you use PeriodicOrbits.jl from the DynamicalSystem ecosystem.

[albertomercurio]

Ok, I see. Thank you for your time.

[oameye]

I have changed the definition in $a_r$ and $a_i$ in #80