|
1 | 1 | # [Stochastic Solver](@id doc-TE:Stochastic-Solver) |
2 | 2 |
|
3 | | -This page is still under construction, please visit [API](@ref doc-API) first. |
| 3 | +When a quantum system is subjected to continuous measurement, through homodyne detection for example, it is possible to simulate the conditional quantum state using stochastic Schrödinger and master equations. The solution of these stochastic equations are quantum trajectories, which represent the conditioned evolution of the system given a specific measurement record. |
4 | 4 |
|
5 | | -## Stochastic Schrodinger equation |
| 5 | +## [Stochastic Schrödinger equation](@id doc-TE:Stochastic-Schrödinger-equation) |
| 6 | + |
| 7 | +The stochastic Schrödinger time evolution of a quantum system is defined by the following stochastic differential equation [Wiseman2009Quantum; section 4.4](@cite): |
| 8 | + |
| 9 | +```math |
| 10 | +d|\psi(t)\rangle = -i \hat{K} |\psi(t)\rangle dt + \sum_n \hat{M}_n |\psi(t)\rangle dW_n(t) |
| 11 | +``` |
| 12 | + |
| 13 | +where |
| 14 | + |
| 15 | +```math |
| 16 | +\hat{K} = \hat{H} + i \sum_n \left(\frac{e_n}{2} \hat{S}_n - \frac{1}{2} \hat{S}_n^\dagger \hat{S}_n - \frac{e_n^2}{8}\right), |
| 17 | +``` |
| 18 | +```math |
| 19 | +\hat{M}_n = \hat{S}_n - \frac{e_n}{2}, |
| 20 | +``` |
| 21 | +and |
| 22 | +```math |
| 23 | +e_n = \langle \psi(t) | \hat{S}_n + \hat{S}_n^\dagger | \psi(t) \rangle. |
| 24 | +``` |
| 25 | + |
| 26 | +Above, ``\hat{H}`` is the Hamiltonian, ``\hat{S}_n`` are the stochastic collapse operators, and ``dW_n(t)`` is the real Wiener increment (associated to ``\hat{S}_n``) which has the expectation values of ``E[dW_n]=0`` and ``E[dW_n^2]=dt``. |
| 27 | + |
| 28 | +The solver [`ssesolve`](@ref) will construct the operators ``\hat{K}`` and ``\hat{M}_n``. Once the user passes the Hamiltonian (``\hat{H}``) and the stochastic collapse operators list (`sc_ops`; ``\{\hat{S}_n\}_n``). As with the [`mcsolve`](@ref), the number of trajectories and the random number generator for the noise realization can be fixed using the arguments: `ntraj` and `rng`, respectively. |
| 29 | + |
| 30 | +## [Stochastic master equation](@id doc-TE:Stochastic-master-equation) |
| 31 | + |
| 32 | +When the initial state of the system is a density matrix ``\rho(0)``, or when additional loss channels are included, the stochastic master equation solver [`smesolve`](@ref) must be used. The stochastic master equation is given by [Wiseman2009Quantum; section 4.4](@cite): |
| 33 | + |
| 34 | +```math |
| 35 | +d \rho (t) = -i [\hat{H}, \rho(t)] dt + \sum_i \mathcal{D}[\hat{C}_i] \rho(t) dt + \sum_n \mathcal{D}[\hat{S}_n] \rho(t) dt + \sum_n \mathcal{H}[\hat{S}_n] \rho(t) dW_n(t), |
| 36 | +``` |
| 37 | + |
| 38 | +where |
| 39 | + |
| 40 | +```math |
| 41 | +\mathcal{D}[\hat{O}] \rho = \hat{O} \rho \hat{O}^\dagger - \frac{1}{2} \{\hat{O}^\dagger \hat{O}, \rho\}, |
| 42 | +``` |
| 43 | + |
| 44 | +is the Lindblad superoperator, and |
| 45 | + |
| 46 | +```math |
| 47 | +\mathcal{H}[\hat{O}] \rho = \hat{O} \rho + \rho \hat{O}^\dagger - \mathrm{Tr}[\hat{O} \rho + \rho \hat{O}^\dagger] \rho, |
| 48 | +``` |
| 49 | + |
| 50 | +The above implementation takes into account 2 types of collapse operators. ``\hat{C}_i`` (`c_ops`) represent the collapse operators related to pure dissipation, while ``\hat{S}_n`` (`sc_ops`) are the stochastic collapse operators. The first three terms on the right-hand side of the above equation is the deterministic part of the evolution which takes into account all operators ``\hat{C}_i`` and ``\hat{S}_n``. The last term (``\mathcal{H}[\hat{S}_n] \rho(t)``) is the stochastic part given solely by the operators ``\hat{S}_n``. |
| 51 | + |
| 52 | + |
| 53 | +## [Example: Homodyne detection](@id doc-TE:Example:Homodyne-detection) |
| 54 | + |
| 55 | +First, we solve the dynamics for an optical cavity at absolute zero (``0K``) whose output is monitored using homodyne detection. The cavity decay rate is given by ``\kappa`` and the ``\Delta`` is the cavity detuning with respect to the driving field. The homodyne current ``J_x`` is calculated using |
| 56 | + |
| 57 | +```math |
| 58 | +J_x = \langle \hat{x} \rangle + dW/dt, |
| 59 | +``` |
| 60 | + |
| 61 | +where ``\hat{x}`` is the operator build from the `sc_ops` as |
| 62 | + |
| 63 | +```math |
| 64 | +\hat{x} = \hat{S} + \hat{S}^\dagger |
| 65 | +``` |
| 66 | + |
| 67 | +```@setup stochastic-solve |
| 68 | +using QuantumToolbox |
| 69 | +
|
| 70 | +using CairoMakie |
| 71 | +CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) |
| 72 | +``` |
| 73 | + |
| 74 | +```@example stochastic-solve |
| 75 | +# parameters |
| 76 | +N = 20 # Fock space dimension |
| 77 | +Δ = 5 * 2 * π # cavity detuning |
| 78 | +κ = 2 # cavity decay rate |
| 79 | +α = 4 # intensity of initial state |
| 80 | +ntraj = 500 # number of trajectories |
| 81 | +
|
| 82 | +tlist = 0:0.0025:1 |
| 83 | +
|
| 84 | +# operators |
| 85 | +a = destroy(N) |
| 86 | +x = a + a' |
| 87 | +H = Δ * a' * a |
| 88 | +
|
| 89 | +# initial state |
| 90 | +ψ0 = coherent(N, √α) |
| 91 | +
|
| 92 | +# temperature with average of 0 excitations (absolute zero) |
| 93 | +n_th = 0 |
| 94 | +# c_ops = [√(κ * n_th) * a'] -> nothing |
| 95 | +sc_ops = [√(κ * (n_th + 1)) * a] |
| 96 | +``` |
| 97 | + |
| 98 | +In this case, there is no additional dissipation (`c_ops = nothing`), and thus, we can use the [`ssesolve`](@ref): |
| 99 | + |
| 100 | +```@example stochastic-solve |
| 101 | +sse_sol = ssesolve( |
| 102 | + H, |
| 103 | + ψ0, |
| 104 | + tlist, |
| 105 | + sc_ops, |
| 106 | + e_ops = [x], |
| 107 | + ntraj = ntraj, |
| 108 | +) |
| 109 | +
|
| 110 | +# plot by CairoMakie.jl |
| 111 | +fig = Figure(size = (500, 350)) |
| 112 | +ax = Axis(fig[1, 1], xlabel = "Time") |
| 113 | +#lines!(ax, tlist, real(sse_sol.xxxxxx), label = L"J_x", color = :red, linestyle = :solid) TODO: add this in the future |
| 114 | +lines!(ax, tlist, real(sse_sol.expect[1,:]), label = L"\langle x \rangle", color = :black, linestyle = :solid) |
| 115 | +
|
| 116 | +axislegend(ax, position = :rt) |
| 117 | +
|
| 118 | +fig |
| 119 | +``` |
| 120 | + |
| 121 | +Next, we consider the same model but at a finite temperature to demonstrate [`smesolve`](@ref): |
| 122 | + |
| 123 | +```@example stochastic-solve |
| 124 | +# temperature with average of 1 excitations |
| 125 | +n_th = 1 |
| 126 | +c_ops = [√(κ * n_th) * a'] |
| 127 | +sc_ops = [√(κ * (n_th + 1)) * a] |
| 128 | +
|
| 129 | +sme_sol = smesolve( |
| 130 | + H, |
| 131 | + ψ0, |
| 132 | + tlist, |
| 133 | + c_ops, |
| 134 | + sc_ops, |
| 135 | + e_ops = [x], |
| 136 | + ntraj = ntraj, |
| 137 | +) |
| 138 | +
|
| 139 | +# plot by CairoMakie.jl |
| 140 | +fig = Figure(size = (500, 350)) |
| 141 | +ax = Axis(fig[1, 1], xlabel = "Time") |
| 142 | +#lines!(ax, tlist, real(sme_sol.xxxxxx), label = L"J_x", color = :red, linestyle = :solid) TODO: add this in the future |
| 143 | +lines!(ax, tlist, real(sme_sol.expect[1,:]), label = L"\langle x \rangle", color = :black, linestyle = :solid) |
| 144 | +
|
| 145 | +axislegend(ax, position = :rt) |
| 146 | +
|
| 147 | +fig |
| 148 | +``` |
0 commit comments