Introduction #

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


One such example is the harmonic oscillator: H^=p^22m+12mω2x^2\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 x^,p^\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 x^,p^\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 x^,p^\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


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 x^,p^\hat{x},\hat{p}^{\dagger}, for example:


Note that for H^\hat{H} to be Hermitian, aa must be real. Looking at the Heisenberg equation of motion for p^\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:


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:


While the Hamiltonian contains quadratic terms, the resulting Heisenberg equations of motion for x^,p^\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

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

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

[r^,r^T]=([x^,x^][x^,p^][p^,x^][p^,p^])=(0ii0).[\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

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

in terms of which the commutator becomes


In general we will have multiple modes r^=(x^1,p^1,x^2,p^2,,x^n,p^n)T\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 [x^j,p^k]=iδjk[\hat{x}_j,\hat{p}_k]=i\delta_{jk}, the commutator will now be




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

j=1nAj=(A1000A2000).\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

j=13Ω1=Ω1Ω1Ω1,=(010000100000000100001000000001000010).\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 r^\hat{\mathbf{r}}, a general second-order Hamiltonian H^\hat{H} can be decomposed as


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

We can see from the above equation that for H^\hat{H} to be Hermitian, HH must be Hermitian and r\vec{r} must be real. We will see in the next section however that any imaginary parts of HH cancel out and don't affect the Hamiltonian, so we may as well take HH 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 HH and r\vec{r} mean.

A single mode #

The simplest case is just a single mode. Letting HH have elements hjkh_{jk} and r\vec{r} elements rjr_j, the Hamiltonian expands as


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

h12x^p^+h12p^x^=(h12+h12)p^x^+ih12,=2Re(h12)p^x^+ih12,\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 Re\mathrm{Re} denotes the real part. Any complex components of h12h_{12} thus only show up as a constant, which doesn't affect the dynamics of the system. For this reason we can just take HH to be a real symmetric matrix without any loss of generality.

The Heisenberg equations of motion are then

x^=i[H^,x^]=h12x^+h22p^+r2,p^=i[H^,p^]=h11x^h12p^r1.\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}

  • r\vec{r} gives a constant driving to the system. The first component r1r_1 acts as a constant force, while r2r_2 a constantly increasing displacement to the position. It can be tempting to think that a positive r2r_2 must add energy while a negative r2r_2 subtracts (or vice-versa for r1r_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.'
  • h12h_{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 h11h_{11} and h22h_{22} mean, we can compare with a harmonic oscillator Hamiltonian: p^22m+12mω2x^2\frac{\hat{p}^2}{2m}+\frac{1}{2}{m\omega^2}\hat{x}^2.

  • h11h_{11} encodes a position-dependent force, which for a harmonic oscillator would be +mω2+m\omega^2. Thus if h11h_{11} is positive there should be oscillations about the origin, while a negative h11h_{11} would push the position further and further away from the origin.
  • h22h_{22} corresponds to 1/m1/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:

r^T=(x^1,p^1,x^2,p^2),=(r^1Tr^2T),\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 r^i=(x^i,p^i)\hat{\mathbf{r}}_i=(\hat{x}_i,\hat{p}_i). We will write the matrix HH in block-diagonal form

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

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

H^=12r^THr+r^Tr,=12(r^1TH11r^1+r^1TH12r^2+r^2TH21r^1+r^2TH22r^2)+r^Tr.\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 r^1TH11r^1\hat{\mathbf{r}}_1^TH_{11}\hat{\mathbf{r}}_1 and r^2TH22r^2\hat{\mathbf{r}}_2^TH_{22}\hat{\mathbf{r}}_2 are our single-mode terms from before.

Let us explore the coupling term 12(r^1TH12r^2+r^2TH21r^1)\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=HH^{\dagger}=H implies H12=H21H_{12}^{\dagger}=H_{21}. We thus have for complex numbers cjkc_{jk}:

H12=(c11c12c21c22);  H21=(c11c21c21c22).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

Re(c11)x^1x^2+Re(c12)x^1p^2+Re(c21)p^1x^2+Re(c22)p^1p^2,\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 cijc_{ij} to be real. We can thus see that the cijc_{ij} couple the position and momenta of different modes to each other.

Equations of motion #

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

Single mode #

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


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

[r^1Tr,r^1]=[x^1r1+p^1r2,(x^1p^1)],=(ir2ir1),=iΩ1r.\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

[12r^1THr^1,r^1]==12[h11x^12+h12x^1p^1+h21p^1x^1+h22p^12,(x^1p^1)],==i2((h12+h21)x^12h22p^12h11x^1+(h12+h21)p^1).\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 HH is symmetric, so h12=h21h_{12}=h_{21}. This gives us

=i(h21x^1h22p^1h11x^1+h12p^1)=i(h21h22h11h12)(x^1p^1),=iΩ1H(x^1p^1).\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 ii from the Heisenberg equations, we derive


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


Then the new mode has equations of motion

q^1=Ω1H(q^1H1r)+Ω1r,=Ω1Hq^1.\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 HH is invertible, which we will briefly discuss at the end of this section.

More modes #

Now suppose we now have multiple modes:

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

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


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

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

where HjkH_{jk} is a two-by-two matrix.

As before we must compute a few commutators.

[r^Tr,r^]=[r^1Tr1++r^nTrn,(r^1r^n)].\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 [r^jT,r^k]=0[\hat{\mathbf{r}}_j^T,\hat{\mathbf{r}}_k]=0 if jkj\ne k, so this becomes

(iΩ1r1iΩ1rn)=iΩr.\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

[12r^THr^,r^]=12[jkr^jTHjkr^k,(r^1r^n)].\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:


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

  • If ll is not equal to jj or kk, then the commutator will be zero.

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


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


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

  • If l=jl=j, then a=kHjkr^k\vec{a}= \sum_kH_{jk}\hat{\mathbf{r}}_k.

Putting all this together we find

12[jkr^jTHjkr^k,r^l]==i2Ω1Hllr^li2klΩ1Hlkr^kijlΩ1r^jTHjl.\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 AA and a vector b\vec{b},


so r^jTHjl=Hjlr^j\hat{\mathbf{r}}_j^TH_{jl}=H_{jl}\hat{\mathbf{r}}_j.

This gives us

12[jkr^jTHjkr^k,r^l]==iΩ1Hllr^li2kΩ1Hlkr^ki2jΩ1Hjlr^j.\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 Hlk=HklH_{lk}=H_{kl}, so the last two terms are equal. We thus have equations of motion

r^l=Ω1Hllr^l+jlΩ1Hljr^j+Ω1rl,\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

r^=ΩHr^+Ωr.\hat{\mathbf{r}}' = \Omega H\hat{\mathbf{r}}+\Omega\vec{r}.

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

q^=ΩHq^.\hat{\mathbf{q}}'=\Omega H\hat{\mathbf{q}}.

H invertible #

We have seen that when the matrix HH is invertible, the linear term proportional to r\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, HH must have full rank. Thus HH 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


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

If we look at HH in two-by-two blocks, the diagonal elements HjjH_{jj} describe the oscillation frequencies and energies of the individual oscillators, and any coupling between their positions and momenta. Other elements HjkH_{jk} couple different modes to each other. The vector r\vec{r} applies a constant 'force', but can be eliminated with a change of variables provided HH 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 HH. In this one we derived a matrix equation of motion:

q^=ΩHq^,\hat{\mathbf{q}}'=\Omega H\hat{\mathbf{q}},

The solution to which can be immediately written down:

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

This looks very similar to the Schrödinger evolution:


with the matrix Ω\Omega playing a role analogous to the complex number ii. Multiplying by eiθ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 x^j,p^j\hat{x}_j,\hat{p}_j subspaces. Thus in some sense, Ω\Omega is a 'high-dimensional ii'. This has a deeper geometric meaning which I hope to investigate in future posts.

If you have any questions or comments, email

References #

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