Quantum trajectories

General Usage
1

Is there a way to implement the quantum trajectory algorithm described in https://qutip.readthedocs.io/en/latest/guide/dynamics/dynamics-monte.html
for a specific model using existing tools of differentialequations.jl without having to resort to dedicated packages?
I can use the ODEProblem to generate one realization, but I cannot find any way to generate an ensemble of distinct realizations. I always get the the same realization of the noise process: prob_func seems attuned to reset initial data and parameters but not a sequence of random numbers that need to be sampled at any time step. Any idea would be greatly appreciated.
Many thanks.

2

Hard to say anything about the noise process issues without seeing your code.

In general, for simulating quantum jump trajectories, I’d check out JumpProcesses.jl. The tutorial on PDMPs looks particularly relevant: Piecewise Deterministic Markov Processes and Jump Diffusion Equations · JumpProcesses.jl. As far as I can tell, it’s a simple wrapper around DifferentialEquations.jl that just makes it easier to set up the jump callback.

3

There is also Parallel Ensemble Simulations · DifferentialEquations.jl for ensembles of trajectories.

Without some code it will be hard diagnosing what doesn’t work for you.

4

Many thanks for your swift answer!

Yes I also used the available PDMPs algorithms, with alterante results. Sometimes they work very well sometines they don’t. My problem is not having a code that works for a specific problem in itself.

As I understand from the documentation, my question comes down to the fact that the underlying logic of the algorithms used in DifferentialEquations.jl and in the other toolboxes appears to be different. Is this that I would like to understand better.

The integration strategy followed by QuTip and (more or less) by QuantumOprics.jl is the following:

  1. Choose a random number r1 in [0,1] representing the probability that a quantum jump occurs.

  2. Choose a random number r2 in [0,1] used to select which collapse operator was responsible for the jump.

  3. Integrate the Schrödinger equation, using the effective Hamiltonian e.g.

G = i H + a C1 C1’ +b C2 C2’

until a time tau such that the norm of the wave function psi satisfies

sum(abs2, psi) = r1

at which time a jump occurs.

  1. The resultant jump projects the system at time tau into one of the renormalized states given by

psi_new = C1 \psi / sqrt(sum(abs2, C1 * psi) )
if
sum(abs2, C1 * psi) / (sum(abs2, C1 * psi) +sum(abs2, C2 * psi) ) >= r2
else
psi_new = C2 \psi / sqrt(sum(abs2, C2 * psi) )

5 Repeat until final time of the simulation

–
On a single realization I can implement this by introducing the condition 3
on the deterministic dissipative evolution. The condition triggers a continuous callback that implements 4.

Problem is how to pass r1 & r2 to the integrator. They must be independent random variables at every time step. Furthermore it must be possible to remake the simulation with different realizations of r1 & r2 for an arbitrary number of trajectories. A RODE algorithm with a custom noise seems perhaps the way to go but the problem is that the noise does not appear in the drift or the diffusion fields but must be only supplied to condition and callback.
Is it just me or does this mean that the exisisting tools cannot implement the QuTip algorithm?
Form this old discussion
[Code for simulating jump diffusion?]
it seems in fact that the DifferentialEquations algorithms have been developed having other applications in mind.

Many thanks!!!

5

sorry the boldfaced row was not meant.

6

Many thanks!

The question is more about the flowchart than a code.
With
prob_func
I can remake a simulation changing initial data or some parameters.

But the r1 and the r2 are neither of the two and do not enter the dynamics as noise process in the documentation. I tried several workarounds. I can share them but the point is that a-priori none of the workarounds corresponds to a documented flowchart.

7

I think you should just use a ContinuousCallback to perform the jumps. The Callback checks the norm and compares it to the threshold r1. If the jump occurs, then the callback projects the wavefunction and also computes a new r1. You can store this r1 in the parameters of the problem. I think this would be canonical way to implement this.

See Event Handling and Callback Functions · DifferentialEquations.jl

8

Let me see I understand your suggested flow chart.
Please bear with me I am a very dummy programmer.

Suppose I have

numf = number time steps

I should have an ODE with

#(parameters) = n1 + n2

n1 = other fixed paramater of the problem which are not relevant for us but need to be there.

n2 = 2 x numf because at each time step I need 2 independent random variables. Let’s denote the full collection
r1(1), …, r1(numf)
and
r2(1), …, r2(numf)

1 I create a function psi!(du,u,p,t)
2 I pass to to the integrator via p the n1 + n2 parameters. Suppose they are ordered as

p = [ 1,…,n1, r1(1), …, r1(numf) , r2(1), …, r2(numf)]

3 At each time step I check the condition on the wave function norm. So if I am at integration step
k
and r1 specifies the norm condition, the continuous callback is invoked if

sum(abs2,psi) == p[n1+k]

4 the continuous call back implements the first jump
(working hypotheses as in the previous post)

psi_new = C1 psi /sqrt(sum(abs2, C1 psi))

if sum(abs2, C1 psi) /(sum(abs2, C1 psi)+sum(abs2, C2 psi)) >= p[n1+numf+k]

5 repeat.

I am not sure whether this is what you are suggesting. I tried it but I encounter the following issues. Suppose I want to integrate up to tf and dt is the mesh i.e. dt = tf / numf

1 AFAIU at step k the integrator interface provides

integrator.t returning t = k dt

but no

“integrator.step” returning k

Now I wrote a code computing k from t but it looks to me an ugly workaround there must be a better way.

Prehaps a clean implementation is then to write an ODEProblem for each dt rather than a single one up to tf. In that case each of the ODE problems would have only

#(parameters) = n1 + 2

Of course it would be simpler had I something like

“integrator.step”

Suggestions?

2 Generation of the ensemble. I should have something like (unless I need to go for distinct ODEProblems for each dt)

function prob_func(prob,i,repeat)
x = rand!(zeros(2 * numf))
c = vcat [p[1], … p[numf], x ]
remake(prob; p = c)
end

to be passed to

ensemble_sse0 = EnsembleProblem(sse_prob0, prob_func=prob_func)

and then e.g.

sol_sse0 = DifferentialEquations.solve(ensemble_sse0, Vern7(), EnsembleSerial(), trajectories = num_traj);

This is the point that nonplusses me the most: I tried it many times but always ending up with the ensemble producing the same realization of the quantum trajectory.

I did not manage to pass correctly to condition and continuous callback the

p[n1+k], p[n1+numf+k]

parameters. Possibly the point is that I should define ODEProblems for each dt
and then remaking all of them one by one.

Any comment greatly appreciated.

Many Thanks!

9

Actually you will need 3 functions in total:

  1. the psi!(du, u, p, t) for computing the normal time evolution of the non-hermitian Hamiltonian
  2. the condition should_jump(u, t, integrator) that checks the norm of the wavefunction has decayed below r1
  3. the perform_jump!(integrator) that actually performs the quantum jump

I think you should generate these r1 and r2 on the fly since you don’t know how many quantum jumps will occur.

Not really: This happens inside the ContinuousCallback. NOT inside du!. DifferentialEquations.jl does this for you - you just set it up.

So a rough sketch of the program would be:

u0 = normalize!(rand(ComplexF64, 64)) # some wavefunction
tspan = (0.0, 15.0)
p = (; 
    jump_operators = [C1, C2, ...] # some jump operators
    H = I(64), # some non-hermitian Hamiltonian including the jump terms
    threshold = rand(), # norm to jump at
)

# wave function update
function psi!(du,u,p,t)
    mul!(du, p.H, u)
end

# jump condition
function should_jump(u, t, integrator)
    norm(u) - integrator.p.threshold
end

# actual jump happening
function perform_jump!(integrator)
    psi = integrator.u
    jump_operators = integrator.p.jump_operators 
    # now find the jump operator that should happen and apply it
    ....
    # set new jump threshhold
    integrator.p.threshold = rand()
end

prob = ODEProblem(psi!, u0, tspan, p)
cb = ContinuousCallback(should_jump, perform_jump!)
sol = solve(prob, Tsit5(), callback = cb)

Edit:
Bonus: Using the EnsembleProblem is quite easy. All that changes is the initial threshold. So you can just do:

prob = ODEProblem(psi!, u0, tspan, p)
cb = ContinuousCallback(should_jump, perform_jump!)
ensemble_prob = EnsembleProblem(prob;
    prob_func = (prob, i, repeat) -> (prob.p.threshold=rand(); prob))
sol = solve(prob, Tsit5(), callback = cb, trajectories=1000) # solve 1000 times

You should probably try and compute your observables for every trajectory instead of saving the full wavefunctions. You can do that by using a SavingCallback. But first get your code working before you add that :slight_smile:

10

Many thanks again. I really appreciate your help!

Your pseudocode is a much neater version what I tried. In particularly, I finally learnt how to use mutable structures in Julia thanks to you. Very grateful.

Unfortunately, once implemented also this flowchart seems to suffer the same problem that has been hauting me for a while. The ensemble ends up generating always the very same trajectory. Furthermore jumps do not occur as frequently as they should be and there is a probability loss.

So either proc_func or proc_func + remake do not do anything useful or there is a problem with the on the fly random number generation or both.

I’ll try to write down a cleaner implememtation with an even simpler model and post the code later.

Many thanks again

11

I think my condition is wrong. It should be norm(u)^2 - integrator.p.threshold. That ofc leads to less jumps.

12

I am not writing the condition in that way and still I get probability loss.
But that is not what concerns me the most.
What concerns me the most is not being able to generate an ensemble consisting of distinct realizations of quantum trajectories.
There is something I am missing on how EnsembleProblem works.
Many thanks.

13

Yes, if you could share a runnable example of a simulation that exhibits your problem, I think we should be able to get to the bottom of this pretty quickly. Sounds like some kind of seeding problem.

You can use a dead simple toy problem like du / dt = -u in place of your Hamiltonian. The important part is that you implement it in the same way you’ve implemented the full problem, with all the same functions, including randomly choosing between two different kinds of jumps every time.

14

Below one code that suffers the problem I was reporting. All realizations correspond to the same quantum trajectory. I must be doing something really silly …

numf = 500       # number of points
tspan=(0.0,tf)

 u0 = [1.0+ im * 0.0, 0.0 + im * 0.0]. # initial datum
par = [0.4,0.5,0.1].   # parameters

# Jump operators
  L1 = [0.0 + im * 0.0 1.0+ im * 0.0; 0.0 + im * 0.0 0.0+ im * 0.0]
  L2 = [0.0 + im * 0.0 0.0 + im * 0.0; 1.0+ im * 0.0 0.0+ im * 0.0]
  H = L1 * L2
  G = L2 * L1

 mutable struct Tuple2{}
    par::Vector{Float64}
    threshold::Float64
end
p = Tuple2(par, rand());

function psi!(du, u, p, t)
    Op = - im * p.par[1] * H  - p.par[2]/2 * G   - p.par[3] /2.0* H 
    mul!(du, Op, u)  
end

#The argument of condition must thus be a function which is zero at the point we want to trigger affect!

function jump_condition!(u,t,integrator)
    sum(abs2,integrator.u)-integrator.p.threshold
end  

function affect!(integrator)
  jump_prob1 = sum(abs2,L1*integrator.u)
  jump_prob2 = sum(abs2,L2*integrator.u)
  tot_jump_prob = jump_prob1 + jump_prob2
  if jump_prob1 / tot_jump_prob >= rand(Float64)
    jump_ampl1 = sqrt(jump_prob1) 
    mul!(integrator.u, L1 ,integrator.u / jump_ampl1 )
  else
     jump_ampl2 = sqrt(jump_prob2) 
     mul!(integrator.u, L2 ,integrator.u / jump_ampl2 )
  end
  integrator.p.threshold = rand(Float64)
end

jump_callback = ContinuousCallback(jump_condition!,affect!);

function prob_func(prob,i,repeat) 
  prob.p.threshold = rand()
   remake(prob; p = Tuple2(par, rand()))
end

num_traj=1000    # trajectories in the ensemble
sse_prob0 = ODEProblem(psi!, u0, (0.0, tf), p,  callback = jump_callback);
ensemble_sse0 = EnsembleProblem(sse_prob0; prob_func = prob_func)
sol_sse0 = DifferentialEquations.solve(ensemble_sse0, Vern7(), EnsembleSerial(), trajectories = num_traj)

Problem is that

sol_sse0[1](100)[1]. i.e. realization 1 at time 100 of u[1]
and e.g.
sol_sse0[100](100)[1]. realization 100 at time 100 of u[1]

are always the same!
Any help greatly appreciated!
Many thanks

15

Please enclose the code in triple backticks ``` like so

    ```
    code
    ```
16

It seems the question is mostly answered, but just FYI, posting in the “Specific Domain - Quantum” category can bring answers to this type of questions faster and ping more folks embedded in that ecosystem: Quantum - Julia Programming Language

17

Thanks done

18

Could it be that this is just caused by the toy example you consider? If I parse that right, then this is essentially a spin-half that randomly switches between up and down. You essentially don’t have any coherent dynamics, since your state will always be an eigenstate of the Hermitian part.
Maybe try this with a system that undergoes Rabi oscillations, i.e. choose the Hermitian part of the Hamiltonian as Sx = [0 1; 1 0] that should change things.

Edit 1: My explanation here is not quite correct. You should still get different states over time due to the random switching and the decay of the norm… I can try the code tomorrow.

Edit 2: Also your code has a few performance problems which don’t matter for this tiny system. I’ll gladly provide some comments once we got it correct if you’d like.

19

I am happy to hear any comments about performance.
But my scope here is really to understand how the algorithm and its Julia implementation work.
Many Thanks!

20

Here’s your problem: You’re referencing integrator.u, not the u that’s provided separately as an argument. This way you fail to ever trigger the callback, so your trajectory doesn’t have any jumps. You should define the condition as follows:

function jump_condition(u, t, integrator)
    return sum(abs2, u) - integrator.p.threshold
end

I highly recommend plotting as a debugging strategy here. It would make the lack of jumps immediately obvious and point you toward the location of the problem. In general, it’s good practice to do some sanity checks and plots of a single trajectory before moving on to ensembles.

Btw., the line prob.p.threshold = rand() in prob_func is redundant. The remake call is sufficient.

(Further explanation of the issue: when the callback mechanism sees a sign flip in the condition, it uses interpolation to locate the exact crossing. The interpolated state values are not saved to integrator.u, but only provided through the argument u, so a state-dependent callback must use the latter.)