Spurious results with steadystate.iterative for non-linear system

[lukasgrunwald]

Considering the steady state of non-linear systems, e.g. a quartic Harmonic oscillator

$$H = \frac{p^2}{ 2} + \Omega^2 \frac{x^2}{ 2} + k x^4 \sim \Omega n + \tilde{k} (a + a^\dagger)^4$$

where we are interested in $\langle n(t) \rangle$, methods steadystate.master and steadystate.eigenvector produce accurate and coinciding results for the steady state, while the iterative solver fails to find the appropriate nullspace. The result of the iterative solver depends strongly on the random seed. I'm not sure if I'm missing something or if this is a bug or something else, but everything seems to work fine for linear systems, viz. without the $x^4$ term

Please find a code snippet below for reproducing these results. (The last section checks if a given density matrix indeed has (approximately) zero eigenvalue

using QuantumOptics
using Random; Random.seed!(2)

basis = FockBasis(10)
f = 0.5
Ω = 0.5
nl = 0.5

n = number(basis)
b = destroy(basis)
bt = create(basis)
H(x::Integer) = dense(Ω * n + f * (b + bt) + nl * (b + bt)^x)
J = [b, ]

# Works for linear systems
H_linear = H(2)
tout, ρ_master = steadystate.master(H_linear, J)
ρ_eig = steadystate.eigenvector(H_linear, J)
ρ_it = steadystate.iterative(H_linear,J)

println("--Linear Systems")
println("n_master = ", real(expect(n, ρ_master[end])))
println("n_eig = ", real(expect(n, ρ_eig)))
println("n_it = ", real(expect(n, ρ_it)))

# Doesn't work for non-linear systems
H_non_linear = H(4)
tout, ρ_master = steadystate.master(H_non_linear, J)
ρ_eig = steadystate.eigenvector(H_non_linear, J)
ρ_it = steadystate.iterative(H_non_linear,J)

println("\n--Non-Linear Systems")
println("n_master = ", real(expect(n, ρ_master[end])))
println("n_eig = ", real(expect(n, ρ_eig)))
println("n_it = ", real(expect(n, ρ_it)))

# Nullspace structure
L = liouvillian(H_non_linear, J).data
M = size(ρ_eig, 1)

vec_master = reshape(ρ_master[end].data, M^2)
vec_eig = reshape(ρ_eig.data, M^2)
vec_it = reshape(ρ_it.data, M^2)

println("\n--Null-space deviation (Should be Δmax ≈ 0)")
println("Δmax_master = ", maximum(abs.(L * vec_master)))
println("Δmax_eigen = ", maximum(abs.(L * vec_eig)))
println("Δmax_iter = ", maximum(abs.(L * vec_it)))

VersionInfo

Julia Version 1.10.5
Commit 6f3fdf7b362 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M3
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 8 default, 0 interactive, 4 GC (on 4 virtual cores)
Environment:
  JULIA_NUM_THREADS = 8
  JULIA_EDITOR = code

and [6e0679c1] QuantumOptics v1.2.1

[david-pl]

Have you tried setting the initial state? You can just pass thenrho0 kwarg to the iterative solver. It may that the default that is used leads to different states depending on the "path" the iterative solver takes (which is chosen using some randomness).

[lukasgrunwald]

Thanks for the quick reply! Yes, I also tried using an initial state (in this case the ground state and 1st excited state), but it didn't help. The state found with the iterative method does not lie in the nullspace of the liouvillian, but it should for any given initial state right?
In this example, the steady state is also unique, viz. there is only one eigenvalue $\lambda$ with $\Re(\lambda) = 0$. So I would expect the iterative solver to at least get close to that state. I'm a bit lost on why it doesn't.