question on mcsolver

[Albert]

I am trying to reproduce the Monte-carlo solver for a cat state (with two-body loss), do you know when I can find the source code for its implementation. 

1. Also, do you know what is the effective implementation of the following (so far, I am still using discretized time step, and check this value in every time step.)

"I: Integrate the Schrödinger equation, using the effective Hamiltonian (1) until a time 

𝜏 such that the norm of the wave function satisfies ⟨𝜓(𝜏)∣∣𝜓(𝜏)⟩=𝑟1

, at which point a jump occurs."

2. In this final time, the wavefunction is not guaranteed to normalized, right? then, should I normalize it?

[Albert]

In this paper, the paragraph below Eq. (11) mention that (https://arxiv.org/pdf/1907.07079.pdf), do you know what exactly is this method? References are greatly appreciated.

"For instance, the popular QuTiP library (Johansson et al., 2012, 2013) uses a logarithmic secant method to numerically solve the equation ⟨ψi(t)|ψi(t)⟩ = r for the time t."

[Éric Giguère]

Qutip solve until it the norm of the wavefunction is smaller than the jump target, then it go back to find the exact time of the jump.

To estimate the location of the jump, it a variant of the bisection method with estimation using logarithm since the norm is expected to go as `~exp(-at)`

https://github.com/qutip/qutip/blob/e010ce4ec150d19cca1265e145b17eaaf6d5480c/qutip/cy/mcsolve.pyx#L186

The state is renormalize when used to compute the density matrix  and at jumps.

[Albert]

Thanks a lot for your quick reply. I now understand the method of "logarithmic secant" method. However, I am still confused by the underlying logic here. It would be great if you can kindly illustrate it:

1. For example, starting from t0 to tf.

2. Now we want to find the first point \tau, such that <\psi(\tau)|\psi(\tau)>=r_1.

3. The secant method need to point. But the available point is only t0. How to get the second point here? -- if a small time step needed here? and how to decide this time step.

[Éric Giguère]

Say we are a t_n. A step to t_n+1is done.

If <\psi(\tau)|\psi(\tau)><r_1, then it use the secant method to find the jump.

The dt between t_n, t_n+1 is set by the integration method.

Many ODE methods have an `dense_output`, which makes going back more efficient. Qutip use one such dense output interval for dt.

[Albert]

Thanks for your quick response!

Can I understand the logic as follow:

1. Divide t0 to tf to N interva, i.e., t0, t1, t2, ....

2. Example: t0 -> t1: Evaluate norm at t_0 and t_1 under non-hermitian Schrödinger equation. If norm at t1 is small then r_1, then we search for the quantum jump point. 

3. Repeat this process for other interval. 

Is this the way qutip implement the mcsolve? If not, can you correct me where I misunderstood?

[Éric Giguère]

Pretty much.

1. The division is not done all at once, but at each step. dt can varies for ODE method that support variable step length.

2.5 When a jump is found, the state is updated according to it.

[Albert]

Thanks a lot for the clarification. Can you explain more on how the following step is implemented, 

"The division is not done all at once, but at each step. dt can varies for ODE method that support variable step length."

So far, I have implemented the method I outlined before, but it is pretty unstable, and slow. Hence, it would be great if you can explain this key step to me.

[Éric Giguère]

That one is not easy, I would suggest to read qutip's source...

It depend on the ODE solver you use. qutip use the scipy.integrate.ode with zvode, it use old but efficient Fortran code and it compute the dt.

But the dense step output is not available from scipy's interface and you need to study scipy's code a little.

Using the new solve_ivp, you could try `events` to find the jump. I have no idea how efficient it is.