\section{Introduction}
Optical heterodyne detection is a commonly used method for measurement of amplitude and phase of spectral components of an optical beam \cite{kuri2003optical, delange1968optical, chtcherbakov2007optical, hall1981optical, jacobs1988optical}. Despite its frequent usage, complete and concise mathematical descriptions of the heterodyne measurement process are difficult to come by, with many sources simply stating a non-generalised final result \cite{keyes2013optical, menzies2005laser}. The lack of a readily available comprehensive description can make it difficult to deeply understand the method, including underlying assumptions, which increases the chance of experimental or analytical errors.

Here, we present a complete and concise derivation of optical heterodyne measurement of a beam consisting of regularly spaced frequency components. Such beams are commonly generated using electro-optic modulators, and ultrafast laser techniques, and where the number components is large, the optical field is referred to as a frequency comb. We attempt to use consistant notation, and where appropriate, give a physical interpretation to the mathematical results. 

At a given position any electric field strength can be written as
\begin{equation}
	E(t) = \frac{1}{2}\mathcal{E}(t)e^{-{\rm{i}}\omega_0t} + \frac{1}{2}\mathcal{E}^*(t)e^{{\rm{i}}\omega_0t}
	\label{a}
\end{equation}
where $\mathcal{E}(t)$ is the complex envelope function, $\mathcal{E}^*$ is its complex conjugate, and $\omega_0$ is the reference angular frequency. The choice of reference frequency is completely arbitrary, and can be positive or negative. However, the reference is \emph{usually} chosen so that the variation in time of $\mathcal{E}(t)$ is as small as possible. In practice, the reference frequency is often taken as the carrier frequency of a beam, so the complex envelope fully represents only the modulations to this carrier (both in phase, amplitude, and any additional frequency components).

The intensity of a given electric field is given by
\begin{equation}
	I(t) = c\varepsilon \frac{1}{T}\int_{-T/2}^{T/2} E^*(t)E(t){\rm{d}}t,
	\label{ba}
\end{equation}
where $T$ is an integer number of periods of the optical cycle, with minimum one period.

Subbing equation \ref{a} into this expression yields:
\begin{align}
	I(t) &= \frac{c\varepsilon}{4T}\int_{-T/2}^{T/2} \left[\mathcal{E}^*(t)\mathcal{E}(t) + \mathcal{E}(t)\mathcal{E}^*(t) + \mathcal{E}^2(t)e^{-{\rm{i}}2\omega_0t} + \mathcal{E}^{*2}(t)e^{{\rm{i}}2\omega_0t}  \right]{\rm{d}}t \label{b0}\\
	     & = \frac{c\varepsilon}{2T}\int_{-T/2}^{T/2} \mathcal{E}^*(t)\mathcal{E}(t){\rm{d}}t + \frac{c\varepsilon}{4T}\int_{-T/2}^{T/2} \left[ \mathcal{E}^2(t)e^{-{\rm{i}}2\omega_0t} + \mathcal{E}^{*2}(t)e^{{\rm{i}}2\omega_0t}  \right]{\rm{d}}t
	\label{b}
\end{align}
If we have chosen $\omega_0$ appropriately such that $\mathcal{E}$ is slowly varying in time, the expression in the first integral can be pulled out of the integral, and the integral evaluates to $T$. The two terms on the right oscillate at twice the optical frequency, and so will integrate to zero. Note this definition of intensity breaks down where the amplitudes vary significantly over the time scale of an optical cycle, such as with ultra-short laser pulses.

Averaging over the fast oscillations gives a convenient expression for intensity written in terms of the complex amplitudes:
\begin{align}
	I(t) &= \frac{c\varepsilon}{2}\mathcal{E}^*(t)\mathcal{E}(t)\label{c0}\\
	     &=\frac{c\varepsilon}{2}|\mathcal{E}(t)|^2.
	\label{c}
\end{align}

Now, when we interfere two beams, we will have the signal $E_{\rm{sig}}(t)$ and the local oscillator $E_{\rm{LO}}(t)$ fields.
The signal can be written as
\begin{equation}
	E_{\rm{sig}}(t) = \frac{1}{2}\mathcal{E}_{\rm{sig}}(t)e^{-{\rm{i}}\omega_0t} + \frac{1}{2}\mathcal{E}_{\rm{sig}}^*(t)e^{{\rm{i}}\omega_0t}.
	\label{d}
\end{equation}

A frequency comb with tooth spacing $\omega_{\rm{SB}}$ (`SB' is signifies `side band'), is written as a sum of fields each with a constant amplitude and a frequency that differs by an integer multiple of $\omega_{\rm{SB}}$. For compactness, a reference frequency can be pulled out of each term, and the total signal field can be written as equation \ref{d}, where the complex envelope is given by:
\begin{equation}
	\mathcal{E}_{\rm{sig}}(t) = \sum_{n}a_{n}e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})}
	\label{e}
\end{equation}
where $n$ can take any integer value (positive, negative or zero), $a_n$ is the real amplitude of the $n^{\rm{th}}$ comb tooth, and $\phi_{n}$ is the corresponding phase of that frequency component (cf. \ref{apxc1}). Note that this is simply a general expression that represents a frequency comb, it is not a recipe for how to generate one. It is well known that modulating either the phase or amplitude of a carrier periodically puts sidebands on that carrier, but here we are simply concerned with how to represent (and then measure) a field that \emph{does} have evenly spaced sidebands.

The full signal electric field is therefore written as:
\begin{equation}
	E_{\rm{sig}}(t) = \frac{1}{2}\sum_{n}a_{n}\left[e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})}e^{-{\rm{i}}\omega_0t} +e^{{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})}e^{{\rm{i}}\omega_0t}\right].
	\label{f}
\end{equation}

In a heterodyne measurement, the signal beam is combined with a local oscillator beam, which is shifted in angular frequency relative to the signal carrier by an amount $\Delta\omega$. The complex envelope can therefore be written as:
\begin{equation}
	\mathcal{E}_{\rm{LO}}(t) = be^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})},
	\label{g}
\end{equation}
where $b$ is the local oscillator real amplitude.

The full local oscillator electric field can then be written as:
\begin{equation}
	E_{\rm{LO}}(t) = \frac{b}{2}\left[e^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})}e^{-{\rm{i}}\omega_0t} + e^{{\rm{i}}(\Delta\omega t + \phi_{\rm LO})}e^{{\rm{i}}\omega_0t}\right].
	\label{h}
\end{equation}

Now, when the signal and local oscillator beams are combines on a beamsplitter, the output of a single port is just the sum of their fields. We will only be interested in a single port, so the phase shift introduced by reflection of a mirror is irrelevant to us and will be ignored (as will the halving of their amplitudes). The result is:
\begin{align}
	E_{\rm{tot}}(t) &= E_{\rm{sig}}(t) + E_{\rm{LO}}(t)\\
	     &= \frac{1}{2}\left[\sum_{n} a_{n}e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})} + be^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right]e^{-{\rm{i}}\omega_0t}
	        +\frac{1}{2}\left[\sum_{n} a_{n}e^{{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})} + be^{{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right]e^{{\rm{i}}\omega_0t}
	\label{i}
\end{align}
where the terms have been grouped by pulling out their common exponential factor.

We can see that equation \ref{i} is in the form of equation \ref{a}, where the complex envelope function is given by the term in the square brackets. We can therefore immediately write down the resulting time varying intensity for the combined field by substituting the envelope and its conjugate into equation \ref{c0}. Note that when calculation the intensity the summations over the signal teeth will multiply and hence the summations must be performed separately. The index label of one summation is therefore changed to $m$.

\begin{align}
	I_{\rm{tot}}(t) &= 	\frac{c\varepsilon}{2} \left[\sum_{n} a_{n}e^{-{\rm{i}}(n\omega_{\rm{SB}}t + \phi_{n})} + be^{-{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right] \times 
					\left[\sum_{m}a_{m}e^{{\rm{i}}(m\omega_{\rm{SB}}t + \phi_{m})} + be^{{\rm{i}}(\Delta\omega t + \phi_{\rm LO})} \right] \label{j0}\\
	     &=  \frac{c\varepsilon}{2}\left\{ 
	     b^2 + \sum_{n, m} a_{n}a_{m}e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]}
	      + \sum_{m}a_{m}be^{{\rm{i}}\left[(m\omega_{\rm{SB}}-\Delta\omega)t + \phi_{m} - \phi_{\rm LO}\right]}
	      + \sum_{n}a_{n}be^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	     \right\}\\
	     &=  \frac{c\varepsilon}{2}\left\{ 
	     b^2 + \sum_{n}a_{n}^2 
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]}
	      + \sum_{n}a_{n}be^{{\rm{i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	      + \sum_{n}a_{n}be^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	     \right\}
	      \label{j1}
\end{align}

So, equation \ref{j1} tells us exactly what the \emph{intensity} of the combined field will be at a given point in time. Clearly the intensity will vary periodically in time so this is commonly referred to as the \emph{beat note}, although in the case of a signal beam which is made up of many comb teeth, there will be multiple frequency components present.

We are often primarily interested in knowing something about these different frequency components, so a common way to analyse the heterodyne signal is using a spectrum analyser to display at the power spectrum of the photocurrent from a photodiode. The photocurrent is directly proportional to intensity, so the generated power spectrum is directly proportional to the power spectrum of the intensity. A spectrum analyser simply Fourier transforms the signal, and then displays the amplitude (or power, which is proportional to the square of the amplitude), so it is useful to see what the analytic form is of the Fourier transform of the time varying intensity, $\breve{I}(\omega)$:
\begin{equation}
	\breve{I}(\omega) = \frac{1}{\sqrt{2\pi}}\int I(t)e^{{\rm{-i}}\omega t}{\rm{d}}t,
	\label{k}
\end{equation}
where the integral is understood to be over $-\infty$ to $\infty$.

Substituting equation \ref{j1} into the Fourier transform yields:

\begin{align}
	\breve{I}_{\rm tot}(\omega) &=    \frac{c\varepsilon}{2\sqrt{2\pi}}  \int        
	    \Big\{
	     b^2 + \sum_{n}a_{n}^2 
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]} \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}be^{{\rm{i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	      + \sum_{n}a_{n}be^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}
	     \Big\}
	     e^{{\rm{-i}}\omega t}{\rm{d}}t  \label{l0}\\
	     &=  \frac{c\varepsilon}{2\sqrt{2\pi}}         
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)\int e^{{\rm{-i}}\omega t}{\rm{d}}t
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}\int e^{{\rm{i}}\left[(m-n)\omega_{\rm{SB}}t + \phi_{m} - \phi_{n}\right]}e^{{\rm{-i}}\omega t}{\rm{d}}t \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b\int e^{{\rm{i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}e^{{\rm{-i}}\omega t}{\rm{d}}t
	      + \sum_{n}a_{n}b\int e^{{\rm{-i}}\left[(n\omega_{\rm{SB}}-\Delta\omega)t + \phi_{n} - \phi_{\rm LO}\right]}e^{{\rm{-i}}\omega t}{\rm{d}}t
	     \Big\} 
	       \label{l2}\\
	     &=  \frac{c\varepsilon}{2\sqrt{2\pi}}         
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)\int e^{{\rm{-i}}\omega t}{\rm{d}}t
	     + \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}(\phi_{m} - \phi_{n})}\int e^{{\rm{i}}(m-n)\omega_{\rm{SB}}t}e^{{\rm{-i}}\omega t}{\rm{d}}t \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}be^{{\rm{i}} (\phi_{n} - \phi_{\rm LO})}\int e^{{\rm{i}}(n\omega_{\rm{SB}}-\Delta\omega)t }e^{{\rm{-i}}\omega t}{\rm{d}}t
	      + \sum_{n}a_{n}be^{{\rm{-i}}(\phi_{n} - \phi_{\rm LO})}\int e^{{\rm{-i}}(n\omega_{\rm{SB}}-\Delta\omega)t}e^{{\rm{-i}}\omega t}{\rm{d}}t
	     \Big\} 
	       \label{l3}	       
\end{align}


The Fourier integrals can be evaluated by noting that:
\begin{equation}
	\frac{1}{\sqrt{2\pi}}\int e^{{\rm i}\alpha t} e^{-{\rm i}\omega t} {\rm d}t= \sqrt{2\pi}\,\delta(\omega - \alpha),
	\label{m}
\end{equation}
where $\delta(\omega)$ is the Dirac delta function.


Substituting this in equation \ref{l2}, taking the corresponding value of $\alpha$ in each integrand gives:
\begin{align}
	\breve{I}_{\rm tot}(\omega) &=   \frac{c\varepsilon\sqrt{2\pi}}{2}         
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)\delta(\omega) \notag\\
	     &\qquad\qquad+ \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}(\phi_{m} - \phi_{n})}\delta(\omega-[m-n]\omega_{\rm{SB}}) \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{i}}(\phi_{n} - \phi_{\rm LO})}\delta(\omega-[n\omega_{\rm{SB}}-\Delta\omega]) \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{-i}}(\phi_{n} - \phi_{\rm LO})}\delta(\omega+[n\omega_{\rm{SB}}-\Delta\omega])
	     \Big\} 
	       \label{n}
\end{align}

Equation \ref{n} describes three combs of teeth the in frequency domain, each with spacing $\omega_{SB}$, and a peak also at zero frequency representing the DC level. The two combs that sum only over $n$ represent actual optical frequencies present in the signal beam. There are two of these combs because each actual unique signal comb tooth has a positive and negative frequency component when the original electric field is represented in a form such as \ref{apxd}. 

The term that sums over $n$ and $m$ represents the beatnotes formed when the signal beam comb teeth beat against \emph{each other}. Different pairs of frequency components with equal frequency spacing will beat at the same frequency, and may add or subtract to that beatnote depending on their phase. For a set of sidebands produced by a phase only electro-optic modulator, the different phases between all frequency components add in a way so this ``cross comb" does not show up at all.

Depending on the number of teeth, the sideband spacing, and the local oscillator frequency offset, it is possible for teeth belonging to all three combs to overlap. It is generally undesirable for the teeth to overlap, since this prevents measurement of the amplitude and phase of individual of the individual frequency components that make up the signal beam frequency comb. Figures \ref{fig1}, \ref{fig2}, \ref{fig3} and \ref{fig4} show how combs can begin to intersect as the number of teeth get larger. Note that $\Delta\omega$ has deliberately been chosen to be an integer plus one third times the sideband spacing, so that even though all three combs become interspersed, the individual comb teeth do not lie on top of one another. The amplitude of the signal comb teeth was was set to be a Gaussian simply as a visual aid.
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_1_teeth.pdf}
    \caption{Beatnote from comb with single tooth.}
    \label{fig1}
\end{figure}
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_3_teeth.pdf}
    \caption{Beatnote from comb with 3 teeth.}
    \label{fig2}
\end{figure}
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_10_teeth.pdf}
    \caption{Beatnote from comb with 10 teeth.}
    \label{fig3}
\end{figure}
\begin{figure}
    \includegraphics[width=\columnwidth]{figures/non_integer_frequency_offset/freq_comb_beatnote_30_teeth.pdf}
    \caption{Beatnote from comb with 30 teeth.}
    \label{fig4}
\end{figure}

Equation \ref{n} gives the Fourier amplitude at a given frequency, but in reality, any measurement we make of the Fourier amplitudes is over some small range of frequencies. In a spectrum analyser, this range is set by the resolution bandwidth RBW, and in a numerical discrete Fourier transform of a time domain signal, it is set by the spacing of the frequency array (which is the inverse of the total length of time sampled). So experimentally, to extract the actual value for the Fourier amplitude at some frequency $\bar{\omega}$, we must integrate over some small frequency range : $\omega=\bar{\omega}\pm\epsilon$, and the resulting quantity will be denoted $\bar{\breve{I}}_{\rm tot}(\bar{\omega})$. Integration over the Dirac deltas means that the resulting values of $\bar{\breve{I}}_{\rm tot}(\bar{\omega})$ will be finite.

\begin{align}
	\bar{\breve{I}}_{\rm tot}(\bar{\omega}) &= \int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}  \breve{I}_{\rm tot}(\omega)  {\rm d}\omega  \label{o1} \\
	    &=\frac{c\varepsilon\sqrt{2\pi}}{2}         
	    \Big\{
	     \left(b^2 + \sum_{n}a_{n}^2\right)
	     \int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega){\rm d}\omega \notag\\
	     &\qquad\qquad+ \sum_{\substack{n,m\\ n\neq m}} a_{n}a_{m}e^{{\rm{i}}(\phi_{m} - \phi_{n})}\int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega-[m-n]\omega_{\rm{SB}}) {\rm d}\omega \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{i}}(\phi_{n} - \phi_{\rm LO})}\int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega-[n\omega_{\rm{SB}}-\Delta\omega]) {\rm d}\omega \notag\\
	      &\qquad\qquad+ \sum_{n}a_{n}b e^{{\rm{-i}}(\phi_{n} - \phi_{\rm LO})}\int_{\bar{\omega}-\epsilon}^{\bar{\omega}+\epsilon}\delta(\omega+[n\omega_{\rm{SB}}-\Delta\omega]) {\rm d}\omega
	     \Big\} 
	       \label{o2}
\end{align}

Equation \ref{o2} provides us with a simple method to determine the amplitude and phase of the optical signal comb teeth by measuring the radio frequency beat note. If the time domain radio signal is measured on an oscilloscope, it can be Fourier transformed numerically to show the radio-frequency combs. Assuming the teeth of the different combs are not overlapped then the Fourier amplitude of at the frequency of the $n$th tooth has a contribution from only a single term of equation \ref{o2}. Under these conditions, the optical amplitude and phase of the $n$th tooth is related to the Fourier transform of the detected intensity simply by:
\begin{equation}
	\bar{\breve{I}}_{\rm tot}(\bar{\omega}_n) = a_{n}b e^{{\rm{i}}(\phi_{n} - \phi_{\rm LO})},
	\label{p}
\end{equation}
so that the optical amplitude of the tooth is:
\begin{equation}
	a_{n} = \frac{|\bar{\breve{I}}_{\rm tot}(\bar{\omega}_n)|}{b},
	\label{q}
\end{equation}
and the optical phase of the tooth is:
\begin{equation}
	\phi_{n} = {\rm Arg}[\bar{\breve{I}}_{\rm tot}(\bar{\omega}_n)]  + \phi_{\rm LO}.
	\label{r}
\end{equation}

It is worth noting that the optical phase of the local oscillator $\phi_{\rm LO}$ will not usually be fixed in a heterodyne measurement. While it is certainly possible to set up some phase locking circuit by either feeding back to the laser or to a piezo controlled mirror, generally the $\phi_{\rm LO}$ will drift over time, as is seen in any unlocked interferometer. Therefore, the recovered value of $\phi_{n}$ will also drift over time. However, so long as $\phi_{\rm LO}$ stays constant over the time scale of a single measurement, then the \emph{difference} between all the $\phi_{n}$s will always be measured to be the same, even though $\phi_{\rm LO}$ may drift.

A final word about the beat note formed by the comb teeth beating against each: just because there is a comb, \emph{does not} mean you will actually observe any modulation of the intensity. Looking at the second term in equation \ref{n} (the term responsible for teeth beating with each other), it can be seen that pairs of $n,m$ with the same difference, do not necessarily have the same phase. If the phase of different pairs is not the same, then it is possible (indeed likely) that subsequent pairs will not all add in a way that increases the magnitude, and may totally cancel each other out. This is the case when a beam is phase modulated using an EOM. If the beam intensity is subsequently measured with only a photodiode, no intensity modulation is seen at all. However beating with a frequency shifted local oscillator reveals that there may be many other frequency components present.


