# Stochastic Schrödinger equation

The stochastic Schrödinger equation has the basic form

\[i\hbar\frac{\mathrm{d}}{\mathrm{d} t} |\Psi(t)\rangle = \left(H + \sum_i H_i^s \xi_i(t)\right) |\Psi(t)\rangle\]

where $H_i^s$ are the terms in the Hamiltonian that are each proportional to an independent noise term $\xi_i(t)$. Note, that by default the noise is assumed to be real white-noise (uncorrelated in time), i.e. $\langle\xi_i(t)\xi_i(t')\rangle = \delta(t-t')$.

The above equation is straightforward to implement with

`stochastic.schroedinger(tspan, ψ0, H, Hs; dt=dt)`

Here, the first three arguments are the same as for `timeevolution.schroedinger`

. The additional argument `Hs`

is either an `AbstractOperator`

or a vector of operators. The first case corresponds to the stochastic Schrödinger equation featuring only a single noise process, while the latter corresponds to a number of noise processes, where the operators $H_i^s$ are the entries in the vector `Hs`

.

A time-dependent stochastic Schrödinger equation can also be implemented using `stochastic.schroedinger_dynamic`

. Instead of the two parts of the Hamiltonian `H`

and `Hs`

, the dynamic version takes two functions as arguments. One, that returns the deterministic part of the dynamic Hamiltonian, while the other returns a vector with the (time- or state-dependent) operators $H_i^s$ as entries.

```
function fdeterm(t, psi)
# Calculate time-dependent stuff
# Update deterministic part of Hamiltonian
H
end
function fstoch(t, psi)
# Calculate time-dependent stuff
# Update stochastic part(s) of Hamiltonian
Hs # This has to be a vector
end
stochastic.schroedinger_dynamic(tspan, ψ0, fdeterm, fstoch; dt=dt)
```

Note, that the solver requires to know the number of noise processes. This number is calculated automatically from the given function `fstoch`

by calculating the output at time `t=0`

. If you want to avoid an initial execution of the function, you can make the solver skip the initial calculation of `fstoch`

. This is done by passing the number of noise processes, i.e. the length of the vector `Hs`

that is returned by `fstoch`

, using the optional argument `noise_processes`

.

For some problems, it can be useful to renormalize the state vector after every time step taken (this can be used to avoid numerical issues for problems where the norm can become very small/large). To do this here, you can simply set the keyword `normalize_state=true`

, e.g.

`stochastic.schroedinger(tspan, ψ0, H, Hs; dt=dt, normalize_state=true)`

To be clear what this renormalization does, it uses a so-called `FunctionCallingCallback`

from **DifferentialEquations.jl's** callback library (see the documentation). So doing the following

```
norm_func(u::Vector{ComplexF64}, t::Float64, integrator) = normalize!(u)
ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func;
func_everystep=true, func_start=false)
stochastic.schroedinger(tspan, ψ0, H, Hs; dt=dt, callback=ncb)
```

is the same as setting `normalize_state=true`

. See also **DifferentialEquations.jl's** callback library.

## Detection schemes with the stochastic Schrödinger equation

Currently, **QuantumOptics.jl** features only one pre-defined measurement scheme to be used with the stochastic Schrödinger equation. Namely, with the function `stochastic.homodyne_carmichael`

, one can obtain the function `fdeterm`

and `fstoch`

needed to describe homodyne detection as derived by H. J. Carmichael. These can then be used with `stochastic.schroedinger_dynamic`

. Note, that this version does not preserve the norm.

Consider a system descibed by the Hamiltonian $H_0$. If it is subject to decay with an operator $C$ and the output from this decay channel is subject to homodyne detection, the problem at hand is completely described by the stochastic state-dependent Hamiltonian [1,2]

\[H = H_0 - \frac{i}{2}C^\dagger C + iC e^{-i\theta}\left(\xi(t) + \langle C e^{-i\theta} + C^\dagger e^{i\theta}\rangle\right).\]

Here, $\theta$ is the phase difference between the signal field (from the decay channel) and the local oscillator. In order to implement this, one needs to split this according to the deterministic and stochastic parts of the Hamiltonian. Additionally, the deterministic part needs to account for the nonlinear term (expectation value). This can be done by, for example,

```
θ = 0.5π
Hs = 1.0im*C*exp(-1.0im*θ)
Y = C*exp(-1.0im*θ) + dagger(C)*exp(1.0im*θ)
CdagC = -0.5im*dagger(C)*C
H_nl(ψ) = expect(Y, normalize(ψ))*Hs + CdagC
fdet_homodyne(t, ψ) = H0 + H_nl(ψ)
fst_homodyne(t, ψ) = [Hs]
stochastic.schroedinger_dynamic(tspan, ψ0, fdet_homodyne, fst_homodyne; dt=dt, normalize_state=true)
```

The above code would calculate a stochastic trajectory for the measurement of the y-quadrature ($\theta = \pi/2$). Since this scheme is quite commonly used, there is a function which defines the (performance optimized) functions `fdeterm`

and `fstoch`

. Hence, one can implement the above example in a much more simple way:

```
θ = 0.5π
fdet_cm, fst_cm = stochastic.homodyne_carmichael(H0, C, θ)
stochastic.schroedinger_dynamic(tspan, ψ0, fdet_cm, fst_cm; dt=dt, normalize_state=true)
```

## Functions

## Examples

## References

[1] Carmichael, Howard. An open systems approach to quantum optics: Lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991. Vol. 18. Springer Science & Business Media, 2009.

[2] Wiseman, H. M., & Milburn, G. J. (1993). Quantum theory of field-quadrature measurements. Physical review A, 47(1), 642.