Introduction #

Many systems are described by second-order Hamiltonians. These have terms linear and quadratic in the operators $\hat{x},\hat{p}$, for example:

$\hat{H}=a\hat{p}^2+b\hat{x}\hat{p}+c\hat{x}.$

One such example is the harmonic oscillator: $\hat{H}=\frac{\hat{p}^2}{2m}+\frac{1}{2}m\omega^2\hat{x}^2$. This describes a particle on a spring, in which case $\hat{x},\hat{p}$ represent the position and momentum.

In quantum optics we use a harmonic oscillator to represent light. The frequency of oscillation is the frequency of the optical field, and the number of vibrational quanta is the number of photons. We then view $\hat{x},\hat{p}$ as two non-commuting continuous variables describing the state of the optical field. They do not represent the position or momentum of a photon, which may not be well defined, since the field could contain an indefinite number of photons with indefinite positions and momenta. Emphasising this somewhat abstract nature, we refer to $\hat{x},\hat{p}$ as quadratures.

We are generally interested in multiple systems interacting with each other. In quantum optics, these may be two different modes of light, or the fields of two coupled cavities. If we index the system with a subscript, a second-order Hamiltonian my look something like

$\hat{H}=a\hat{p}_1\hat{p}_2+b\hat{x}_1^2+c\hat{x}_2.$

The evolution of our systems will depend on if the system obeys fermionic or bosonic statistics. We will consider bosons here, where different sub-systems are allowed to be in the same state. Light is a boson, which is why you can have multiple photons confined to the same cavity mode.

In this post we will show that a general second-order bosonic Hamiltonian can be described by a real symmetric matrix, and a real vector. This is useful for many reasons:

• We can reduce questions about operators to questions in matrix algebra. For example, diagonalising a general Hamiltonian becomes a problem of diagonalising a symmetric matrix.
• We can better understand how many degrees of freedom there are in a general Hamiltonian, and find transformations which make the equations easier to solve.
• Evolution in the infinite-dimensional Hilbert-space of the harmonic oscillator can be replaced with a differential equation in terms of real vectors and matrices.

Second-order Hamiltonians #

A first-order Hamiltonian is one which is linear in $\hat{x},\hat{p}^{\dagger}$, for example:

$\hat{H}=a\hat{x}.$

Note that for $\hat{H}$ to be Hermitian, $a$ must be real. Looking at the Heisenberg equation of motion for $\hat{p}$ will show that this describes a constant force on the system. Unfortunately there isn't a lot that we can describe using only linear terms, so a bit more complexity is required.

Quadratic Hamiltonians contain only squares of terms, for example:

$\hat{H}=\frac{\hat{p}_1^2}{2m}+\frac{\hat{p}_2^2}{2m}+\frac{1}{2}m\omega^2(\hat{x}_1-\hat{x}_2)^2.$

This describes two point masses connected by a spring of frequency $\omega$. A quadratic Hamiltonian can described two coupled oscillators, but not a constant force acting on one of them.

A second-order Hamiltonian contains both linear and quadratic terms:

$\hat{H}=\frac{\hat{p}_1^2}{2m}+\frac{\hat{p}_2^2}{2m}+\frac{1}{2}m\omega^2(\hat{x}_1-\hat{x}_2)^2+a\hat{x}_1.$

While the Hamiltonian contains quadratic terms, the resulting Heisenberg equations of motion for $\hat{x},\hat{p}$ will be linear. Describing nonlinear phenomena is thus still out of reach.

Parameterisation #

In this section we will follow [§3.1,Ser17] and introduce our parameterisation. Let us begin with a single mode, and define

$\hat{\mathbf{r}}=\left(\begin{array}{c}\hat{x} \\ \hat{p}\end{array}\right).$

Using the relation $[\hat{x},\hat{p}]=i$ (this is where the 'bosonic' assumption comes in), this has commutator:

$[\hat{\mathbf{r}},\hat{\mathbf{r}}^T]=\left(\begin{array}{cc}[\hat{x},\hat{x}] & [\hat{x},\hat{p}] \\ [\hat{p},\hat{x}] & [\hat{p},\hat{p}]\end{array}\right)=\left(\begin{array}{cc} 0 & i \\ -i & 0 \end{array}\right).$

Motivated by this we will define the matrix

$\Omega_1=\left(\begin{array}{cc} 0 & 1 \\ -1 & 0\end{array}\right),$

in terms of which the commutator becomes

$[\hat{\mathbf{r}},\hat{\mathbf{r}}^T]=i\Omega_1.$

In general we will have multiple modes $\hat{\mathbf{r}}=(\hat{x}_1,\hat{p}_1,\hat{x}_2,\hat{p}_2,\ldots,\hat{x}_n,\hat{p}_n)^T$ (with the transpose 'T' indicating this is a column vector). Using $[\hat{x}_j,\hat{p}_k]=i\delta_{jk}$, the commutator will now be

$[\hat{\mathbf{r}},\hat{\mathbf{r}}^T]=i\Omega,$

where

$\Omega=\oplus_{j=1}^n\Omega_1.$

The symbol '$\oplus$' denotes the 'direct sum', a cousin of the tensor product. If $\{A_i\}$ is a collection of matrices, the direct sum is the new matrix formed by placing these along the diagonal:

$\oplus_{j=1}^nA_j=\left(\begin{array}{ccc} A_1 & 0 & 0 \\ 0 & A_2 & 0 \\ 0 & 0 & \ddots\end{array}\right).$

Then for example

$\begin{aligned}\oplus_{j=1}^3\Omega_1 &= \Omega_1\oplus\Omega_1\oplus\Omega_1, \\ &=\left(\begin{array}{cccccc} 0 & 1 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & -1 & 0 \end{array}\right).\end{aligned}$

In terms of this $\hat{\mathbf{r}}$, a general second-order Hamiltonian $\hat{H}$ can be decomposed as

$\hat{H}=\frac{1}{2}\hat{\mathbf{r}}^TH\hat{\mathbf{r}}+\hat{\mathbf{r}}^T\vec{r},$

• $H$ is a real, symmetric matrix.
• $\vec{r}$ is a vector of real numbers.

We can see from the above equation that for $\hat{H}$ to be Hermitian, $H$ must be Hermitian and $\vec{r}$ must be real. We will see in the next section however that any imaginary parts of $H$ cancel out and don't affect the Hamiltonian, so we may as well take $H$ to be real. The commutation relations that we worked out will turn out to be helpful later on.

Some examples #

Let's go through some examples to get some intuition about $H$ and $\vec{r}$ mean.

A single mode #

The simplest case is just a single mode. Letting $H$ have elements $h_{jk}$ and $\vec{r}$ elements $r_j$, the Hamiltonian expands as

$\hat{H}=\frac{1}{2}\left(h_{11}\hat{x}^2+h_{12}\hat{x}\hat{p}+h_{21}\hat{p}\hat{x}+h_{22}\hat{p}^2\right)+r_1\hat{x}+r_2\hat{p}.$

For $\hat{H}$ to be Hermitian, we must have $h_{11},h_{22}$ real and $h_{21}=h_{12}^*$. The cross terms then simplify as

$\begin{aligned}h_{12}\hat{x}\hat{p}+h_{12}^*\hat{p}\hat{x}=\left(h_{12}+h_{12}^*\right)\hat{p}\hat{x}+ih_{12}, \\ =2\mathrm{Re}(h_{12})\hat{p}\hat{x}+ih_{12},\end{aligned}$

where $\mathrm{Re}$ denotes the real part. Any complex components of $h_{12}$ thus only show up as a constant, which doesn't affect the dynamics of the system. For this reason we can just take $H$ to be a real symmetric matrix without any loss of generality.

The Heisenberg equations of motion are then

$\begin{aligned}\hat{x}' = i[\hat{H},\hat{x}] &= h_{12}\hat{x}+h_{22}\hat{p}+r_2,\\ \hat{p}' = i[\hat{H},\hat{p}] &= -h_{11}\hat{x}-h_{12}\hat{p}-r_1.\end{aligned}$

• $\vec{r}$ gives a constant driving to the system. The first component $r_1$ acts as a constant force, while $r_2$ a constantly increasing displacement to the position. It can be tempting to think that a positive $r_2$ must add energy while a negative $r_2$ subtracts (or vice-versa for $r_1$ due to the sign difference), but this isn't true. The energy of an oscillator depends on the absolute value of the position and momentum, not the sign.'
• $h_{12}$ couple the position and momentum quadratures to each other, leading to oscillations between them. These notes and this stackexchange answer show how such a term can arise when considering a charged particle in a magnetic field in the Landau gauge.

To figure out what $h_{11}$ and $h_{22}$ mean, we can compare with a harmonic oscillator Hamiltonian: $\frac{\hat{p}^2}{2m}+\frac{1}{2}{m\omega^2}\hat{x}^2$.

• $h_{11}$ encodes a position-dependent force, which for a harmonic oscillator would be $+m\omega^2$. Thus if $h_{11}$ is positive there should be oscillations about the origin, while a negative $h_{11}$ would push the position further and further away from the origin.
• $h_{22}$ corresponds to $1/m$. Thus the smaller this is, the faster the particle should respond due to a reduced mass. A negative value would lead to exotic behaviour corresponding to a negative mass.

Note that none of these terms give damping. This is because damping results in a loss of energy, while Hamiltonians conserve energy. More formally, Hamiltonian evolution preserves volumes in phase space, while damping contracts volumes.

Two modes #

Let's add a second mode:

$\begin{aligned}\hat{\mathbf{r}}^T &= (\hat{x}_1,\hat{p}_1,\hat{x}_2,\hat{p}_2), \\ &=(\begin{array}{cc}\hat{\mathbf{r}}_1^T & \hat{\mathbf{r}}_2^T \end{array}),\end{aligned}$

where we have defined $\hat{\mathbf{r}}_i=(\hat{x}_i,\hat{p}_i)$. We will write the matrix $H$ in block-diagonal form

$\left(\begin{array}{cc}H_{11} & H_{12} \\ H_{21} & H_{22}\end{array}\right),$

where each $H_{ij}$ is a two by two matrix. The Hamiltonian is then

$\begin{aligned}\hat{H} &= \frac{1}{2}\hat{\mathbf{r}}^TH\mathbf{r}+\hat{\mathbf{r}}^T\vec{r}, \\ &=\frac{1}{2}\left(\hat{\mathbf{r}}^T_1H_{11}\hat{\mathbf{r}}_1+\hat{\mathbf{r}}^T_1H_{12}\hat{\mathbf{r}}_2+\hat{\mathbf{r}}^T_2H_{21}\hat{\mathbf{r}}_1+\hat{\mathbf{r}}^T_2H_{22}\hat{\mathbf{r}}_2\right)+\hat{\mathbf{r}}^T\vec{r}.\end{aligned}$

The terms $\hat{\mathbf{r}}_1^TH_{11}\hat{\mathbf{r}}_1$ and $\hat{\mathbf{r}}_2^TH_{22}\hat{\mathbf{r}}_2$ are our single-mode terms from before.

Let us explore the coupling term $\frac{1}{2}\left(\hat{\mathbf{r}}_1^TH_{12}\hat{\mathbf{r}}_2+\hat{\mathbf{r}}_2^TH_{21}\hat{\mathbf{r}}_1\right)$. The requirement $H^{\dagger}=H$ implies $H_{12}^{\dagger}=H_{21}$. We thus have for complex numbers $c_{jk}$:

$H_{12}=\left(\begin{array}{cc}c_{11} & c_{12} \\ c_{21} & c_{22}\end{array}\right);\; H_{21}=\left(\begin{array}{cc}c_{11}^* & c_{21}^* \\ c_{21}^* & c_{22^*}\end{array}\right).$

Using this to expand the coupling gives

$\mathrm{Re}(c_{11})\hat{x}_1\hat{x}_2+ \mathrm{Re}(c_{12})\hat{x}_1\hat{p}_2+ \mathrm{Re}(c_{21})\hat{p}_1\hat{x}_2+ \mathrm{Re}(c_{22})\hat{p}_1\hat{p}_2,$

so again we lose nothing in requiring $c_{ij}$ to be real. We can thus see that the $c_{ij}$ couple the position and momenta of different modes to each other.

Equations of motion #

We will now find equations of motion for $\hat{\mathbf{r}}$ directly in terms of $H$ and $\vec{r}$.

Single mode #

Again let's start with a single mode $\hat{\mathbf{r}}_1=(\hat{x}_1,\hat{p}_1)^T$. The Heisenberg equations are

$\hat{\mathbf{r}}_1'=i[\hat{H},\hat{\mathbf{r}}_1],$

where we interpret both the left and right hand sides as vectors with two components. The simplest term to compute is

$\begin{aligned}\left[\hat{\mathbf{r}}^T_1\vec{r},\hat{\mathbf{r}}_1\right] &= \left[\hat{x}_1r_1+\hat{p}_1r_2,\left(\begin{array}{c}\hat{x}_1 \\ \hat{p}_1\end{array}\right)\right], \\ &=\left(\begin{array}{c}-ir_2 \\ ir_1\end{array}\right), \\ &=-i\Omega_1\vec{r}.\end{aligned}$

The other term is

$\begin{aligned}&\left[\frac{1}{2}\hat{\mathbf{r}}_1^TH\hat{\mathbf{r}}_1,\hat{\mathbf{r}}_1\right] \\&\phantom{=}=\frac{1}{2}\left[h_{11}\hat{x}_1^2+h_{12}\hat{x}_1\hat{p}_1+h_{21}\hat{p}_1\hat{x}_1 +h_{22}\hat{p}_1^2,\left(\begin{array}{c}\hat{x}_1 \\ \hat{p}_1\end{array}\right)\right], \\ &\phantom{=}=\frac{i}{2}\left(\begin{array}{c}-(h_{12}+h_{21})\hat{x}_1-2h_{22}\hat{p}_1 \\ 2h_{11}\hat{x}_1+(h_{12}+h_{21})\hat{p}_1\end{array}\right).\end{aligned}$

To proceed further we must use the fact that $H$ is symmetric, so $h_{12}=h_{21}$. This gives us

$\begin{aligned}=i\left(\begin{array}{c}-h_{21}\hat{x}_1-h_{22}\hat{p}_1 \\ h_{11}\hat{x}_1 +h_{12}\hat{p}_1\end{array}\right)&=i\left(\begin{array}{cc}-h_{21} & -h_{22} \\ h_{11} & h_{12}\end{array}\right)\left(\begin{array}{c}\hat{x}_1 \\ \hat{p}_1\end{array}\right),\\ &=-i\Omega_1H\left(\begin{array}{c}\hat{x}_1 \\ \hat{p}_1\end{array}\right).\end{aligned}$

Multiplying by the factor of $i$ from the Heisenberg equations, we derive

$\hat{\mathbf{r}}_1'=\Omega_1H\hat{\mathbf{r}}_1+\Omega_1\vec{r}.$

We can use this to find a change of coordinates that removes the linear term. Let's define a new mode $\hat{\mathbf{q}}_1=(\hat{q}_1,\hat{p}_1)$ by

$\hat{\mathbf{q}}_1=\hat{\mathbf{r}}_1+H^{-1}\vec{r}.$

Then the new mode has equations of motion

$\begin{aligned}\hat{\mathbf{q}}_1' &=\Omega_1H\left(\hat{\mathbf{q}}_1-H^{-1}\vec{r}\right)+\Omega_1\vec{r}, \\ &=\Omega_1H\hat{\mathbf{q}}_1.\end{aligned}$

This transformation is only possible if $H$ is invertible, which we will briefly discuss at the end of this section.

More modes #

Now suppose we now have multiple modes:

$\hat{\mathbf{r}}^T=\left(\begin{array}{ccc}\hat{\mathbf{r}}_1^T & \cdots & \hat{\mathbf{r}}_n^T\end{array}\right),$

where $\hat{\mathbf{r}}_j=(\hat{x}_j,\hat{p}_j)^T$. We will similarly decompose

$\vec{r}=(\vec{r}_1,\cdots,\vec{r}_2),$

where $\vec{r}_j$ contains two elements, and $\vec{r}_1$ corresponds to $\vec{r}$ from the single-mode case. The matrix $H$ will be in block form

$H=\left(\begin{array}{ccc}H_{11} & \cdots & H_{1n} \\ \vdots & \ddots & \vdots \\ H_{n1} &\cdots & H_{nn}\end{array}\right),$

where $H_{jk}$ is a two-by-two matrix.

As before we must compute a few commutators.

$\left[\hat{\mathbf{r}}^T\vec{r},\hat{\mathbf{r}}\right]=\left[\hat{\mathbf{r}}_1^T\vec{r}_1+\cdots+\hat{\mathbf{r}}_n^T\vec{r}_n,\left(\begin{array}{c}\hat{\mathbf{r}}_1 \\ \vdots \\ \hat{\mathbf{r}}_n\end{array}\right)\right].$

We have $[\hat{\mathbf{r}}_j^T,\hat{\mathbf{r}}_k]=0$ if $j\ne k$, so this becomes

$\left(\begin{array}{c}-i\Omega_1\vec{r}_1 \\ \vdots \\ -i\Omega_1\vec{r}_n\end{array}\right)=-i\Omega\vec{r}.$

For the other term

$\left[\frac{1}{2}\hat{\mathbf{r}}^TH\hat{\mathbf{r}},\hat{\mathbf{r}}\right] = \frac{1}{2}\left[\sum_{jk}\hat{\mathbf{r}}^T_jH_{jk}\hat{\mathbf{r}}_k,\left(\begin{array}{c}\hat{\mathbf{r}}_1 \\ \vdots \\ \hat{\mathbf{r}}_n\end{array}\right)\right].$

Let's look at a single element:

$\frac{1}{2}\left[\sum_{jk}\hat{\mathbf{r}}^T_jH_{jk}\hat{\mathbf{r}}_k,\hat{\mathbf{r}}_l\right].$

There are several cases to consider. The easier two are:

• If $l$ is not equal to $j$ or $k$, then the commutator will be zero.

• If $j=k=l$ then we have the single-mode case:

  $\frac{1}{2}\left[\hat{\mathbf{r}}_l^TH_{ll}\hat{\mathbf{r}}_l\right]=-i\Omega_1H_{ll}\hat{\mathbf{r}}_l$.

Then there are two more complicated cases, when $l$ is equal to one of $j$ or $k$. We can handle these by noting that for an arbitrary vector $\vec{r}$, we have already shown

$\left[\hat{\mathbf{r}}^T\vec{a},\hat{\mathbf{r}}\right]=\left[\vec{a}^T\hat{\mathbf{r}},\hat{\mathbf{r}}\right]=-i\Omega_1\vec{a}.$

• If $l=k$, then we can take $\vec{a}= \sum_j\hat{\mathbf{r}}_j^TH_{jk}$.

• If $l=j$, then $\vec{a}= \sum_kH_{jk}\hat{\mathbf{r}}_k$.

Putting all this together we find

$\begin{aligned}&\frac{1}{2}\left[\sum_{jk}\hat{\mathbf{r}}^T_jH_{jk}\hat{\mathbf{r}}_k,\hat{\mathbf{r}}_l\right] \\ &\phantom{=}=-\frac{i}{2}\Omega_1H_{ll}\hat{\mathbf{r}}_l-\frac{i}{2}\sum_{k\ne l}\Omega_1 H_{lk}\hat{\mathbf{r}}_k-i\sum_{j\ne l}\Omega_1\hat{\mathbf{r}}^T_{j}H_{jl}.\end{aligned}$

Let's see how the last two terms simplify by writing them in terms of their components. For a symmetric matrix $A$ and a vector $\vec{b}$,

$(\vec{b}^TA)_k=\sum_jb_jA_{jk}=\sum_jA_{kj}b_j=(A\vec{b})_k,$

so $\hat{\mathbf{r}}_j^TH_{jl}=H_{jl}\hat{\mathbf{r}}_j$.

This gives us

$\begin{aligned}&\frac{1}{2}\left[\sum_{jk}\hat{\mathbf{r}}^T_jH_{jk}\hat{\mathbf{r}}_k,\hat{\mathbf{r}}_l\right] \\ &\phantom{=}=-i\Omega_1H_{ll}\hat{\mathbf{r}}_l-\frac{i}{2}\sum_{k}\Omega_1 H_{lk}\hat{\mathbf{r}}_k-\frac{i}{2}\sum_j\Omega_1 H_{jl}\hat{\mathbf{r}}_j.\end{aligned}$

The Hamiltonian matrices satisfy $H_{lk}=H_{kl}$, so the last two terms are equal. We thus have equations of motion

$\hat{\mathbf{r}}_l'=\Omega_1H_{ll}\hat{\mathbf{r}}_l+\sum_{j\ne l}\Omega_1 H_{lj}\hat{\mathbf{r}}_j+\Omega_1\vec{r}_l,$

which is equivalent to

$\hat{\mathbf{r}}' = \Omega H\hat{\mathbf{r}}+\Omega\vec{r}.$

As before if $H$ is invertible we can introduce a new variable $\hat{\mathbf{q}}=\hat{\mathbf{r}}+H^{-1}\vec{r}$, in terms of which we have simpler equations of motion:

$\hat{\mathbf{q}}'=\Omega H\hat{\mathbf{q}}.$

H invertible #

We have seen that when the matrix $H$ is invertible, the linear term proportional to $\vec{r}$ can be eliminated from the equations of motion by a change of variables. It is worth reflecting on what this means physically. Unfortunately I haven't had much of a chance to do this, but will try and revisit it in a week or so. Feel free to send me an email if you have any thoughts.

From the equations of motion, we can see that for the resulting modes to be linearly independent, $H$ must have full rank. Thus $H$ will only fail to be invertible if some of the modes do not evolve with time, or if the evolution spans a lower-dimensional subspace.

Conclusion #

We've seen that a general second-order bosonic Hamiltonian can be written as

$\hat{H}=\frac{1}{2}\hat{\mathbf{r}}^TH\hat{\mathbf{r}}+\hat{\mathbf{r}}^T\vec{r}.$

For $\hat{H}$ to be Hermitian, we required $\vec{r}$ real and $H$ Hermitian. However since the complex parts of $H$ don't affect the dynamics, we lose nothing in requiring it to be real and symmetric.

If we look at $H$ in two-by-two blocks, the diagonal elements $H_{jj}$ describe the oscillation frequencies and energies of the individual oscillators, and any coupling between their positions and momenta. Other elements $H_{jk}$ couple different modes to each other. The vector $\vec{r}$ applies a constant 'force', but can be eliminated with a change of variables provided $H$ is invertible.

With this parameterisation, we can work entirely in terms of finite-dimensional real matrices and vectors — quite an achievement given the infinite-dimensional Hilbert space. In a later post we might investigate how the eigenmodes can be found by diagonalising $H$. In this one we derived a matrix equation of motion:

$\hat{\mathbf{q}}'=\Omega H\hat{\mathbf{q}},$

The solution to which can be immediately written down:

$\hat{\mathbf{q}}(t)=\exp(\Omega Ht)\hat{\mathbf{q}}(0).$

This looks very similar to the Schrödinger evolution:

$|\psi(t)\rangle=\exp(-i\hat{H}t)|\psi(0)\rangle,$

with the matrix $\Omega$ playing a role analogous to the complex number $i$. Multiplying by $e^{i\theta}$ causes rotation in the two-dimensional complex plane. Similarly, the block-matrix form for $\Omega$ shows that it induces independent rotations among the individual two-dimensional $\hat{x}_j,\hat{p}_j$ subspaces. Thus in some sense, $\Omega$ is a 'high-dimensional $i$'. This has a deeper geometric meaning which I hope to investigate in future posts.

If you have any questions or comments, email me@ruvi.blog.

References #

[Ser17] Serafini, A. (2017). Quantum continuous variables: a primer of theoretical methods. CRC press.

---

• Next: Solving RC circuits with differential operators
• Previous: Geometric interpretation of continuous entropy