Dynamics and QuantumOptics.jl - Stefan Krastanov

This package should be thought of as a domain specific linear algebra wrapper.

Operator basics

using QuantumOpticsBase

cutoff = 25

F = FockBasis(cutoff)

Ψ = fockstate(F,5)

Ket(dim=26)
  basis: Fock(cutoff=25)
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
     ⋮
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

S = SpinBasis(1//2)
s = spindown(S)

Ket(dim=2)
  basis: Spin(1/2)
 0.0 + 0.0im
 1.0 + 0.0im

G = GenericBasis(100)
basisstate(G,2)

Ket(dim=100)
  basis: Basis(dim=100)
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
     ⋮
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

C = F ⊗ S

sparse(dm(Ψ⊗s))

Operator(dim=52x52)
  basis: [Fock(cutoff=25) ⊗ Spin(1/2)]
⎡⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎦

projector(Ψ⊗s)

Operator(dim=52x52)
  basis: [Fock(cutoff=25) ⊗ Spin(1/2)]
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
    ⋮                             ⋱                ⋮       
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im

using BenchmarkTools
P = projector(Ψ)
Ps = sparse(P)

@benchmark P*Ψ

BenchmarkTools.Trial: 10000 samples with 310 evaluations per sample.
 Range (min … max):  270.926 ns …  23.414 μs  ┊ GC (min … max): 0.00% … 97.55%
 Time  (median):     275.613 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   302.371 ns ± 298.936 ns  ┊ GC (mean ± σ):  6.10% ±  8.88%

  █▁▄▂        ▁                                                 ▁
  ████▆▆▄▄▃▃▅▃██▆▆▃▃▃▃▁▁▁▁▁▃▁▁▁▁▁▁▃▁▁▄▃▃▃▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▃▅▅▄▄▄▆ █
  271 ns        Histogram: log(frequency) by time       1.22 μs <

 Memory estimate: 528 bytes, allocs estimate: 3.

@benchmark Ps*Ψ

BenchmarkTools.Trial: 10000 samples with 967 evaluations per sample.
 Range (min … max):   80.960 ns …  10.027 μs  ┊ GC (min … max):  0.00% … 97.32%
 Time  (median):      84.429 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   109.826 ns ± 175.276 ns  ┊ GC (mean ± σ):  15.45% ± 13.19%

  █▅               ▂▁                                           ▁
  ███▆▃▄▄▄▄▄▁▁▁▁▄▃███▆▄▁▃▁▁▁▃▁▁▃▁▁▁▁▁▁▁▃▅▄▅▄▄▅▁▄▄▆▄▆▆▅▅▆▄▄▄▅▆▅▅ █
  81 ns         Histogram: log(frequency) by time        733 ns <

 Memory estimate: 528 bytes, allocs estimate: 3.

We can use in-place or mutating functions to improve performance. This is one of the first things to think about optimizing when looking for speedups.

using LinearAlgebra: mul!

buffer = copy(Ψ)
@benchmark mul!(buffer, Ps, Ψ)

BenchmarkTools.Trial: 10000 samples with 990 evaluations per sample.
 Range (min … max):  43.921 ns … 124.113 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     44.639 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   45.082 ns ±   2.256 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▂▅▇██▇▅▃▁▁                                          ▁▁▁▂    ▂
  ▇███████████▇▄▁▅▅▇▆▇▄▁▃▁▃▁▁▁▁▃▁▁▃▁▁▃▁▃▁▁▁▁▃▁▁▁▁▄▇▁▁▅▇██████▆ █
  43.9 ns       Histogram: log(frequency) by time      52.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Composite spaces

We can create composite spaces and operators with the '\otimes' symbol representing the tensor product.

psi = fockstate(F,4)⊗spinup(S)

Ket(dim=52)
  basis: [Fock(cutoff=25) ⊗ Spin(1/2)]
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
     ⋮
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

create(F)⊗identityoperator(S)

Operator(dim=52x52)
  basis: [Fock(cutoff=25) ⊗ Spin(1/2)]
⎡⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⎦

This does not scale well! We have a function called embed to address this problem

embed(F⊗S, 1, create(F))

Operator(dim=52x52)
  basis: [Fock(cutoff=25) ⊗ Spin(1/2)]
⎡⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⎦

We can take a partial trace with ptrace()

l = spindown(S)
h = spinup(S)

bell = (l⊗l + h⊗h) / √2

ptrace(bell,2)

Operator(dim=2x2)
  basis: Spin(1/2)
 0.5+0.0im  0.0+0.0im
 0.0+0.0im  0.5+0.0im

Coherent states

using QuantumOptics

QuantumOptics depends on DifferentialEquations, so has a much longer precompile time compared to QuantumOpticsBase, which defines the methods for working with quantum states and operators without any time evolution.

We will work with the quantum harmonic oscillator and look at coherent states in phase space.

cutoff = 32
B = FockBasis(cutoff)
ω = 1.0
a = destroy(B)

H = ω*a'*a

Operator(dim=33x33)
  basis: Fock(cutoff=32)
⎡⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⎦

Now we can make a coherent state

α = 4
ψ = coherentstate(B, α)

Ket(dim=33)
  basis: Fock(cutoff=32)
 0.00033546262790251 + 0.0im
 0.00134185051161005 + 0.0im
 0.00379532638439241 + 0.0im
 0.00876493083876579 + 0.0im
 0.01752986167753158 + 0.0im
 0.03135836987770328 + 0.0im
 0.05120800357722345 + 0.0im
 0.07741922434367954 + 0.0im
 0.10948731705523689 + 0.0im
 0.14598308940698249 + 0.0im
                     ⋮
  0.1198756849197267 + 0.0im
 0.09590054793578136 + 0.0im
 0.07523057927519837 + 0.0im
 0.05791252692776968 + 0.0im
 0.04377775544177429 + 0.0im
  0.0325173004055702 + 0.0im
 0.02374727858840358 + 0.0im
 0.01706054857273676 + 0.0im
 0.01206362958654464 + 0.0im

a*ψ ≈ α*ψ

false

TODO: what is going wrong here?

using CairoMakie
quad = -10:0.1:10

w = wigner(ψ,quad,quad)

fig = Figure()
ax = Axis(fig[1,1], aspect = DataAspect())
heatmap!(quad,quad,w)
fig

Schrodinger dynamics

ts = 0:0.1:3*2*π

_, ψs = timeevolution.schroedinger(ts,ψ,H)

x = (a' + a)/√2
p = 1im*(a' - a)/√2

fig = Figure()
ax = Axis(fig[1,1], aspect = DataAspect())
lines!(ax,ts,real.(expect.((x,),ψs)))
lines!(ax,ts,real.(expect.((p,),ψs)))
fig

lines(norm.(ψs))

First, the naive way to create a time-dependent Hamiltonian

function Hdynamic(t,psi)
    return H+p*4*sin(8*ω*t)
end

_, ψs = timeevolution.schroedinger_dynamic(ts,ψ,Hdynamic)

fig = Figure()
ax = Axis(fig[1,1], aspect = DataAspect())
lines!(ax,ts,real.(expect.((x,),ψs)))
lines!(ax,ts,real.(expect.((p,),ψs)))
fig

SciMLOperators has an interesting way to create time-dependent operators in an allocation-efficient way.

using BenchmarkTools
Hlazy = H + TimeDependentSum(t-> 4*sin(8*ω*t),p)
@benchmark _, ψs = timeevolution.schroedinger_dynamic(ts, ψ, Hlazy)

BenchmarkTools.Trial: 241 samples with 1 evaluation per sample.
 Range (min … max):  19.485 ms … 49.194 ms  ┊ GC (min … max): 0.00% … 55.77%
 Time  (median):     20.418 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   20.775 ms ±  2.812 ms  ┊ GC (mean ± σ):  1.70% ±  6.14%

  █▁ ▃ ▆▄                                                      
  ██▆████▇▇█▆▅▄▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▅
  19.5 ms      Histogram: log(frequency) by time      35.5 ms <

 Memory estimate: 3.94 MiB, allocs estimate: 161486.

And the old way:

@benchmark _, ψs = timeevolution.schroedinger_dynamic(ts,ψ,Hdynamic)

BenchmarkTools.Trial: 110 samples with 1 evaluation per sample.
 Range (min … max):  40.350 ms … 106.050 ms  ┊ GC (min … max): 10.19% … 20.43%
 Time  (median):     41.595 ms               ┊ GC (median):    11.40%
 Time  (mean ± σ):   45.719 ms ±  13.218 ms  ┊ GC (mean ± σ):  12.99% ±  4.37%

  █▇▃                                                           
  ███▄▁▄▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▄▄▁▁▁▄▁▁▁▁▁▁▁▄▁▁▁▁▁▄▄▁▁▄▁▄▁▁▁▁▁▁▁▁▁▁▁▄ ▄
  40.3 ms       Histogram: log(frequency) by time       106 ms <

 Memory estimate: 121.67 MiB, allocs estimate: 439236.

Trying out different solvers

using OrdinaryDiffEqLowOrderRK: DP5
using OrdinaryDiffEqTsit5: Tsit5
using OrdinaryDiffEqVerner: Vern8

A good discussion on Schrodinger dynamics and solver tolerance

@benchmark _, ψs = timeevolution.schroedinger_dynamic(ts, ψ, Hlazy; alg =
DP5())

BenchmarkTools.Trial: 237 samples with 1 evaluation per sample.
 Range (min … max):  19.589 ms … 53.518 ms  ┊ GC (min … max): 0.00% … 60.07%
 Time  (median):     21.073 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.143 ms ±  3.330 ms  ┊ GC (mean ± σ):  2.00% ±  6.90%

  █   ▆                                                        
  █▃▄▄█▅▃▅▃▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂ ▂
  19.6 ms         Histogram: frequency by time        40.4 ms <

 Memory estimate: 3.94 MiB, allocs estimate: 161486.

@benchmark _, ψs = timeevolution.schroedinger_dynamic(ts, ψ, Hlazy; alg =
Tsit5())

BenchmarkTools.Trial: 228 samples with 1 evaluation per sample.
 Range (min … max):  20.604 ms … 51.843 ms  ┊ GC (min … max): 0.00% … 56.61%
 Time  (median):     21.303 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.995 ms ±  3.216 ms  ┊ GC (mean ± σ):  1.87% ±  6.63%

  █▂   █▁                                                      
  ██▆▄▆██▄▄▄▅▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▅
  20.6 ms      Histogram: log(frequency) by time      38.9 ms <

 Memory estimate: 4.19 MiB, allocs estimate: 172048.

@benchmark _, ψs = timeevolution.schroedinger_dynamic(ts, ψ, Hlazy; alg =
Vern8())

BenchmarkTools.Trial: 293 samples with 1 evaluation per sample.
 Range (min … max):  15.716 ms … 44.843 ms  ┊ GC (min … max): 0.00% … 61.06%
 Time  (median):     16.465 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   17.093 ms ±  2.582 ms  ┊ GC (mean ± σ):  2.18% ±  6.81%

  █  ▁                                                         
  █▇▆██▄▃▃▃▆▅▃▅▃▃▃▃▂▁▂▁▂▂▁▁▁▁▁▂▁▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▂ ▂
  15.7 ms         Histogram: frequency by time        26.6 ms <

 Memory estimate: 4.12 MiB, allocs estimate: 130192.

Master equation evolution

γ = 0.1
_, ρs = timeevolution.master(ts,ψ,H,[√γ*a])

fig = Figure()
ax = Axis(fig[1,1], aspect = DataAspect())
lines!(ax,ts,real.(expect.((x,),ρs)))
lines!(ax,ts,real.(expect.((p,),ρs)))
fig

TODO: Lindblad is nonlinear with respect to collapse operators?