Simulation problem for "Arrow of Time for Continuous Quantum Measurement"

Dear Prof.Justin Dressel, i am using QuantumBayesian.jl to reproduce the Fig. 2 in paper "Arrow of Time for Continuous Quantum Measurement" and encounter some problems.

I use $H=\Omega * \sigma_y$ and continuously measure the $\sigma_z$. I find that the statistics of the result trajectories does not agree with the Gaussian distribution given in this paper. I do not know where is my problem.

Here is the code i use

using QuantumBayesian
using QuantumOptics
using Statistics

q = SpinBasis(1//2)

z = spinup(q)
zm = spindown(q)
x = (z + zm)/sqrt(2)
Z = sigmaz(q) |> dense
Y = sigmay(q) |> dense
X = sigmax(q) |> dense

τ = 3.0        # Measurement collapse timescale
Ω  = 4π / τ    # Rabi frequency
T = (0.0, 2τ)  # Time duration of simulation;
η = 1.0        # Measurement efficiency

dt = 1e-3      # Simulation timestep (coarse to show method precision);

init = dm(x)

# Bloch expectation values
fs = collect(ρ -> expect(ρ, l)|>real for l in [X, Y, Z]);

# Trajectory function
function solve(inc, init, dt, T, fs; readout=false, n=1)
    if n==1
        out = trajectory(inc, init, (0.0, 0.1), fs..., dt=dt, points=10, verbose=false) #JIT throw-away
        @time out = trajectory(inc, init, T, fs..., dt=dt)
    else
        out = ensemble(1, inc, init, (0.0, 0.1), fs..., dt=dt, points=10, verbose=false) #JIT throw-away
        @time out = ensemble(n, inc, init, T, fs..., dt=dt)
    end

    return out
end;

# Hamiltonian
H(t) = Ω*Y
# Stochastic monitoring (unit efficiency)
SM = [(Z, τ, η)];

num_traj = 2000
out = solve(meas(dt, H, mclist=SM), init, dt, T, fs, readout=true, n=num_traj);
lnR = 2*dt/τ * sum(out[3][2:end-1,:] .* out[4][1:end,:], dims=1)[1,:];

In my code, i find both the mean and the variance of lnR disagree with the results in the paper $N(\frac{3T}{2\tau},\frac{2T}{\tau})$. I am not sure i do this all right?

Many thanks!

justindressel / QuantumBayesian.jl

Simulation problem for "Arrow of Time for Continuous Quantum Measurement" #13