Solving with matrix exponential

[mabuni1998]

@Krastanov
As we have talked about, we can also solve the differential equation system by continuously applying the unitaryU_n = exp(-i delta t H) to every time-bin. Since our method is only precise up until dt, this should work fairly well and the potential performance gain is enormous. We will be probably sacrificing numerical stability and precision though.

To get a crasp of the performance gain I made the following example where we solve the evolution of a single-photon scattering on a cavity. Due to the low number of excitations, it is sufficient to expand U_n up to first order in t_n. Whether this will always hold is not clear to me.

using QuantumOptics
using WaveguideQED

δ = 0
γ = 1
times = 0:0.1:10
dt = times[2] - times[1]

#Create operators for two photons interacting with cavity
bc = FockBasis(1)
bw = WaveguideBasis(1,times)
a = destroy(bc)
ad = create(bc);
n = ad*a ⊗ identityoperator(bw)
#$w†a and a†w efficient implementation$
wda = emission(bc,bw)
adw = absorption(bc,bw)

#Standard Hamiltonain
H = δ*n + im*sqrt(γ/dt)*(adw-wda)

#Unitary
U = identityoperator(n) - im*dt*(δ*n + im*sqrt(γ/dt)*(adw-wda)) - 1/2*dt^2*( γ/dt*(adw*wda)  )

#Define unitary solver (veeeeeery simple)
function evolve_unitary(psi::Ket,U,times)
    out = copy(psi)
    tmp = copy(psi)
    for i in eachindex(times)
        set_waveguidetimeindex!(U,i)
        mul!(out,U,tmp,1,0)
        tmp.data .= out.data
    end
    out

end

#Define input onephoton state shape
σ = 1
t0 = 5
ξfun(t,σ,t0) = complex(sqrt(2/σ)* (log(2)/pi)^(1/4)*exp(-2*log(2)*(t-t0)^2/σ^2))
ξvec = ξfun.(times,σ,t0)

ψ_cw = onephoton(bw,ξvec)
psi = fockstate(bc,0) ⊗  ψ_cw 

#Solve
@btime ψ = waveguide_evolution(times, psi, H)
@btime ψ = evolve_unitary(psi,U,times)

This gives:

28.266 ms (134309 allocations: 5.53 MiB)
444.100 μs (2832 allocations: 120.38 KiB)

We are gaining like 60 orders of magnitude in performance, and this will be even more apparent for larger systems! To allow for an automatic expansion of the Hamiltonian, you mentioned some tricks we could use. How should we go about constructing the unitary automatically and symbolically?