Stochastic Differential Equations, Quantum Phase Space, and Julia¶

$$ \def\julia{\texttt{julia}} \def\dashint{{\int\!\!\!\!\!\!-\,}} \def\infdashint{\dashint_{\!\!\!-\infty}^{\,\infty}} \def\D{\,{\rm d}} \def\E{{\rm e}} \def\dx{\D x} \def\dt{\D t} \def\dz{\D z} \def\C{{\mathbb C}} \def\R{{\mathbb R}} \def\CC{{\cal C}} \def\HH{{\cal H}} \def\I{{\rm i}} \def\qqqquad{\qquad\qquad} \def\qqfor{\qquad\hbox{for}\qquad} \def\qqwhere{\qquad\hbox{where}\qquad} \def\Res_#1{\underset{#1}{\rm Res}}\, \def\sech{{\rm sech}\,} \def\vc#1{{\mathbf #1}} $$

Dr. Ashton Bradley

ashton.bradley@otago.ac.nz
http://amoqt.otago.ac.nz

Lecture 2: SDE's and Quantum Phase Space¶

• Weiner increments
• Stochastic differential equations
• Introduction to phase space mappings for bosons

References

• An algorithmic introduction to numerical simulation of stochastic differential equations, D. J. Higham, SIAM Review, 43, 525-546 (2001)
• The Quantum World of Ultracold Atoms and Light: Book I: Foundations of Quantum Optics, C. W. Gardiner and P. Zoller, Imperial College Press, London (2014)
• Quantum Noise - A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics, C. W. Gardiner, P. Zoller, Springer-Verlag Berlin Heidelberg (2004)
• Quantum Optics, D. F. Walls, G. J. Milburn, Springer-Verlag Berlin Heidelberg (2008)

Wiener paths¶

We can sample a discretization of the Wiener process, built up from Wiener increments:

Our first aim is to see how to construct such individual paths, $W(t)$.

Wiener process¶

We start with an idealized Gaussian Langevin noise $\xi(t)$, with moments $\langle\xi(t)\rangle=0$, and $\langle \xi(t)\xi(t')\rangle=\delta(t-t')$.

This is sometimes referred to as white noise since the correlation spectrum involves all frequencies:

$$\langle \xi(t)\xi(t')\rangle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty d\omega\; e^{i\omega (t-t')}.$$

We can view this as a useful idealization of a physical noise with a (large) spectral bandwidth.

We use $\xi(t)$ to define a Wiener process as

The Wiener process has moments \begin{align} \langle W(t)\rangle&=\int_0^t\langle\xi(t')\rangle = 0,\\ \langle W(t)^2\rangle&=\int_0^t dt'\int_0^t dt'\langle\xi(t')\xi(t)\rangle = t, \end{align}

so $W(t)$ is Gaussian, with mean zero and variance $t$, and probability density

\begin{align} p(w,t|0,0)dw &\equiv \textrm{Prob.} \quad w<W(t)<w+dw\notag\\ &=\frac{1}{\sqrt{2\pi t}}e^{-w^2/2t}dw \tag{Diffusion} \end{align}

$W(t)$ is the fundamental Gaussian diffusion process from which all others in these lectures are constructed.

Wiener increments¶

Consder the time interval $[t_0,t_N]$, broken into discrete times $t_0<t_1<t_2<\dots <t_{n-1}< t_n<\dots t_N$, and define $\Delta t_n\equiv t_n-t_{n-1}=n\Delta t$, where $\Delta t=t_0+(t_N-t_0)/N$, and $n=0,1,\dots, N$. The Wiener increments are defined as

\begin{align} \Delta W_n &\equiv W(t_n)-W(t_{n-1})\notag\\ &=\int_{t_{n-1}}^{t_n}dt'\;\xi(t') \end{align}

The form of $p(w,t|0,0)$ shows that the $\Delta W_n$ are all statistically independent, with probability density

$$ p(\Delta w_n)=\frac{1}{\sqrt{2\pi\Delta t }}\exp{\left(-\frac{(\Delta w_n)^2}{2\Delta t }\right)}.$$

Hence they are independent Gaussian random numbers with zero mean and varance $\Delta t$:

Note that this means $\Delta W_n\sim \sqrt{\Delta t}$, and hence a small differential (Wiener increment) cannot be understood in terms of regular calculus (e.g. standard Taylor expansion).

Langevin equation¶

Consider the fairly general form the the Langevin equation

$$\frac{dx(t)}{dt}=a(x(t),t)+b(x(t),t)\xi(t).$$

How can we solve this numerically? It turns out that the way we discretize this equation will fundamentally alter the numerical solution. Thus it appears that there is no unique solution. However, stochastic calculus saves us here by making a rigorous connection between different choices of discretization and the FPE; since the FPE is our link to a physical system, we will always find a well-defined numerical problem to solve.

SDE's according to Ito¶

Discrete form¶

Using the Wiener increments, we define a discretized SDE for an unknown field $x(t)$ in Ito form as follows:

and we thus set things up to give separate treatment to the regular $\Delta t$ term (an ODE term that will turn out to be drift in the FPE), and the $\sqrt{\Delta t}$ term (diffusion in the FPE).

In discretization, the terms in the Langevin equation are evaluated in the following way:

\begin{align} a_n \equiv a(x_n,t_n)\notag\\ b_n \equiv b(x_n,t_n)\notag \end{align}

Importantly, in the Ito SDE, the coefficients for the $x_{n+1}$ equation only depend on the value of the field at the previous time, $x_n$.

Ito stochastic integral¶

We define Ito stochastic integration with respect to the Wiener increments as

where the crucial feature of the Ito form is that the noise is non-anticipating: The noises $\Delta W_j$ and the integrand $G(t_{j-1})$ are independent over each small time increment $\Delta t$. Here we have introduced the notation $(I)$ to denote an Ito integral.

This definition is precisely the choice of discretization appearing in the discrete Ito SDE. Provided the SDE converges, it may be written as the ___Ito stochastic differential equation___

a shorthand for the solution that we may write as the ___Ito stochastic integral equation___

Terminology

• When $b[x(t),t]\equiv b$, the noise term is referred to as ___additive noise___.
• When $b[x(t),t]$ depends on $x(t)$, the noise term is referred to as ___multiplicative noise___.

Example: Integrate the pure diffusion SDE

$$dx(t) = dW(t)$$

for $x(t_0=0)=0$.

We find $x(t) = (I)\int_{0}^t dW(t')=\lim_{n\to\infty}\sum_{j=1}^n\Delta W_j=W(t)-W(0)$. Since $W(0)=0$, we have

$$x(t)=\int_0^t dW(t)=W(t),$$

the Wiener process.

Ito calculus¶

Since the discrete increments $\Delta W_n$ are always independent of the field $x(t)$, the rules for manipulating Ito SDE's can be reduced to a rather simple Ito calculus: there are two infinitesimals $dt$, and $dW(t)$. As shown above, $dW^2$ is of order $dt$. The rules of ordinary calculus are now supplemented by Ito rules:

in any manipulation, we must keep all terms up to order $dt$, accounting for the fact that $dW$ is also of order $\sqrt{dt}$. Thus, a consistent Taylor expansion to first order in $dt$ requires all second order terms in $dW$ to be accounted for.

Example: $dW dt=O(dt^{3/2})$ would not appear at first order in $dt$.

Independence of Ito noise¶

Since $x(t)$ only depends on $W(t')$ for $t'<t$, the increment $dW(t)$ is statistically independent of $f(x(t))$ for any function $f$. Hence

$$\langle f(x(t))dW(t)\rangle = \langle f(x(t))\rangle \langle dW(t)\rangle = 0$$

Change of variables and FPE¶

To see how this all fits together, we can now use Ito rules to carry out a change of variables and understand the link between the Ito SDE and the Fokker-Planck equation.

Ito's change of variables¶

Given $x(t)$ that obeys the Ito SDE, for any function $f(x)$, we can Taylor expand as usual

$$df(x)=f'(x)dx + \tfrac{1}{2}f''(x)dx^2+\dots,$$

but now to find a new SDE for $df$, we substitute the SDE and use Ito's rules, keeping terms up to order $dt$, to give ___Ito's forumula:___

The second part of the new drift term is called the Ito correction, and came from keeping terms of order $dW^2$.

Exercise¶

Consisder the Brownian motion equation introduced by Langevin (last lecture). Using an Ito interpretation of the noise (choice of discretization), the Langevin equation becomes the Ito SDE

$$(I)du = \underbrace{-\gamma u}_{a(u,t)}dt+\underbrace{\sqrt{f}}_{b}dW(t).$$

Work out the mean and variance for $u(t)$, the velocity of the Brownian particle.

Solution:

First, we take the average over realizations of the noise, acting on both sides of the equation

$$\langle du\rangle = d\langle u\rangle = -\gamma\langle u\rangle dt,\tag{$\langle dW(t)\rangle=0$}$$

and so $\langle u(t)\rangle = \langle u(0)\rangle e^{-\gamma t}$; the initial velocity is damped by the frictional force.

Second, we can change variables using the function $g(u)=u^2$ to get a new SDE for $dg(u)=d(u^2)$. Using Ito's formula, we have $g'=2u$, $g''=2$. Identifying $a$ and $b$, Ito's formula reads

$$ dg=\left[g'(u)(-\gamma u)+\tfrac{1}{2}g''(u)(\sqrt{f})^2\right]dt + g'(u)\sqrt{f}dW(t)$$

or

$$d(u^2)=(-2\gamma u + f)dt + 2u\sqrt{f}dW(t)$$

and we see an important "Ito correction" ($f$) appearing in the drift term due to the noise. Taking the noise-average, we then have an ODE

$$d\langle u^2\rangle=(-2\gamma\langle u^2\rangle +f)dt$$

and the solution

$$\langle u(t)^2\rangle = \langle u(0)^2\rangle e^{-2\gamma t}+\frac{f}{2\gamma}.$$

The steady state $\langle u^2\rangle_s=f/2\gamma$ depends crucially on the Ito correction. In this case the result we have found using Ito calculus (a particular choice of discretization of the the Langevin equation) agrees with our formal integration of the Langevin equation (last lecture). There is an interesting reason behind this: for additive noise, the result is insensitive to the particular discretization choice. Shortly we will discuss an alternative discretization due to Stratonovich [denoted $(S)$], and you should be aware that for multiplicative noise the two SDEs give different results [yes, there is a way to map between $(I)$ and $(S)$].

Fokker-Planck equation¶

We can now understand a fundamental link between SDE and FPE.

The solution $x(t)$ of the SDE is a Markov process: only the initial condition and functional form of $a[x(t),t]$ and $b[x(t),t]$ are required to determine the solution. We can hence describe the solution with a conditional probability $p(x,t|x_0,t_0)$. Consider now the time evolution of an arbitrary function $f(x):$

$$\frac{d}{dt}\langle f(x)\rangle = \int dx\;f(x)\frac{\partial}{\partial t}p(x,t|x_0,t_0)$$

According to Ito's formula

\begin{align} \frac{d}{dt}\langle f(x)\rangle &= \frac{\langle df(x)\rangle}{dt} = \langle a(x,t)f'(x)+\tfrac{1}{2}b(x,t)^2f''(x)\rangle\notag\\ &=\int dx\; \left\{ a(x,t)f'(x)+\tfrac{1}{2}b(x,t)^2f''(x)\right\}p(x,t|x_0,t_0)\notag\\ &=\int dx\;\left\{ -\frac{\partial}{\partial x}\left(a(x,t)p(x,t|x_0,t_0)\right)+\tfrac{1}{2}\frac{\partial^2}{\partial x^2}\left(b(x,t)^2p(x,t|x_0,t_0)\right)\right\}f(x) \end{align}

Since $f(x)$ is arbitrary, we can equate the previous equation with our last expression and find the ___Fokker-Planck equation___

is equivalent to a stochastic process governed by the ___Ito SDE___

SDE's according to Stratonovich¶

For completeness, let us briefly link to a different version of stochastic integration due to Stratonovich.

Stratonovich proposed an alternate dicretization of the Langevin equation:

based upon a Runge-Kutta midpoint rule. The stochastic integral is now defined as

$$(S)\int_{t_0}^t G(t)dW(t)\equiv \lim_{n\to\infty}\sum_{j=1}^n \tfrac{1}{2}\left[G(t_{j-1})+G(t_j)\right]\Delta W_j\tag{Stratonovich},$$

where the crucial point here is the appearance of $G(t_j)$ in the integrand (or $b_{n+1}$ in the discretization): the noise and the integrand are no longer statistically independent! The solution for $G(t)$ must be found by some form of implicit method, self-consistent with a particular realization of the noise.

This difference has the effect (that we make no effort to show!) that the SDE

obeys the chain rule of ordinary calculus for change of variables

$$df(x)=f'(x)dx,$$

and for a Stratonovich SDE solution $\langle x(t)dW(t)\rangle\neq 0$.

The main advantage of the Stratonovich form is that it is more stable to numerically integrate.

The Stratonovich SDE often arises from a physical system as the limit of a non-white noise, for which there is then an equivalent Ito SDE. [See Stochastic Methods for rules for converting between $(I)$ and $(S)$ forms of an SDE].

Summary¶

Ito

• Useful analytically since field and noise are independent
• Must take account of $dW^2$ terms, at order $dt$

Stratonovich

• Often arises physically as the limit of a non-white noise process
• More stable to integrate numerically than Ito.
• Differentials obey regular chain rule.

Quantum phase space¶

Our aim: mapping Bose field dynamics to stochastic diffusion.

Why? because Hilbert space for bosons gets very large very quickly unless we can utilize an efficient basis.

Open Quantum Systems¶

Bose-Einstein condensates¶

Stochastic Projected Gross-Pitaevskii Equation (truncated Wigner): 1D harmonic trap.

Our approach is to use the coherent state basis to effect the mapping

$$\dot{\rho}(t)\longrightarrow \partial_tP(\alpha,\alpha^*,t)\longrightarrow d\alpha(t)$$

and then solve the resulting SDE, enabling calculation of operator averages via averages over ensembles of independent paths of the SDE.

When the SDE formulation is available, it has a much better scaling with system size than the full Hilbert space.

Key ideas¶

Consider a single mode described by annhilation operator $a$, with non-vanishing Bose commutation relation $[a,a^\dagger]=1$. The coherent states, eigenstates of annihilation $a|\alpha\rangle = \alpha|\alpha\rangle$, with complex eigenvalue $\alpha=\alpha_r+i\alpha_i$ and useful expansion in number states

$$|\alpha\rangle = e^{-|\alpha|^2/2}\sum_{n=0}^\infty \frac{\alpha^n}{\sqrt{n!}}|n\rangle,$$

are overcomplete:

$$\mathbb{1}=\frac{1}{\pi}\int d^2\alpha \;|\alpha\rangle\langle\alpha|$$

Exercise: If you haven't done this before, prove it by using the number state expansion and evaluating the integral $d^2\alpha\equiv d\alpha_rd\alpha_i$.

Note there is no right eigestate of $a^\dagger$, however, the action of $a^\dagger$ on $|\alpha\rangle$ can be written in terms of a differential operator.

We will focus on the simplest example of the mapping to phase space, using the Glauber-Sudarshan $P$-function. There are other phase-space distributions known as the Wigner function ($W$), and the Husimi-$Q$ function.

P-function¶

We expand the density matrix $\rho(t)$ describing an open quantum system in terms of a distribution $P(\alpha,\alpha^*,t)$. We suppress time for now since the expansion applies at each instant:

Provide such an expansion exists, we can try to use it to map the quantum dynamics

$$\dot{\rho}=\underbrace{-\frac{i}{\hbar}[H,\rho]}_{\textrm{Hamiltonian} = H(a,a^\dagger)}+\underbrace{{\cal L}(\rho)}_{\textrm{Lindblad dissipator}}$$

to an equation of motion for $P(\alpha,\alpha^*)$, and then to a stochastic differential equation.

Operator averages¶

First we should ask: what information is available from the $P$-function?

Using properties of the trace, the $P$-fuction gives:

\begin{align} \langle (a^\dagger)^na^m\rangle &= \textrm{tr}\left[\rho(a^\dagger)^na^m\right]=\textrm{tr}\int d^2\alpha\;|\alpha\rangle\langle\alpha|P(\alpha,\alpha^*)(a^\dagger)^na^m\notag\\ &=\int d^2\alpha\;P(\alpha,\alpha^*)(\alpha^*)^n\alpha^m\equiv \overline{((\alpha^*)^n\alpha^m)}_P. \end{align}

Hence we find that moments of $P$ give normally ordered operator averages:

and this will later give us a very clear link to ensemble averages over the paths of an equivalent SDE.

Operator mappings¶

Based upon the form of the expansion, we can focus on the action of operators on $|\alpha\rangle\langle\alpha|$. For example $a|\alpha\rangle\langle\alpha|=\alpha|\alpha\rangle\langle\alpha|$, and hence

$$a\rho\longrightarrow \alpha P$$

While this first mapping is easy, a bit more work is needed to show (you should verify this):

$$a^\dagger|\alpha\rangle\langle\alpha|=\left(\alpha^*+\frac{\partial}{\partial \alpha}\right)|\alpha\rangle\langle\alpha|.$$

Note that $\alpha$ and $\alpha^*$ are formally independent, representing two independent real degrees of freedom $\alpha_r$, $\alpha_i$. Pluggin the result into the definition, and integrating by parts to shift the differential operator onto $P$, we find the mapping

$$a^\dagger\rho\longrightarrow \left(\alpha^*-\frac{\partial}{\partial \alpha}\right) P.$$

The other mappings follow similarly.

Summary¶

We can summarize the $P$-function mappings as

In these lectures we will see how to use these mappings to solve problems in Quantum Optics as stochastic diffusion problems.

We can now use the coherent state basis to effect the mapping

$$\dot{\rho}(t)\longrightarrow \partial_tP(\alpha,\alpha^*,t)\longrightarrow d\alpha(t)$$

and then solve the resulting SDE, enabling calculation of operator averages via averages over ensembles of independent paths of the SDE.

There are limitations to such an approach:

• For highly nonlinear systems, the mapping to $P$ may lead to 3rd order or higher derivatives, for which the mapping to SDE cannot be carried out.

• In some cases, a strict FPE cannot be mapped to an SDE, if there is no equivalent diffusion. Formally the SDE mapping requires the diffusion matrix to be positive semidefinite.

Coherent state basis: from P to +P¶

If we use the completeness relation to change bases in the usual way, we find

$$\rho=\mathbb{1}\rho\mathbb{1}=\frac{1}{\pi}\int d^2\alpha|\alpha\rangle\langle \alpha|\rho \frac{1}{\pi}\int d^2\beta\;|\beta\rangle\langle\beta|\\ =\int d^2\alpha\int d^2\beta\;\frac{\langle \alpha|\rho|\beta\rangle}{\pi^2}|\alpha\rangle\langle\beta|$$

Or in other words, the $P$ function is a special case of this expansion, chosen to be diagonal in coherent states.

• The $P$-function has limitations that mean it is not always positive, or non-singular.
• The doubled-phase space expansion leads to the $+P$-representation. This may always be chosen positive.
• Some caution is still required when using the $+P$, since it has its own limitations; this is a large topic in itself. As a simple heuristic, the practicality of the $+P$ usually is determined by the the relative strength of damping to nonlinearity: ___sufficient damping stabilizes the $+P$___.

Next: Workshop 1¶

• Introduction to $\julia$