How do I create a Quantum Time Operator with DifferentialEquations?

Specific Domains Quantum
1

For the special case of unitary matrices, my quantum mechanics books tells me to use Time Evolution Operators.
β(t) = exp(-im*H*t/h)*β(0)

As far I know, I just need to use Matrix Exponential.
First, decompose H = ψ*λ*inv(ψ)

Second, we create a time evolution operator U = exp(-im*H*t) = ψ*exp(-im*λ*t)*inv(ψ)

If one sets a small enough time interval δt : δU = exp(-im*H*δt) = ψ*exp(-im*λ*δt)*inv(ψ)

Then, I can evolve iteratively the solution βₜ₊₁ = δU*βₜ of my problem.

using LinearAlgebra
using Random
Random.seed!(1256)
N = 2

H = [1 1;-1 1]./sqrt(2)
β₀= rand(ComplexF64, N);

maxTimeSteps, t_max = 1000, 10
t = range(1, stop=t_max, length=maxTimeSteps);
δₜ = t[2]-t[1];

@time begin
    λ, ψ = eigen(H);
  
    δU = ψ*exp(-im*Diagonal(λ*δₜ))*inv(ψ)

    βₜ = []
    push!(βₜ,  deepcopy(β₀));
    for idx_t = 2:maxTimeSteps
        push!(βₜ, δU*βₜ[idx_t-1]  ) 
    end
end

I would like to solve this problem with DifferentialEquations, but I don’t know how to properly write down. It should be simple, but I’m taking time.
My bet was to use du/dt = a*u as input, because its solution is u(t)=exp(a*t)

using DifferentialEquations
function ff(du, u, p,t) 
    du[:] .= -im*H*u
end
tspan = (0.0,t_max);
prob = ODEProblem(ff,β₀,tspan);
@time sol = solve(prob, Tsit5(), reltol=δₜ, abstol=δₜ);

I can check the difference of the last results are far away of each other.

julia> abs.(βₜ[end] .- sol.u[end])
2-element Vector{Float64}:
 424.06945981919637
 424.06825305239903

Is there any one with physics background that already knows where my silly mistake is :tired_face:?

1 Like
2

See the examples which use the exponential integrators and Magnus methods.

2 Likes
3

Chris is right about Magnus propagators being the method of choice for QM. See also this earlier post: How to efficiently exponentiate a complex (symmetric) tridiagonal matrix?

Are you sure you want this Hamiltonian matrix though?

This is non-Hermitian, which is very untypical in Markovian (=normal) QM.

2 Likes
4

I knew it. I did some silly mistake, the wrong Hamiltonian was the first.
The second was to initialize the first method at t_0=1, and the other with t_0=0.

BUT I found major mistake :

I setup reltol=δₜ, abstol=δₜ, and in this example, δₜ=0.10010, which is a terrible tolerance. Just using default values, or not adaptive methods, everything worked as expected.

(...)
β₀ = [1.0+0im,1.0+0im] ./sqrt(2)

(...)
@time sol1 = solve(prob, Tsit5(), adaptive=true, dt=δₜ, dense=true, reltol=δₜ, abstol=δₜ);
@time sol2 = solve(prob, Tsit5(), adaptive=true, dt=δₜ, dense=true);
@time sol3 = solve(prob, Tsit5(), adaptive=false, dt=δₜ);


using Plots
fig1 = plot(t,  [abs2(βₜ[i][1]) for i=1:length(t)], label="exp")
plot!(t,  [abs2(sol1(t[i])[1]) for i=1:length(t)], label="sol1", lw=4, linestyle=:dot)

fig2 = plot(t,  [abs2(βₜ[i][1]) for i=1:length(t)], label="exp")
plot!(t,  [abs2(sol2(t[i])[1]) for i=1:length(t)], label="sol2", lw=4, linestyle=:dot)


fig3 = plot(t,  [abs2(βₜ[i][1]) for i=1:length(t)], label="exp")
plot!(t,  [abs2(sol3(t[i])[1]) for i=1:length(t)], label="sol3", lw=4, linestyle=:dot)

plot(fig1, fig2, fig3)

toDiscourse

Anyway, using the matrix exponential still faster approach if I don’t need to guarantee stability.

2 Likes