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

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$$?

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.

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)$$.