How to determine elements of the covariance matrix from homodyne measurements?

Question

Suppose we have a source which can repeatedly create an $n$-mode Gaussian state with covariance matrix $\sigma$. How can I use homodyne measurements to completely determine $\sigma$? If the quadratues are $\left( \hat q_1, \hat p_1, ..., \hat q_n, \hat p_n \right)$, then assume that the first moments $\langle \hat q_i \rangle$, $\langle \hat p_i \rangle$, where $i$ runs from $1$ to $n$, are zero.

From what I understand, homodyne measurement can be used to make measurements of the quadratures $\hat q_i$ and $\hat p_i$. Hence, we can determine terms like $\langle \hat q_1^2 \rangle$ by repeatedly measuring $\hat q_1$ etc. with homodyne measurement, squaring all the outcomes and taking the average of all of them. But what about cross terms, like $\langle \hat q_i \hat p_j + \hat q_j \hat p_i \rangle$?

Answer 1

The trick is that homodyne measurements allow you to measure not just $\hat{x}_i(\theta)=\hat{q}_i$ or $\hat{p}_i$ but also the value of any quadrature operator $\hat{q}_i\cos\theta-\hat{p}_i\sin\theta$. The value of $\theta$ is set by the phase of the local oscillator used for homodyne measurement (i.e., in homodyne detection one is really measuring expectation values of operators of the form $\hat{a}_i\hat{b}^\dagger+\hat{a}_i^\dagger\hat{b}$ with the state in mode $b$ being a coherent state $|\beta\rangle$ whose phase sets the value of $\theta$).

For example, taking $\theta=-\pi/4$, we can determine $$\langle \hat{x}_i(-\pi/4)^2\rangle\propto\langle \hat{q}_i^2\rangle+\langle \hat{p}_i^2\rangle+\langle \hat{q}_i\hat{p}_i+\hat{p}_i\hat{q}_i\rangle,$$ so knowledge of $\langle \hat{q}_i^2\rangle$, $\langle \hat{p}_i^2\rangle$, and $\langle \hat{x}_i(-\pi/4)^2\rangle$ allows you to determine $\langle \hat{q}_i\hat{p}_i+\hat{p}_i\hat{q}_i\rangle$.

As for cross terms between modes, one simply interferes some of those modes before performing homodyne detection on a single mode. For example, we can interfere modes 1 and 2 at a 50:50 beam splitter to yield the transformation $\hat{q}_1\to (\hat{q}_1+\hat{q}_2)/\sqrt{2}$ such that a homodyne measurement in mode 1 gives us access to $\langle \hat{q}_1^2+ 2\hat{q}_1 \hat{q}_2+\hat{q}_2^2\rangle$ ($\hat{q}_1$ and $\hat{q}_2$ commute), while interfering them on a symmetric beam splitter yields $\hat{q}_1\to (\hat{q}_1-\hat{p}_2)/\sqrt{2}$, so a homodyne measurement on mode 1 gives access to $\langle \hat{q}_1^2+ 2\hat{q}_1 \hat{p}_2+\hat{p}_2^2\rangle$ ($\hat{q}_1$ and $\hat{p}_2$ commute).

By varying the interferometric setup among the modes prior to homodyne detection and the phase of the local oscillator $\theta$ for the homodyne detection, one can build up all of the require terms of the covariance matrix. For non-Gaussian states, for which the covariance matrix does not contain all of the information about the state, more values of $\theta$ and interferometric setups will suffice to determine all of the information about the state; however, the amount of measurements required may be quite large or infinite for a generic state.

Answer 2

There is no fundamental difference between the diagonal terms and the cross terms, you just compute a different function of the measurement outcomes. However, you need more measurements.

Let's consider a two-mode state for simplicity. The covariance matrix is:

$$ \sigma = \begin{bmatrix} \Delta \hat{q}_a^2 & \Delta \hat{q}_a\hat{p}_a & \Delta \hat{q}_a\hat{q}_b & \Delta \hat{q}_a\hat{p}_b \\ \Delta \hat{p}_a\hat{q}_a & \Delta \hat{p}_a^2& \Delta \hat{p}_a\hat{q}_b & \Delta \hat{p}_a\hat{p}_b \\ \Delta \hat{q}_b\hat{q}_a & \Delta \hat{q}_b\hat{p}_a & \Delta \hat{q}_b^2 & \Delta \hat{q}_b\hat{p}_b \\ \Delta \hat{p}_b\hat{q}_a & \Delta \hat{p}_b\hat{p}_a & \Delta \hat{p}_b\hat{q}_b & \Delta \hat{p}_b^2 \end{bmatrix} $$

Here $\Delta \hat{q}\hat{p} = \frac{1}{2}\langle \{\hat{q},\hat{p}\} \rangle - \langle \hat{q} \rangle\langle \hat{p} \rangle$.

In a homodyne measurement you're continuously measuring the quadratures and get a time trace of the eigenvalues, for example: $q_a(t)$.

As you point out, the diagonal terms in the covariance matrix given by the variances of the quadratures: $\Delta q_a^2 = \langle \hat{q}_a^2 \rangle - \langle \hat{q}_a \rangle^2$ can be found by simply computing the variance of the recorded time trace $q_a(t)$.

The two-mode off-diagonal terms in the covariance matrix quantify the correlations between the modes $a$ and $b$, and these can also be computed from local measurements of the quadratures, but in this case two different parallel measurements on two modes, instead of a single measurement on one mode. For example:

$$ \Delta \hat{q}_a\hat{p}_b = \frac{1}{2}\langle \{\hat{q}_a,\hat{p}_b\} \rangle - \langle \hat{q}_a \rangle\langle \hat{p}_b \rangle = \frac{1}{2}\langle \hat{q}_a \hat{p}_b + \hat{p}_b \hat{q}_a \rangle -\langle \hat{q}_a \rangle\langle \hat{p}_b \rangle. $$

Because the modes $a$ and $b$ are orthogonal the quadrature operators commute, and the expression reduces to

$$ \Delta \hat{q}_a\hat{p}_b = \langle \hat{q}_a \hat{p}_b \rangle - \langle \hat{q}_a \rangle\langle \hat{p}_b \rangle. $$

To measure this term one does two parallel homodyne measurements, one of $\hat{q}$ on mode $a$ and one of $\hat{p}$ on mode $b$. The first term in the r.h.s. is then found by taking the mean value of the product of the time-traces $q_a(t)$ and $p_b(t)$ from the two parallel homodyne measurements on the two modes.

Finally, there are terms containing two different quadrature operators of the same mode, such as $\Delta \hat{q}_a \hat{p}_a$. These can be measured by switching to a different basis. First note that:

$$ \hat{q}_a\hat{p}_a + \hat{p}_a\hat{q}_a = i(\hat{a}^{\dagger}\hat{a}^{\dagger} - \hat{a}\hat{a}). $$

Then introduce the rotated quadratures: $$ \hat{x}_a = \frac{e^{i\frac{\pi}{4}}\hat{a}^{\dagger} + e^{-i\frac{\pi}{4}}\hat{a}}{\sqrt{2}}, \qquad \hat{y}_a = \frac{e^{-i\frac{\pi}{4}}\hat{a}^{\dagger} + e^{i\frac{\pi}{4}}\hat{a}}{\sqrt{2}}. $$

It's straightforward to verify that: $$ \hat{x}_a^2 - \hat{y}_a^2 = i(\hat{a}^{\dagger}\hat{a}^{\dagger} - \hat{a}\hat{a}) = \hat{q}_a\hat{p}_a + \hat{p}_a\hat{q}_a. $$

The diagonal terms in the blocks on the main diagonal can therefore be found as: $$ \Delta \hat{q}_a \hat{p}_a = \langle \hat{q}_a \hat{p}_a + \hat{p}_a \hat{q}_a \rangle - \langle \hat{q}_a \rangle \langle \hat{p}_a \rangle =\langle \hat{x}_a^2 \rangle - \langle \hat{y}_a^2 \rangle - \langle \hat{q}_a \rangle \langle \hat{p}_a \rangle, $$ which can be evaluated by taking the variances of the time traces $x_a(t)$ and $y_a(t)$, and the means of $q_a(t)$ and $p_a(t)$.