\section{Introduction}

In high-precision measurements, our understanding of physics is predominantly limited by quantum noise, arising due to the fundamental quantum fluctuations of the probing fields~\cite{Braginsky, Caves1980, Gardiner1991, Clerk2010a}. This is particularly true for laser interferometric gravitational-wave detectors~\cite{Adhikari2014} where the quantum shot noise dominates at high frequencies due to the positive dispersion of the arm cavities~\cite{Miao}. Quantum and classical noises are also limiting factors in quantum optomechanical  experiments~\cite{ Chen2013, Aspelmeyer2014} and searches for new physics using an interferometer~\cite{Sikivie1983, Derocco2018}. To achieve a maximal signal-to-noise ratio, it is essential to engineer the frequency-dependent response of the measurement devices depending on the frequency content of the signal being measured. For example, advanced gravitational wave detectors are tuned to have maximum sensitivity in a frequency range containing the binary black hole inspiral waveform, however not all of the binary neutron star inspiral waveform is observed~\cite{TheLIGOScientificCollaboration2019}. Quantum filters are designed to engineer this response. As illustrated in Fig.~\ref{fig:quantum-filters}, there are three ways that the measurement device can be augmented with quantum filters. First, the input filter, coupling the noise input to the probe degrees of freedom, shapes how the quantum fluctuations enter the device. Next, the coherent feedback filter, coupled to the probe degrees of freedom and input-output fields, modifies the dynamics of the probe~\cite{James2008,Mabuchi2008a,Hamerly2012,Jacobs2014}. This can enhance the response to the signal of interest when the quantum system is converted into a probe coupled to a classical signal. Finally, the output filter, coupling the probe degrees of freedom to the readout port, modifies the response of readout to the detector's the output field. As a simple example, an optical Fabry-Perot cavity is used as an 
input filter to in implementing frequency-dependent squeezed light~\cite{Kimble2000,Oelker2016c,Schnabel2017}.

Until now, formulating a physical realisation of a given quantum filter with a desired frequency response required a combination of intuition and prior experience, making more complicated frequency responses difficult to engineer. We adopt the general formalism for describing linear stochastic quantum networks and the synthesis of such networks, recently developed by the quantum control community~\cite{Gough2007, James2008, Gough2009, Gough2009a, Tezak2012a, Combes2016, Nurdin2009, Nurdin2010a,Nurdin2010b,Nurdin2016, Grivopoulos2016,Nurdin2017,Grivopoulos2017,Petersen2018}. This allows us in this paper to develop a formalism for systematically realising quantum filters for high-precision measurements directly from their frequency-domain transfer matrices. Although all of the components of the formalism existed previously this is the first time the entire process of going from the transfer matrix to physical realisation has been written down, and the concept of using the formalism to produce quantum filters for high-precision measurement is a totally original contribution, as well as the method of transforming to an unrealisable state-space to a realisable one.

\begin{figure}[!b]
    \centering
    \includegraphics[width=\linewidth]{graphics/quantum-filters-flowchart.pdf}
    \caption{Flowchart illustrating the different places quantum filters can be used within a quantum measurement device. We consider a generic device consisting of a probe (e.g.\, a mirror-endowed test mass or an atomic ensemble) coupled to some classical signal, which receives an input (e.g.~non-classical squeezed light) and whose output field is measured by the readout scheme (e.g.homodyne readout).}
    \label{fig:quantum-filters}
\end{figure}

The starting point of this approach is to convert the frequency-domain transfer matrix into a
state-space representation. However, such a mapping is not unique, and in fact infinitely many different state-space representations exist.  
We apply the concept of physical realisability, which was first introduced in~\cite{James2008} and further discussed in~\cite{Gough2009}\cite[Chapter 2]{Nurdin2017}. It tells us whether a given time-domain state-space representation of the system obeys quantum mechanics.
Requiring that the state-space is physically realisable sets a significant constraint on the range of this mapping. However, when a state-space realisation satisfies certain conditions, following \cite{Shaiju2012}, we can transform that state-space representation of the frequency-domain transfer matrix to one that is physically realisable. 
Applying the network synthesis theory can then lead to the physical setup of the quantum filter. 

This approach has powerful implications on how both passive and active quantum filters are designed, making the realisation of filters with arbitrarily complicated frequency responses a possibility. Since in principle we can view the entire measurement device as a many degrees-of-freedom quantum filter, this approach also provides a new paradigm for designing optimal quantum measurement devices. The outline of 
this paper goes as follows. In 
Section\,\ref{sec:approach}, we
present the mathematical details of this 
approach. In Section\,\ref{sec:unstable-filter}, 
we apply it to find a new optical 
realisation of the 
active unstable filter, the original proposal
of which was based upon optomechanics. 
We show that the optical loss is the limiting factor for the new 
realisation. In Section\,\ref{sec:discussion}, we summarise our result and provide 
an outlook of applying this approach to the
design of the optimal measurement devices. 

\section{Direct approach}
\label{sec:approach}

We now provide the details of the approach. The process to find a physical realisation, e.g.~an optical layout and its associated parameters, from a given set of transfer functions is general to multi-input multi-output lossless linear quantum systems; losses and other noise sources can be added later by augmenting the system description. Our starting point is the frequency-domain transfer function matrix, which is the square matrix that relates the frequency-domain system outputs $\mathbf{y}(s)$ to its inputs $\mathbf{u}(s)$:
\begin{align}
    G(s) &= C(-sI-A)^{-1}B + D,\label{eq:tf}\\ \mathbf{y}_i(s) &= \sum_j G_{ij}(s) \mathbf{u}_j(s),\nonumber
\end{align}
where $(A, B, C, D)$ are the system matrices as defined below, and $I$ is the identity matrix, and we assume that the number of inputs is equal to the number of outputs. Here the Laplace transform is defined as $f(s) = \int_{0^-}^{\infty} e^{+st} f(t) \mathrm{d}t$, with the lower bound at $t = 0^-$ so that an impulse can be added at $t = 0$. For a given transfer matrix a non-unique state-space representation can be found of the form~\cite{Luenberger1967, Ackermann1971, Kailath1980, Antoniou1988, Note5}:
\begin{align}
	\dot{\mathbf{x}} &= A\,\mathbf{x} + B\,\mathbf{u}, \label{eq:ss-1}\\
	\mathbf{y} &= C\,\mathbf{x} + D\,\mathbf{u}, \label{eq:ss-2}
\end{align}
which is a \emph{non-unique} time-domain representation of the system's dynamics. Here the quantity $\mathbf{x} \in \mathbb{L}^{2 n\times 1}$ ($\mathbb{L}$ being the space of linear operators on the relevant Hilbert space $\mathcal{H}$) is a vector of conjugate operator pairs representing the internal $n$ degrees of freedom of the system, $\mathbf{u} \in \mathbb{L}^{2m\times 1}$ is the vector of $m$ system inputs, and $\mathbf{y} \in \mathbb{L}^{2m\times 1}$ is the vector of $m$ system outputs. Note that in a quantum mechanical state-space, two conjugate operators are used to represent each individual degree of freedom of the system hence the factors of $2$. In the context of quantum optomechanics, $\hat{x}$ represents the cavity/oscillator eigenmodes for the cavities and mechanical oscillators in the system, while $\hat{u}$ and $\hat{y}$ are continuous Bosonic fields in free space~\cite{Blow1990,Nurdin2009}.
The dynamical matrix $A \in \mathbb{C}^{2n\times2n}$ describes the internal dynamics of the system, the input matrix $B \in \mathbb{C}^{2n\times2m}$ describes the coupling of the input into the system, the output matrix $C \in \mathbb{C}^{2m\times2n}$ describes the coupling of the system to the output, and the direct feedthrough matrix $D \in \mathbb{C}^{2m\times2m}$ describes the coupling of the input directly to the output. $(A, B, C, D)$ are together called the \emph{system matrices} and fully describe the linear dynamics of the system.

The system is called \emph{physically realisable} (and a corresponding physical realisation can be designed) if, in the Heisenberg picture evolution of the system, the commutation relations are preserved~\cite{James2008}:
\begin{equation}
	\forall i, j\ \ \ \mathrm{d}[\mathbf{x}_i,\mathbf{x}_j] = 0,\ \ [\mathbf{y}_i(t), \mathbf{y}^\dagger_j(t')] = \delta(t-t')\delta_{ij},
	\label{eq:physical-realisability-0}
\end{equation}
where the differential is treated using the quantum It\^o rule, meaning that the cross-products of the differentials of the operators must be calculated~\cite{Hudson1984,Parthasarathy1992,Bouten2006}.
The conditions on the system matrices for all such evolutions to preserve these commutation relations are found by using Eqs.\,\eqref{eq:ss-1} and\,\eqref{eq:ss-2} to calculate the increment of the system state $\mathrm{d}\mathbf{x}_i$ in Eq.\,\eqref{eq:physical-realisability-0} for an infinitesimal time period $\mathrm{d}t$. For an $n$ degree-of-freedom system described using complex mode operators (as are usually used in quantum optics) such that $\mathbf{x} = (\hat{a}_1, \hat{a}_1^\dagger;\dots;\hat{a}_n,\hat{a}_n^\dagger)^T$ with $m$ inputs and outputs described by $\mathbf{u} = (\hat{u}_1, \hat{u}_1^\dagger;\dots;\hat{u}_m,\hat{u}_m^\dagger)^T$ and similarly for $\mathbf{y}$, the constraints on the system matrices are given by
\begin{align}
	A J + J A^\dagger + B J_m B^\dagger &= 0, 
	\label{eq:physical-realizability-1}\\
	J C^\dagger + B J_m D^\dagger &= 0, \label{eq:physical-realizability-2}\\
	\ \ D J_m D^\dagger &= J_m,
	\label{eq:physical-realizability-3}
\end{align}
where $J = \text{diag}(1, -1; \dots; 1, -1)\in \mathbb{R}^{2n \times 2n}$ and $J_m = \text{diag}(1, -1; \dots; 1, -1)\in \mathbb{R}^{2m \times 2m}$\footnote{As discussed in Section.~\ref{appendix:r-matrix} of this paper, this matrix takes a different form when using Hermitian observable quadrature operators.}. See Appendix A of~\cite{Gough2009} for a proof of these constraints. So now we have a restriction on the possible system matrices that can lead to a physically realisable system.

Now we consider how to generate such a physically-realisable state-space model from the system's transfer matrix. The conventional procedure for transforming the transfer matrix to a minimal state-space model is outlined in Refs.\,~\cite{Luenberger1967, Kailath1980}. Such a state-space model constructed from a pole-zero form transfer matrix is minimal if the number of internal degrees of freedom (e.g.~the number of pairs of conjugate ladder operators describing the system state) is equal to the highest polynomial order in the frequency $s$ among all of the transfer functions in the transfer matrix. Generally such a procedure will lead to system matrices $(A',B',C',D')$ that do not satisfy Eqs.~\eqref{eq:physical-realizability-1} and~\eqref{eq:physical-realizability-2} and therefore cannot be physically realised. 

Here we show a method allowing us to transform a minimal realisation $(A',B',C',D')$ to a physically realisable counterpart $(A, B, C, D)$, given that the transfer matrix ${G}(s)$ obeys a condition that will be given in Eq.~\eqref{eq:j-j-unitary} and the state-space realisation satisfies the conditions given in \cite[Theorem 3]{Shaiju2012} (to be recalled below). The transformation is achieved by looking for a Hermitian matrix $X$ that obeys the constraints:
\begin{align}
	A' X + X (A')^\dagger + B' J_m (B')^\dagger &= 0,
	\label{eq:physical-realizability-X1} \\
	X (C')^\dagger + B' J_m (D')^\dagger &= 0.
	\label{eq:physical-realizability-X}
\end{align}
This matrix $X$ can be written in the form of a similarity transformation $X = TJT^\dagger$ for some non-singular matrix $T$. Substituting this into Eqs.~\eqref{eq:physical-realizability-1} and~\eqref{eq:physical-realizability-2} we see that the conditions are satisfied after making the transformations,
\begin{equation}
A = T^{-1} A' T,\ B = T^{-1} B',\ C = C' T,\ D = D',
\label{eq:applying-similarity-transformation}
\end{equation}
to find the physically realisable state-space $(A, B, C, D)$. The existence of $X$, i.e. $T$, is guaranteed by the 
symplectic condition imposed on any physically realisable transfer matrix ${G}(s)$ and direct-feed matrix $D$\,\cite{Shaiju2012}:
\begin{equation}
    {G}^\dag(s^*)J_m{G}(-s) = J_m, 
    \label{eq:j-j-unitary}
\end{equation}
and fulfilment of the conditions of \cite[Theorem 3]{Shaiju2012}, which says that $\lambda_i(A)+\lambda_j(A)^* \neq 0$ for any pair of eigenvalues $\lambda_i(A)$ and $\lambda_j(A)$ of $A$, and the feedthrough matrix $D$ is unitary and satisfies \eqref{eq:physical-realizability-3}.

Now that we have shown how to obtain the physically realisable state-space model $(A, B, C, D)$ from the transfer matrix obeying Eq.~\eqref{eq:j-j-unitary} we can infer the physical realisation. We describe the realisation chiefly using the generalised open oscillator~\cite{Nurdin2009} formalism. This is a general formalism describing open quantum systems with arbitrary internal linear dynamics and input-output couplings, providing a language for describing and analysing systems with internal degrees of freedom coupled to external continuum fields, such as quantum measurement devices and quantum filters. As shown in~\cite{James2008} for an $n$ degree-of-freedom system, when the direct-feed matrix $D$ is symplectic and unitary (i.e.~it satisfies \eqref{eq:physical-realizability-3} and $D^{\dag}D = DD^{\dag}=I$), there is a one-to-one correspondence between the system matrices $(A, B, C, D)$ and the generalised open oscillator which is parameterized by a triplet $(S, \hat{L}, \hat{H})$~\cite{Nurdin2009, Tezak2012a, Combes2016}. Here, the scattering matrix $S \in \mathbb{C}^{m \times m}$ describes the transformation of the input fields through a passive network, i.e.~any passive pre-processing of the system's input fields. The coupling operator $\hat{L} = K \mathbf{x}$ where $K \in \mathbb{C}^{m\times 2n}$ describes the coupling between the input and output fields and the internal degrees of freedom, e.g.~equivalent to the usual input-output Langevin equations~\cite{Chen2013} when the input-output fields are coupled to the internal fields by a mirror. The Hamiltonian $\hat{H}$ describes the free evolution of the internal system dynamics as if the system were closed. The relation between the system matrices and the generalised open oscillator parameters is given by,
\begin{align}
\label{eq:SKR}
\begin{split}
	S_{kl} &= D_{2k-1,2l-1}\,, \quad 
	\hat{L} = [\begin{array}{cc} I & 0 \end{array}]P\, C\, \mathbf{x},
	\\ \hat{H} &= \frac{i}{4} \hbar\,\mathbf{x}^\dag  \left(J A - A^\dagger J\right) \mathbf{x},
\end{split}
\end{align}
where $I$ is the identity matrix and $0$ is the null matrix, both of dimension $m$, $P$ is the permutation matrix that maps $\mathbf{u} = (\hat{u}_1,\hat{u}_1^{\dag};\ldots;\hat{u}_m,\hat{u}_m^{\dag})^T$ to $(\hat{u}_1,\hat{u}_2,\ldots,\hat{u}_m,\hat{u}_1^{\dag},\hat{u}_2^{\dag},\ldots,\hat{u}_m^{\dag})^T$, and $\hat{H}$ is derived in Section.~\ref{appendix:r-matrix}. The total Hamiltonian is then
\begin{equation}
    \hat{H}_\text{tot} = \hat{H} + i\hbar[\begin{array}{cc} \hat{L}^{\dag} & -\hat{L}^{T} \end{array}] \mathbf{u},
\end{equation}
where the input fields $\mathbf{u}$ are pre-processed by a static passive network described by $S$. Now we have achieved the full physical Hamiltonian describing the system starting from the transfer matrix describing the frequency-domain input-output behaviour.



In the case where only one of the internal degrees of freedom is coupled to the input-output fields, it could be straightforward to construct the physical realisation by inspection, as in the  illustrative example of the unstable filter
discussed in the next section. However, having one internal degree of freedom is not always the case. Systems consisting of more than one internal degree-of-freedom can first be sub-divided into separate one degree-of-freedom systems coupled via direct interaction Hamiltonians via the main synthesis theorem proved in~\cite{Nurdin2009}. These systems can then be systematically realised by connecting the individual one degree-of-freedom systems in series, and overlapping them accordingly, giving a systematic way to construct the physical realisation regardless of complexity~\footnote{Note that the approach is also entirely general to optomechanical systems, provided that the dynamics can be linearised.}. 
The outline of such a general approach to constructing the physical realisation given an $n$ degree-of-freedom generalised open oscillator goes as follows:
\begin{enumerate}
    \item First, the main synthesis theorem is used to split the $n$ degree-of-freedom oscillator into $n$ one degree-of-freedom oscillators which are connected in series, i.e.~the output of each oscillator is fed into the input of the next for example via a beamsplitter which is known as the series product~\cite{Gough2007}, and also a direct interaction Hamiltonian is produced coupling the oscillators. Often the series product connection is not needed as only one of the internal modes is coupled to the external continuum, i.e.~the operator $\hat{L}$ is only non-zero for one of the internal modes. The task is then to realise each of these one degree-of-freedom oscillators and the direct interaction Hamiltonian.
    \item For each one degree-of-freedom oscillator we do the following. First the scattering matrix can be realised as a static passive linear network using only beamsplitters and mirrors. Then, the general coupling operator of the form $\hat{L} = \alpha \hat{a} + \beta \hat{a}^\dag$ can be realised by indirectly coupling the mode $\hat{a}$ to the external continuum fields $\hat{u}$ and $\hat{y}$ via an auxiliary mode $\hat{b}$, which has sufficiently fast dynamics with coupling rate $\gamma$ such that it can be adiabatically eliminated from the final input-output relation. This auxiliary mode is coupled to the main mode via non-linear crystal (two-mode squeezing process) for the $\beta \hat{a}^\dagger$ term, and via a beamsplitter for the $\alpha \hat{a}$ term. These are related to the physical parameters via $\alpha = - \epsilon_2^* \sqrt{2/\gamma}$ where $\epsilon_2 = 2 \Theta e^{-i\Phi}$ where $\Theta$ is the beamsplitter mixing angle and $\Phi$ is the relative phase detuning introduced by the beamsplitter, and $\beta = \epsilon_1 \sqrt{2/\gamma}$ where $\epsilon_1$ is the effective pump intensity, shown to be equal to $c r / (2L)$ in Section.~\ref{appendix:single-pass}, where $r$ is the single-pass squeezing factor and $L$ is the cavity length. The resulting interaction Hamiltonian is $\hat{H}_{ab} = \hbar (\epsilon_1 \hat{a}^\dagger \hat{b}^\dagger + \epsilon_1^* \hat{a} \hat{b}) + \hbar (\epsilon_2 \hat{a}^\dagger \hat{b} + \epsilon_2^* \hat{a} \hat{b}^\dagger)$. Finally, the Hamiltonian $\hat{H}$ can be realised in the most general case as a detuned DPA (degenerate parametric amplifier), which can be implemented as a detuned cavity with a $\chi^{(2)}$ non-linear crystal with a pump frequency twice the laser carrier frequency $\omega_0$. Specifically to realise the Hamiltonian $\hat{H} = \hbar \Delta \hat{a}^\dagger \hat{a}^{\vphantom{\dagger}} + \hbar (\epsilon (\hat{a}^\dag)^2 + \epsilon^* \hat{a}^2)$ we use a cavity with resonant frequency $\omega_\text{cav} = \omega_0 + \Delta$ where $\omega_0$ is the laser carrier frequency and a non-linear crystal with effective pump intensity $\epsilon = c r / (2L)$ where again $r$ is the single-pass squeezing factor and $L$ is the cavity length.
    \item To implement the interaction Hamiltonian between each one degree-of-freedom oscillator we overlap the relevant internal modes of each oscillator via a non-linear crystal and/or beamsplitter depending on the interaction. The interaction Hamiltonian between modes $\hat{a}_k$ and $\hat{a}_l$ can be written in the form $\hat{H}_{kl} = \hbar(\epsilon_2 \hat{a}^\dag_k \hat{a}_l + \epsilon_2^* \hat{a}_k \hat{a}^\dag_l + \epsilon_1 \hat{a}^\dag_k \hat{a}^\dag_l + \epsilon_1 \hat{a}_k \hat{a}_l)$ where in this case the effective pump intensity is $-2i\epsilon_1$, and again $\epsilon_2 = 2 \Theta e^{-i\Phi}$ where $\Theta$ is the mixing angle and $\Phi$ is the relative phase difference between the two modes. In the unstable filter example below we do not need to do this as we already only have one internal degree of freedom.
\end{enumerate}

\begin{figure}[!b]
    \centering
    \includegraphics[width=0.95\linewidth]{graphics/figure1_paper.pdf}
    \caption{Flowchart showing the steps in constructing the physical realisation of a quantum filter; an active filter is used as an illustration. }
    \label{fig:realisation_steps}
\end{figure}

In summary, from the input-output transfer matrix we have developed the physical parameters describing the system. The procedure is summarised in Fig.\,\ref{fig:realisation_steps}.

\section{Illustrative example: an unstable filter}
\label{sec:unstable-filter}

To demonstrate the power of this approach, we will go beyond passive optical cavities by considering a non-trivial active filter for beating the universal gain-bandwidth product limit in resonant detection schemes~\cite{Wicht1997, Muller2000, Wise2004, Pati2007, Yum2013, Ma2015, Miao2015, David2015, Miao, Page2018, Bentley2019, Zhou2015, Zhou2018, Page2019, Shimazu2019}, specifically, the so-called unstable filter~\cite{Miao2015} which has a broadband anomalous (negative) dispersion. The gain-bandwidth product, also known as the Mizuno limit~\cite{JunMizuno}, states that the integral of the squared frequency-domain signal transfer function for a resonant detector is bounded purely by the energy stored in the detector. Therefore, when considering the quantum shot noise due to the input quantum vacuum, the detection bandwidth and peak sensitivity are inversely proportional which cannot be surpassed by changing any physical parameter other than the power in the detector. One way this can be surpassed is by using the broadband anomalous dispersion of the aforementioned unstable filter to partially negate the positive dispersion of the detector's signal cavity as first discussed in~\cite{Miao2015}. The filter also is unusual because it seemingly violates the Kramers-Kronig relations which imply that a stable anomalous dispersion filter without absorption (i.e.~with unity gain over the range of the anomalous dispersion) violates causality, however since this system is dynamically unstable this restriction does not apply~\cite{Kronig1926, Toll1956a, Doyle1990, Hirschorn2009}. As discussed in the procedure above we start with the frequency-domain transfer function of the unstable filter,
\begin{equation}
    G(s) = \frac{s - s_0}{s + s_0},
    \label{eq:transfer-function}
\end{equation}
where $s \equiv i \omega $ and $s_0 = \gamma_\text{neg}$ is a characteristic frequency quantifying the anomalous (negative) dispersion. An an example will now infer a physical realisation for this device using the above procedure.

First we note that since the transfer function is first order in frequency $s$, only one internal degree of freedom will be required for the minimal state-space realisation. Therefore the system state vector $\hat{x}$ will have two elements: $\mathbf{x}_1 = \hat{a}, \mathbf{x}_2 = \hat{a}^\dagger$ describing a single cavity mode, and similarly for the vectors $\mathbf{u}$ and $\mathbf{y}$ describing the input and output modes respectively. As a transfer-function 
matrix, Eq.~\eqref{eq:transfer-function} can be written as
\begin{equation}
    \mathbf{G}(s) = \frac{s-s_0}{s+s_0} \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix},
\end{equation}
which can be verified to satisfy the constraint \eqref{eq:j-j-unitary} and therefore a corresponding physical realisation can be found. To simplify the notation, we define a dimensionless $s$ (and the corresponding time) which is normalised with respect to $s_0/2=\gamma_{\rm neg}/2$ (a factor of $2$ for convenience), namely $s\rightarrow (s_0/2)s$ where $s$ is now dimensionless.
A corresponding minimal but not physically realisable state-space model is given by,
\begin{align}
	\begin{bmatrix} \dot{\hat{a}} \\ \dot{\hat{a}}^\dagger \end{bmatrix}
	&=
	\begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}
	\begin{bmatrix} \hat{a} \\ \hat{a}^\dagger \end{bmatrix}
	+
	\begin{bmatrix} \hat{u} \\ \hat{u}^\dagger \end{bmatrix},\\
	\begin{bmatrix} \hat{y} \\ \hat{y}^\dagger \end{bmatrix}
	&=
	\begin{bmatrix} 4 & 0 \\ 0 & 4 \end{bmatrix} 
	\begin{bmatrix} \hat{a} \\ \hat{a}^\dagger \end{bmatrix}
	+
	\begin{bmatrix} \hat{u} \\ \hat{u}^\dagger \end{bmatrix}. 
\end{align}

The matrix $X$ that solves Eqs.~\eqref{eq:physical-realizability-X1} and~\eqref{eq:physical-realizability-X} is given by $X = -J / 4$, which can be written in the form $X = TJT^\dag$ with the matrix $T$ which transforms the above state-space model to the physically realisable one is given by,
\begin{equation}
	T = \frac{1}{2} \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}.
\end{equation}
The resulting dimensionless state-space model can be found by applying the similarity transformation as shown in Eq.~\eqref{eq:applying-similarity-transformation},
\begin{equation}
	\tilde{A} = \begin{bmatrix}2 & 0 \\ 0 & 2\end{bmatrix},\ 
	\tilde{B} = \begin{bmatrix}0 & 2 \\ -2 & 0\end{bmatrix},\ 
	\tilde{C} = \begin{bmatrix}0 & -2 \\ 2 & 0\end{bmatrix},\ \tilde{D} = I,
	\label{eq:state-space}
\end{equation}
which obey Eqs.\,\eqref{eq:physical-realizability-1} and \eqref{eq:physical-realizability-2} by construction and therefore is a physically realisable state-space model. Reversing the normalisation process gives us,
\begin{align}
	A &= \begin{bmatrix}s_0 & 0 \\ 0 & s_0\end{bmatrix},\ 
	B = \begin{bmatrix}0 & \sqrt{2 s_0} \\ -\sqrt{2 s_0} & 0\end{bmatrix},\nonumber\\
	C &= \begin{bmatrix}0 & -\sqrt{2 s_0} \\ \sqrt{2 s_0} & 0\end{bmatrix},\ D = \tilde{D},
	\label{eq:state-space-with-dimensions}
\end{align}

Eq.\,\eqref{eq:SKR} can now be used to calculate the scattering matrix, input-output coupling, and internal Hamiltonian for the unstable filter. We have
\begin{equation}
    S = I,\ \hat{L} = -\sqrt{2 s_0}\hat{a}^\dag,\ \hat{H}=0.
\end{equation}
This implies that there is no input scattering with $S = I$, and $\hat{L} = -2\hat{a}^\dagger$, and there is no detuning or internal squeezing of the cavity mode as $\hat{H}=0$.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{graphics/unstable-filter-integrated-in-setup.pdf}
    \caption{Diagram showing where the filter realisation (highlighted by the shaded box) would be integrated into a standard Michelson inteferometer, using a scheme similar to that proposed in Ref.\,\cite{Miao2015}, which improves the signal response at high frequencies via the negative dispersion compensating the positive dispersion of the arm cavities as discussed above.}
    \label{fig:integrated-in-setup}
\end{figure}

Since $\hat{H}$ and $S$ are trivial to implement we now implement the coupling operator $\hat{L} = -\sqrt{2 s_0}\hat{a}^\dag$ as discussed above. In this case we have $\alpha = 0$, and so just have the auxiliary mode coupled to the main cavity mode via a non-linear crystal. This auxiliary mode will later be adiabatically eliminated, however it makes the physical realisation more feasible as coupling two cavity modes via a parametric oscillator is more experimentally durable. Therefore we construct the physical realisation shown in Fig.\,\ref{fig:realisation_steps}, which can be integrated into an interferometer as shown in Fig.\,\ref{fig:integrated-in-setup}. 

The realisation simply consists of two tuned cavities (the main mode $\hat{a}$ and the auxiliary mode $\hat{b}$) coupled via a $\chi^{(2)}$ non-linear crystal, labelled OPO (optical parametric oscillator), pumped by a classical pump field, labelled pump. One of the cavities is coupled to the external fields. Specifically, we have
\begin{align}
	\hat{H}_{ab} & = -\hbar \sqrt{s_0\,\gamma} (\hat{a}^\dagger \hat{b}^\dagger + \hat{a}\hat{b})\,, 
	\label{eq:interaction-hamiltonian-2}\\
	\hat H_{\rm ext} & = - i \hbar \sqrt{\gamma} 
	(\hat b\,\hat c_{\rm ext}^{\dag} - \hat b^{\dag} \hat c_{\rm ext})\,.
	\label{eq:interaction-hamiltonian-external-continuum}
\end{align}
The interaction Hamiltonian $\hat H_{ab}$ describes the coupling of both cavity modes $\hat{a}$ and $\hat{b}$ via the OPO. As shown in Section.~\ref{appendix:single-pass}, the coupling rate $\sqrt{s_0\,\gamma}$ is equal to $rc/(2L_b)$, where $r$ is the  single-pass squeezing factor of the crystal and $L_b$ is the length of the auxiliary  cavity. As an order of magnitude estimate for implementation in a laser interferometer with arm length of $L_\text{arm} = \SI{4}{\kilo\metre}$ (where $s_0 \equiv \gamma_\text{neg} = c/L_\text{arm}$~\cite{Miao2015}), the required squeezing factor is
\begin{equation}
    r = 7.7\times 10^{-5} \sqrt{\frac{T_b}{\SI{100}{\ppm}}} \sqrt{\frac{L_b}{\SI{24}{\centi\metre}}\vphantom{\frac{T_b}{\SI{10}{\ppm}}}} \sqrt{\frac{\SI{4}{\kilo\metre}}{L_\text{arm}}\vphantom{\frac{T_b}{\SI{10}{\ppm}}}}.
\end{equation}
The Hamiltonian $\hat H_{\rm ext}$ describes the coupling
between the auxiliary mode $\hat b$ and the external
continuum field $\hat c_{\rm ext}$, which is related 
to the input and output operators via $\hat u\equiv \hat c_{\rm ext}(t=0_-)$ and $\hat y\equiv \hat c_{\rm ext}(t=0_+)$~\cite{Walls2008, Chen2013}. The coupling rate $\gamma$ is defined as $T_b c / (4 L_b)$
where $T_b$ is the input mirror transmissivity.
The negative dispersion transfer function shown in Eq.~\eqref{eq:transfer-function} can then be recovered 
by solving the resulting Heisenberg equations of motion 
in the frequency domain, and then applying the approximation $\gamma \gg \omega$, the so-called ``resolved-sideband regime'', which effectively adiabatically eliminates $\hat{b}$~\cite{Nurdin2009}.

In Section.~\eqref{appendix:losses}, we include the 
effect of optical loss for the realistic implementation. 
We found that the noise contribution from the auxiliary cavity loss is insignificant compared to the contribution from the $\hat{a}$ cavity loss. The resulting 
input-output relation including the optical loss is given by
\begin{equation}
    \hat{y}(s) \approx \frac{\omega + i(\gamma_a^\epsilon + s_0)}{\omega + i(\gamma_a^\epsilon - s_0)} \hat{u}(s) + \frac{2 \sqrt{s_0\,\gamma_a^\epsilon}}{\omega + i(\gamma_a^\epsilon - s_0)} \hat{n}_a^\dagger(s),
\end{equation}
where $\gamma_a^\epsilon = \epsilon_a c / (4 L_a)$ with $\epsilon_a$ being the total optical loss in the $\hat{a}$ cavity and $L_a$ being the cavity length, and $\hat{n}_a$ is 
the corresponding vacuum noise process. The distortion of the transfer function due to $\gamma_a^\epsilon$ is on the order of $\gamma_a^\epsilon / s_0$, while the noise term is on the order of $\sqrt{\gamma_a^\epsilon / s_0}$ and is therefore more significant. 

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{graphics/required_losses.pdf}
    \caption{Required total $\hat{a}$ cavity loss $\epsilon_a$ as a function of ratio of $\hat{a}$ cavity length to arm cavity length $L_a / L_\text{arm}$ for the cases where the noise power contribution at $\omega = 0$ due to $\hat{n}_a$ is a tenth of that of the signal power (blue line) and a hundred (orange line).}
    \label{fig:required_losses}
\end{figure}

The above input-output relation takes the same form as the optomechanical case~\cite{Miao2015} if we view $\hat n_a$ 
as the thermal noise of the mechanical oscillator. 
In contrast, in this case the loss $n_a$ is sourced by the quantum vacuum and so it only has vacuum fluctuations, equivalent to a mechanical oscillator at environmental temperature $T_\text{env} = 0$. Therefore the strict thermal requirements of the optomechanical unstable filter are avoided. Instead vacuum fluctuations are injected due to losses in the mirrors and the non-linear crystal. The required loss to achieve low noise as a function of $\hat{a}$ cavity length is shown in Fig.~\ref{fig:required_losses}. As we can see, given an interferometer arm length of $L_\text{arm} = \SI{4}{\kilo\metre}$, a loss per unit length of $\epsilon_a / L_a = \SI{2.5}{\ppm\per\metre}$ is required to achieve a $1/10$ noise contribution, which is already achievable with state-of-the-art optics for the free space part of the cavity~\cite{Isogai2013,Oelker2016a}. The loss of a \SI{1}{\centi\metre} crystal is between around 15--40 \si{\ppm} \cite{Bogan2015}, so taking a middle value of around $\SI{25}{\ppm}$ then to achieve a loss per unit length of $\SI{2.5}{\ppm\per\metre}$ the cavity length must be at least \SI{10}{\metre}. Further assuming losses of \SI{10}{\ppm} for each side of the anti-reflective coating of the crystal drives the minimum cavity length up to \SI{18}{\metre}.

\section{Discussion}
\label{sec:discussion}

In general the physical realisation produced may have unstable internal dynamics, as is the case for the unstable filter above. Since the single-mode approximation is used, the stabilising controller previously derived in~\cite{Miao2015} can be used for each unstable degree of freedom.

In addition to 
realising quantum filters with a known transfer function, this approach can also be used to design the optimal high-precision measurement devices, where the optimality is based upon the quantum Cram\'{e}r-Rao bound\,~\cite{Helstrom1967,Holevo2011,Braunstein1996,Giovannetti2011,Tsang2011,Miaoa}. We can view the
entire measurement device as a $N$ degree-of-freedom quantum filter, and then tune the filter parameters so as to minimise the quantum Cram\'{e}r-Rao bound. Therefore, we can construct the most sensitive possible $n$ degree-of-freedom
measurement device. This opens up a new paradigm of designing and optimising measurement devices and is worthy of being further explored.

\section{Acknowledgements}

We would like to thank Rana Adhikari, Denis Martynov, Naoki Yamamoto, LIGO AIC, and QNWG for fruitful discussions. 
J.B. is supported by STFC and School of Physics and
Astronomy at the University of Birmingham. J.B. and H.M.
acknowledge the additional support from the Birmingham
Institute for Gravitational Wave Astronomy.
H.M. has also been supported by UK STFC Ernest Rutherford 
Fellowship (Grant No. ST/M005844/11). Y.C. is supported by the Simons Foundation (Award Number 568762), and the National Science Foundation, through Grants PHY-1708212 and PHY-1708213. 

