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.761 ns …  27.236 μs  ┊ GC (min … max): 0.00% … 97.88%
 Time  (median):     274.900 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   306.462 ns ± 355.662 ns  ┊ GC (mean ± σ):  6.43% ±  8.58%

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

 Memory estimate: 528 bytes, allocs estimate: 3.

@benchmark Ps*Ψ

BenchmarkTools.Trial: 10000 samples with 966 evaluations per sample.
 Range (min … max):   80.233 ns …  10.894 μs  ┊ GC (min … max):  0.00% … 98.28%
 Time  (median):      84.122 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   116.372 ns ± 241.468 ns  ┊ GC (mean ± σ):  16.26% ± 12.52%

  █▄▁            ▂▃▁                                            ▁
  ███▆▅▃▃▁▁▁▁▃▃▃▃███▆▄▃▃▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▄▄▄▄▄▄▄▄▄▄▄▅▄▄▄▄▄▄▅▄▆▆▆ █
  80.2 ns       Histogram: log(frequency) by time        808 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 989 evaluations per sample.
 Range (min … max):  46.011 ns … 225.010 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     48.990 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   49.498 ns ±   2.800 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▂▁▁▃▄▆▇██▇▇▅▄▃▃▂▁▁  ▂                       ▁▂▂▂▁▁▁   ▃
  ▄▅▁▅▆▆▆█████████████████████▅▅▆▃▃▁▁▁▁▁▃▁▁▁▁▁▄▁▃▁▄▅▇███████▇▇ █
  46 ns         Histogram: log(frequency) by time      57.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: 238 samples with 1 evaluation per sample.
 Range (min … max):  18.902 ms … 50.949 ms  ┊ GC (min … max): 0.00% … 59.58%
 Time  (median):     20.622 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.007 ms ±  2.957 ms  ┊ GC (mean ± σ):  1.79% ±  6.44%

  ▄█  ▅ █  ▄▂                                                  
  ██▄████▅███▄▂▃▆▃▁▃▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂ ▃
  18.9 ms         Histogram: 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: 112 samples with 1 evaluation per sample.
 Range (min … max):  38.850 ms … 110.757 ms  ┊ GC (min … max): 10.62% … 24.80%
 Time  (median):     40.539 ms               ┊ GC (median):    11.34%
 Time  (mean ± σ):   44.687 ms ±  13.574 ms  ┊ GC (mean ± σ):  13.26% ±  4.40%

  ██▇▃                                                          
  ████▁▁▁▁▁▅▁▁▁▁▅▁▁▁▁▁▁▁▁▁▁▁▁▅▁▁▁▁▁▁▁▁▁▅▁▁▁▅▁▁▁▅▁▅▅▁▅▁▁▅▁▁▁▁▁▅ ▅
  38.8 ms       Histogram: log(frequency) by time      97.1 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: 235 samples with 1 evaluation per sample.
 Range (min … max):  19.082 ms … 57.405 ms  ┊ GC (min … max): 0.00% … 60.90%
 Time  (median):     20.746 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.320 ms ±  3.562 ms  ┊ GC (mean ± σ):  1.95% ±  6.85%

  ▃  ▁█   ▁                                                    
  █▃▄██▆▅██▃▄▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂ ▃
  19.1 ms         Histogram: frequency by time        39.3 ms <

 Memory estimate: 3.94 MiB, allocs estimate: 161486.

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

BenchmarkTools.Trial: 231 samples with 1 evaluation per sample.
 Range (min … max):  19.965 ms … 51.499 ms  ┊ GC (min … max): 0.00% … 60.96%
 Time  (median):     21.663 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.695 ms ±  3.451 ms  ┊ GC (mean ± σ):  2.02% ±  7.08%

  █▆   █▅▁                                                     
  ██▇▆▇████▄▄█▄▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▄▄▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▁▁▁▁▁▁▁▁▁▁▁▄ ▆
  20 ms        Histogram: log(frequency) by time      38.4 ms <

 Memory estimate: 4.19 MiB, allocs estimate: 172048.

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

BenchmarkTools.Trial: 298 samples with 1 evaluation per sample.
 Range (min … max):  15.437 ms … 47.735 ms  ┊ GC (min … max): 0.00% … 64.67%
 Time  (median):     16.004 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   16.820 ms ±  2.967 ms  ┊ GC (mean ± σ):  2.33% ±  7.20%

  █▇█▂ ▁ ▅▇▅                                                   
  ████▄█▄████▄▄▄▄▁▁▄▄▄▁▁▁▄▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▆
  15.4 ms      Histogram: log(frequency) by time      29.2 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?