Graph-based analysis of nonreciprocity in coupled-mode systems

Reciprocity is the symmetry of a physical system with respect to the exchange of a source and detector. For example, in systems which demonstrate nonreciprocal behavior the transmission loss and/or phase delay depends on the direction of propagation [1]. In general, reciprocal symmetry can be violated in multi-mode devices and the most well-known examples are ferrite-based circulators and isolators [2, 3]. Such devices are useful in several practical contexts: circulators for instance are commonly employed in both optical [4] and cryogenic microwave systems [5–7] to reduce reflections and, more importantly, to protect from amplifier backaction on the device-under-test. Likewise, in superconducting quantum information, circulators are often used to direct the amplified reflected signal from lumped-element superconducting parametric amplifiers.

Since the most common schemes for producing nonreciprocity involve the use of high magnetic fields, they are often incompatible with the goal of integrating several measurement circuits into smaller and smaller physical volumes while also limiting the magnetic flux near sensitive superconducting circuits [7]. In the superconducting quantum information and microwave engineering fields, this has motivated a push to understand how nonreciprocity can be generated using alternative methods [7–11] and this general pursuit has been paralleled by similar efforts in optics [12–17]. Although these different ideas are embodied in different physical implementations, they all share a common mathematical description based in coupled-mode theory [18]. In this paper, we develop a scheme for determining the scattering matrix of an arbitrary coupled-mode system in terms of directed graphs. In doing so, we show that nonreciprocity is a consequence of asymmetric interference between different connecting paths in a graph, as well as judiciously-placed dissipation in ancillary modes. We emphasize that while a direct solution of the coupled-mode equations can always be found by means of standard mathematical techniques, the reverse problem of synthesizing a multi-mode system with desired properties (gain, impedance match, directionality) can be much more difficult. With the approach we describe in this work, the synthesis of a multi-mode device is translated into the process of building a directed graph whose edges are subject to specific conditions. Since this approach is agnostic to physical implementation, it may benefit similar efforts based on optical, mechanical, and hybrid systems. As many of these systems are parametrically coupled systems, we frame this work in the context of frequency conversion and amplification processes.

This paper is organized as follows—in section 2 we introduce normalization and matrix conventions for extending the coupled-mode equations of motion to several modes. This facilitates the definition of reciprocity in parametrically coupled systems with three or more modes. In section 3 we visualize the matrix representation of the equations of motion using directed graphs and show how the scattering matrix can be computed by use of subgraphs which connect the input/output modes. We then use this picture to discuss the mechanism for nonreciprocity in general parametrically-coupled systems. In sections 4 and 5 we show how directed graphs can be used as both a synthesis and analysis tool. As examples, we synthesize a novel three-mode parametric amplifier in which gain, directionality, and impedance match are present (section 4). Finally, we analyze the isolation properties of a complicated system, the dc-SQUID, by casting it as a ten-mode scattering problem (section 5) based on earlier work by Kamal et al [19]. We conclude the paper with a discussion of how this graphical approach might be effectively applied in future work.

Resonant and parametrically coupled modes are typically described by the Langevin equations of motion [20]. Since in this work we are primarily focused on the topology of the coupling in a given system, we utilize a normalized matrix form of these equations and build a graph representation of multimode coupling dynamics. This approach has an additional advantage in that it allows one to deal with systems with several modes in different physical oscillator incarnations as one may encounter in hybrid electro- and opto-mechanical systems. While the transition from the time-domain equations of motion to their matrix representation is given in appendix A, in this section we present only the matrix variables that are necessary for building the graphs.

We begin by considering a set of N_r coupled resonators, described by a set of N_m internal mode amplitudes. In general, we allow the number of modes to exceed the number of resonators ${{N}_{m}}\geqslant {{N}_{r}}$ such that several modes may reside within a single resonator. This system can be a set of parametrically or resonantly coupled oscillators, but we are not going to assume a particular physical implementation for now, in order to derive results that are as general as possible. We will call ${{\omega }_{j}}$ and ${{\gamma }_{j}}$ the natural oscillation frequency and total dissipation rate for the resonator in which mode j resides. The complex coupling rate between modes j and k will be denoted by g_jk. In order to describe both mode frequency conversion and amplification, we divide the modes into two subsets of p internal mode amplitudes ${{b}_{1...p}}$ and $q={{N}_{m}}-p$ internal mode amplitudes $b_{p+1...{{N}_{m}}}^{\dagger }$ in the frequency basis, with each set of modes corresponding to the driven response at frequencies $\omega _{1...p}^{s}$ and $-\omega _{p+1:{{N}_{m}}}^{s}$. We assume that ${{g}_{jk}}=g_{kj}^{*}$ in the case of frequency conversion between modes j and k (for $1\leqslant j,k\leqslant p$ or $p+1\leqslant j,k\leqslant {{N}_{m}}$) and ${{g}_{jk}}=-g_{kj}^{*}$ in the case of parametric amplification (for $1\leqslant j\leqslant p$ and $p+1\leqslant k\leqslant {{N}_{m}}$, or $1\leqslant k\leqslant p$ and $p+1\leqslant j\leqslant {{N}_{m}}$). For example, the coupling matrix for a two-mode frequency converter or resonantly coupled oscillator system corresponds to the case $p=2,q=0$, while a conventional parametric amplifier couples two modes, one at a positive frequency (p = 1) and the other at a negative frequency (q = 1) [21]. With this general prescription for defining our mode basis, we can perform an input/output analysis of any system of coupled resonators and calculate the corresponding scattering matrix connecting the vector of input fields bⁱⁿ to the vector of output fields b^out (see appendix A). The scattering matrix can in general be expressed as

$$\begin{eqnarray}&&S={\rm i}\frac{1}{{{\gamma }_{M}}}K{{M}^{-1}}K-\mathbb{I},\end{eqnarray} \tag{ 1 }$$

where M is a ${{N}_{m}}\times {{N}_{m}}$ matrix given by

$$\begin{eqnarray}M=\left[ \begin{array}{ccccccccccccccc} {{\Delta }_{1}} & \cdots & {{\beta }_{1,p}} & {{\beta }_{1,p+1}} & \cdots & {{\beta }_{1,{{N}_{m}}}} \\ \vdots & \ddots & \vdots & \vdots & {} & \vdots \\ \beta _{1,p}^{*} & \cdots & {{\Delta }_{p}} & {{\beta }_{p,p+1}} & \cdots & {{\beta }_{p,{{N}_{m}}}} \\ -\beta _{1,p+1}^{*} & \cdots & -\beta _{p,p+1}^{*} & -\Delta _{p+1}^{*} & \cdots & -\beta _{p+1,{{N}_{m}}}^{*} \\ \vdots & {} & \vdots & \vdots & \ddots & \vdots \\ -\beta _{1,{{N}_{m}}}^{*} & \cdots & -\beta _{p,{{N}_{m}}}^{*} & -{{\beta }_{p+1,{{N}_{m}}}} & \cdots & -\Delta _{{{N}_{m}}}^{*} \\ \end{array} \right]\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{{\Delta }_{j}}=\frac{{{1}_{{}}}}{{{\gamma }_{M}}}\left( (\omega _{j}^{s}-{{\omega }_{j}})+{\rm i}\frac{{{\gamma }_{j}}}{2} \right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\beta }_{jk}}=\frac{{{g}_{jk}}}{2{{\gamma }_{M}}},\end{eqnarray} \tag{ 4 }$$

where we have defined an overall normalization prefactor

$$\begin{eqnarray}&&{{\gamma }_{M}}=\sqrt[{{N}_{m}}]{\mathop{\prod }\limits_{\ell =1}^{{{N}_{m}}}{{\gamma }_{\ell }}},\end{eqnarray} \tag{ 5 }$$

and environmental coupling matrix

$$\begin{eqnarray}K=\left[ \begin{array}{ccccccccccccccc} \sqrt{\gamma _{1}^{{\rm ext}}} & {} & 0 \\ {} & \ddots & {} \\ 0 & {} & \sqrt{\gamma _{{{N}_{m}}}^{{\rm ext}}} \\ \end{array} \right].\end{eqnarray} \tag{ 6 }$$

In general ${{\gamma }_{j}}\geqslant \gamma _{j}^{{\rm ext}}$, with the equal sign if the internal dissipation rate of mode j is zero. The diagonal elements of M are complex normalized detunings between the driven response frequencies and the natural resonator frequencies and include the mode dissipation rates, while the off-diagonal elements are normalized coupling coefficients. In principal, the normalizations above are unnecessary, but they allow one to cast the Langevin equations of motion (see appendix A) in a simple matrix form (as M, the 'Langevin matrix'), emphasizing the underlying structure of the connections between modes. We show below that this structure, as revealed in the graph representation of M, allows one to draw immediate conclusions about the (non)reciprocity in multimode coupled systems.

2.1. A formal definition of reciprocity

In order to proceed, we must formally define reciprocity within the context of several modes, possibly oscillating at different frequencies (as found in parametrically coupled systems). The scattering matrix (equation (1)) describes a fully reciprocal system (transmission between any two modes is reciprocal) if it obeys, in a suitable mode basis, the general constraint [22, 23]

$$\begin{eqnarray}&&{{S}^{T}}={{U}_{\phi }}SU_{\phi }^{\dagger }.\end{eqnarray} \tag{ 7 }$$

Here the matrix ${{U}_{\phi }}$ corresponds to a set of phaseshifts in the mode basis

$$\begin{eqnarray}{{U}_{\phi }}=\left[ \begin{array}{ccccccccccccccc} {{{\rm e}}^{{\rm i}{{\phi }_{1}}}} & {} & 0 \\ {} & \ddots & {} \\ 0 & {} & {{{\rm e}}^{{\rm i}{{\phi }_{{{N}_{m}}}}}} \\ \end{array} \right],\end{eqnarray} \tag{ 8 }$$

where we can assume without loss of generality that $\mathop{\sum }\limits_{\ell =1}^{N}{{\phi }_{\ell }}=0$. The matrix ${{U}_{\phi }}$ therefore corresponds to an arbitrary redefinition of the mode phases. This is similar to the concept of a common mode phaseshift, but generalized to several modes.

A system that violates condition (7) is nonreciprocal and is characterized by an asymmetric phase and/or amplitude transmission coefficient between at least one pair of modes. We note that equation (7) differs from the definition commonly found in textbooks [3], ${{S}^{T}}=S$. In fact, equation (7) is valid for more general multimode parametric systems in which the mode frequencies can differ. As such, it is a gauge-invariant definition of nonreciprocity² [22, 23]. Since the environmental coupling matrix K in equation (6) is diagonal, it follows from equation (1) that reciprocity is guaranteed if and only if the Langevin matrix, M, obeys the same similarity transformation as in equation (7)

$$\begin{eqnarray}&&{{M}^{T}}={{U}_{\phi }}MU_{\phi }^{\dagger },\end{eqnarray} \tag{ 9 }$$

for some set of phaseshifts $\{{{\phi }_{\ell }}\}$. Above, we noted that all resonant and parametric systems considered here are constrained to couplings which are related by pairwise (anti)-conjugation, i.e., ${{\beta }_{jk}}=\pm \beta _{kj}^{*}$. Placing this constraint in equation (9) we obtain

$$\begin{eqnarray}&&{{\phi }_{j}}-{{\phi }_{k}}=-2\angle {{\beta }_{jk}}+{{n}_{jk}}\pi ,\end{eqnarray} \tag{ 10 }$$

where $\angle {{\beta }_{jk}}$ is the phase of the coupling coefficient ${{\beta }_{jk}}$ between modes j and k and n_jk is an integer. The system of equation (10) has at most a set of ${{N}_{m}}({{N}_{m}}-1)/2$ independent relations (if all off-diagonal elements of M are nonzero) in ${{N}_{m}}-1$ unknown phases $\{{{\phi }_{\ell }}\}$. The system is therefore overdetermined and the reciprocity condition (equation (7)) is, for parametrically coupled systems, satisfied only for specific choices of coupling phases.

The conditions for equation (10) to have a unique solution, and therefore for the system to be reciprocal, can be visualized by drawing a directed graph representing the matrix M (see figure 1). The graph has N_m vertices and we associate a weighted directed branch ('edge') between vertices j and k to every element ${{\beta }_{jk}}$ in the matrix M. The diagonal elements ${{\Delta }_{j}}$ and $-\Delta _{k}^{*}$ of M are the weights of self-loops in the graph. In other words, all of the coupling terms in M are shown as directed branches, and all mode detuning and dissipation is represented by self-loops in the graph.

Zoom In Zoom Out Reset image size

Within this representation, the conditions for a unique solution for the entire system can be obtained by summing equation (10) along every closed loop L, yielding a simple result

$$\begin{eqnarray}&&\mathop{\sum }\limits_{{{\beta }_{jk}}\in L}\angle {{\beta }_{jk}}={{k}_{L}}\pi ,\end{eqnarray} \tag{ 11 }$$

where k_L is an integer. Equation (11), reminiscent of Kirchoff's first law in electrical circuits, is analogous to the reciprocity condition in Jaynes–Cummings lattices described in [11, 24]. We stress, however that equation (11) is applicable to any linear coupled-mode system bearing arbitrarily complex connections and is not confined to the nearest-neighbor, planar graph couplings shown here. Moreover, in the case of parametrically coupled devices in particular, the vertices of our graphs represent modes that may physically correspond to energy that is differentiated by both space and frequency. As such, the graphs are a very general representation of a coupled mode system and can describe coupling of energy at several different frequencies as well as physical ports.

In order for a system to be nonreciprocal, the condition (11) must be violated. This, however, only expresses a minimal requirement for phase nonreciprocity. We must further define a condition for which a system will also demonstrate amplitude nonreciprocity. Specifically, we may ask which conditions allow the amplitude isolation ratio $\mid {{S}_{jk}}/{{S}_{kj}}\mid$ to be different than 1. From (1) it follows that the scattering coefficient S_jk scales with the corresponding element of the inverse Langevin matrix, $M_{jk}^{-1}$. The problem of calculating the scattering coefficients is therefore reduced to the problem of inverting the Langevin matrix which, as we show below, can also be conveniently cast in terms of directed graphs.

The inverse Langevin matrix can be expressed in terms of the transpose of its cofactor matrix [25]

$$\begin{eqnarray}&&{{M}^{-1}}=\frac{{{C}^{T}}}{\mid M\mid }.\end{eqnarray} \tag{ 12 }$$

The scattering coefficient is then given by

$$\begin{eqnarray}&&{{S}_{jk}}={\rm i}\frac{\sqrt{\gamma _{j}^{{\rm ext}}\gamma _{k}^{{\rm ext}}}}{{{\gamma }_{M}}}\frac{{{C}_{kj}}}{\mid M\mid }-{{\delta }_{jk}},\end{eqnarray} \tag{ 13 }$$

where ${{\delta }_{jk}}$ is the Kronecker delta function. The elements of the cofactor matrix C_kj are equal to the determinant of the matrix obtained from M after removing the kth row and the jth column (i.e., the kj minor), multiplied by ${{(-1)}^{j+k}}$ [25].

The determinant of M and the cofactor C_kj can be easily calculated from the graph in figure 1(a) by the following procedure (see [26, 27]). Given a graph, we call a subgraph G a 'permutation subdigraph' if every vertex in G has one and only one outgoing edge and one and only one incoming edge [27]. The total weight w(G) of a subgraph G is the product of the weights of all of its edges times ${{(-1)}^{c+{{N}_{m}}}}$, where c is the number of cycles in G (by cycle we mean a sequence of distinct connected vertices, where the first and last vertex coincide). The determinant of the matrix M is the sum of weights of every possible permutation subdigraph of the graph associated to M [26, 27], as shown in figure 1(b).

The cofactor matrix C_kj is the coefficient of the M_kj element of the matrix M in the Laplace expansion of the determinant $\mid M\mid$. Therefore, in order to calculate C_kj, we need to compute only the permutation subdigraphs in the expansion of $\mid M\mid$ that contain the kj edge, with the kj edge removed. In other words, for every path p from vertex k to vertex j, we consider the collection of all the permutation subdigraphs G^p_r of the graph of M containing p, with the kj edge removed as in figure 1(c). The cofactor C_kj is then given by

$$\begin{eqnarray}&&{{C}_{kj}}=-\mathop{\sum }\limits_{p,r}w(G_{r}^{p}),\end{eqnarray} \tag{ 14 }$$

where p varies over all the paths connecting vertex j to vertex k and r varies over all the subdigraphs containing the path p. The minus sign is due to the fact that, by removing the kj edge, we reduced the number of cycles in the subdigraph by 1. An example of this procedure is outlined in figure 1(b). The graphs $\{G_{r}^{p}\}$ in the expression for the cofactor elements have a simple physical interpretation: they represent the possible scattering mechanisms that connect modes j and k.

Since, from (13), the ratio between the scattering coefficients ${{S}_{jk}}/{{S}_{kj}}$ is equal to the ratio between the corresponding cofactors ${{C}_{kj}}/{{C}_{jk}}$, the problem of analyzing isolation in a given multi-mode system can be reduced to evaluating the differences between the specific graphs that connect j to k and vice versa.

As a more concrete example, we apply this procedure to a simple three-mode system in which all modes are coupled via frequency conversion (${{\beta }_{jk}}=\beta _{kj}^{*}$) and are therefore described by the following matrix (see figure 2),

$$\begin{eqnarray}M=\left[ \begin{array}{ccccccccccccccc} {{\Delta }_{1}} & {{\beta }_{12}} & {{\beta }_{13}} \\ \beta _{12}^{*} & {{\Delta }_{2}} & {{\beta }_{23}} \\ \beta _{13}^{*} & \beta _{23}^{*} & {{\Delta }_{3}} \\ \end{array} \right],\end{eqnarray} \tag{ 15 }$$

where ${{\Delta }_{j}}$ and ${{\beta }_{jk}}$ are normalized detunings and couplings as before. The isolation between ports 1 and 2 can be calculated by use of the graphs representing the cofactors (see figures 2(b) and (c)) yielding the expression

$$\begin{eqnarray}&&{{I}_{12}}=\frac{{{S}_{12}}}{{{S}_{21}}}=\frac{{{C}_{21}}}{{{C}_{12}}}=\frac{-\beta _{21}^{*}{{\Delta }_{3}}+{{\beta }_{32}}{{\beta }_{13}}}{-{{\beta }_{21}}{{\Delta }_{3}}+\beta _{32}^{*}\beta _{13}^{*}}=\frac{-2g_{21}^{*}(\omega _{3}^{S}-{{\omega }_{3}}+{\rm i}{{\gamma }_{3}}/2)+{{g}_{32}}{{g}_{13}}}{-2{{g}_{21}}(\omega _{3}^{S}-{{\omega }_{3}}+{\rm i}{{\gamma }_{3}}/2)+g_{32}^{*}g_{13}^{*}},\end{eqnarray} \tag{ 16 }$$

where, in the last line, we have removed all normalizations.

Zoom In Zoom Out Reset image size

From (16) we observe that if the total dissipation of mode '3', ${{\gamma }_{3}}=0$, then ${{S}_{12}}=S_{21}^{*}$ and there is no isolation between modes 1 and 2. In this case mode 3 serves as an internal lossless mode and the system behaves as a simple reciprocal two-mode system. Here, dissipation plays a crucial role in violating reciprocity, as it breaks the symmetry between scattering elements which are otherwise related by complex conjugation. Assuming finite dissipation in mode 3, we find the condition for maximum isolation ${{I}_{12}}=0$,

$$\begin{eqnarray}&&\frac{{{g}_{13}}{{g}_{32}}}{2g_{21}^{*}}=\omega _{3}^{s}-{{\omega }_{3}}+{\rm i}{{\gamma }_{3}}/2.\end{eqnarray} \tag{ 17 }$$

In other words, to achieve maximum isolation, the sum of the phases of the couplings along the loop connecting the three modes, ${{\phi }_{13}}+{{\phi }_{21}}+{{\phi }_{32}}$, must be equal to the detuning angle, ${{{\rm tan} }^{-1}}[{{\gamma }_{3}}/(2(\omega _{3}^{s}-{{\omega }_{3}}))]$. At zero detuning in mode 3, $\omega _{3}^{s}={{\omega }_{3}}$, this loop phase must be $\pi /2$ for perfect isolation. Equivalently, one can interpret this isolation condition, equation (17), as interference between the two possible coupling paths connecting modes 1 to 2 in the graph (figure 2). Although we focus on the isolation here, condition (17) also yields unit transmission in the forward direction, $\mid {{S}_{21}}\mid =1$, for $\mid {{g}_{12}}\mid {{\gamma }_{3}}=1$. In fact, one can continue to build a three-mode circulator by finding similar constraints for the 2↔3 and 1↔3 coupling.

Traditional single-stage lumped-element parametric amplifiers, such as the one shown in figure 3(a), operate in reflection-mode such that a circulator is needed in order to achieve unidirectionality. In [28] a single-stage amplifier with two parametric pumps at the same frequency was introduced, but in addition to forward gain it also had unity reverse gain. The graph methodology introduced here allows us to revisit the problem of gain and isolation from another perspective. In particular, our discussion above identified the importance of introducing a multi-edge loop to obtain nonreciprocity. Here we analyze a closed-loop amplifier which we term the 'delta' amplifier, shown schematically in figure 3(b) both as a device concept and as its corresponding graph. This amplifier can provide forward gain and zero reverse gain as discussed below.

Zoom In Zoom Out Reset image size

As drawn, the delta amplifier requires three independent pumps to generate three pairs of coupling edges $\{{{\beta }_{jk}}\}$. Two of these pairs create amplification (i.e., ${{\beta }_{23}}=-\beta _{32}^{*}$, ${{\beta }_{12}}=-\beta _{21}^{*}$), while the third pair couples through frequency conversion (${{\beta }_{31}}=\beta _{13}^{*}$). One concept for how this might be implemented is shown schematically in figure 3(b) as three resonators that share a common current anti-node in which a flux-driven SQUID is operated as a parametrically modulated inductor. In principle, one might drive two of these couplings using the same pump frequency³ [29], but for clarity we discuss only the case in which each of the three couplings is driven independently, which has the advantage of increased flexibility due to the increased number of degrees of freedom. Regardless of particular implementation, we can recognize the presence of a closed loop connecting all three modes, as well as the integration of amplification as a coupling process. From this bare description we recognize that with three independent pumps we are constrained to optimizing only three of the nine possible scattering elements. In laboratory use, we can choose to optimize for amplification in the 'forward' direction ($\mid {{S}_{21}}\mid \gt 1$), isolation in the 'reverse' direction (${{S}_{12}}\sim 0$), and low input reflection coefficient (input match, i.e., ${{S}_{11}}\sim 0$) to prevent unwanted reflections of signals propagating away from a device-under-test.

We begin by writing down the expressions for the forward and reverse transmission from equation (13),

$$\begin{eqnarray}&&{{S}_{21}}={\rm i}\frac{\sqrt{\gamma _{1}^{{\rm ext}}\gamma _{2}^{{\rm ext}}}}{\sqrt[3]{{{\gamma }_{1}}{{\gamma }_{2}}{{\gamma }_{3}}}}\frac{\beta _{12}^{*}{{\Delta }_{3}}+{{\beta }_{31}}{{\beta }_{23}}}{\mid M\mid },\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{S}_{12}}={\rm i}\frac{\sqrt{\gamma _{1}^{{\rm ext}}\gamma _{2}^{{\rm ext}}}}{\sqrt[3]{{{\gamma }_{1}}{{\gamma }_{2}}{{\gamma }_{3}}}}\frac{-{{\beta }_{12}}{{\Delta }_{3}}-\beta _{23}^{*}\beta _{31}^{*}}{\mid M\mid },\end{eqnarray} \tag{ 19 }$$

where $\mid M\mid$ is the determinant of the coupling matrix represented in the graph shown in figure 3(b). Assuming ideal isolation, ${{S}_{12}}=0$, we find an expression similar to equation(17),

$$\begin{eqnarray}&&\frac{\beta _{23}^{*}\beta _{31}^{*}}{{{\beta }_{12}}}=-{{\Delta }_{3}},\end{eqnarray} \tag{ 20 }$$

which provides an overall constraint on the coupling amplitudes and phases as a function of the complex detuning of mode 3. Under this primary constraint, we then set the input reflection coefficient to zero,

$$\begin{eqnarray}&&{{S}_{11}}={\rm i}\frac{\gamma _{1}^{{\rm ext}}}{\sqrt[3]{{{\gamma }_{1}}{{\gamma }_{2}}{{\gamma }_{3}}}}\frac{\mid {{\beta }_{23}}{{\mid }^{2}}-\Delta _{2}^{*}{{\Delta }_{3}}}{\mid M\mid }-1=0,\end{eqnarray} \tag{ 21 }$$

giving a constraint on the determinant

$$\begin{eqnarray}&&\mid M\mid ={\rm i}\frac{\gamma _{1}^{{\rm ext}}}{\sqrt[3]{{{\gamma }_{1}}{{\gamma }_{2}}{{\gamma }_{3}}}}(\mid {{\beta }_{23}}{{\mid }^{2}}-\Delta _{2}^{*}{{\Delta }_{3}}).\end{eqnarray} \tag{ 22 }$$

Taken together, these constraints yield an expression for the forward transmission gain

$$\begin{eqnarray}&&\sqrt{G}=\mid {{S}_{21}}\mid =\sqrt{\frac{\gamma _{2}^{{\rm ext}}}{\gamma _{1}^{{\rm ext}}}}\frac{2\mid {{\beta }_{12}}{\rm Im}({{\Delta }_{3}})\mid }{\mid \mid {{\beta }_{23}}{{\mid }^{2}}-\Delta _{2}^{*}{{\Delta }_{3}}\mid }.\end{eqnarray} \tag{ 23 }$$

If external coupling dominates, $\{{{\gamma }_{k}}\}=\{\gamma _{k}^{{\rm ext}}\}$, in the limit of small detunings $(\omega _{k}^{s}-{{\omega }_{k}})\ll {{\gamma }_{k}}$ the gain reduces to

$$\begin{eqnarray}&&\sqrt{G}\approx \sqrt{\frac{{{\gamma }_{2}}}{{{\gamma }_{1}}}}\frac{2\mid {{g}_{12}}\mid {{\gamma }_{3}}}{\mid {{g}_{23}}{{\mid }^{2}}-{{\gamma }_{2}}{{\gamma }_{3}}}.\end{eqnarray} \tag{ 24 }$$

Note that this form gives the gain at a fixed frequency where the input match is perfect (equation (21)). By substituting equations (18) and (20) in equation (21) at resonance we can determine the three coupling rates. In the limit of high gain, $\mid {{g}_{jk}}\mid \approx \sqrt{{{\gamma }_{j}}{{\gamma }_{k}}}$, with the exact value depending on the gain G.

We have calculated the resulting S-parameters for the delta amplifier in figure 4 under various constraints on the pump detunings for the case of uniform dissipation rates ${{\gamma }_{k}}=\gamma$. In practice, one might fix the pump frequencies, locking the three detunings together, and measure the S-parameters shown in figure 4(d). For this example calculation, we selected values for the pump amplitudes to fix the nominal power gain to ∼20 dB while obtaining significant isolation. The device bandwidth is determined according to equation (18) that fixes the gain-bandwidth product for specific dissipation rates. We note here that the device bandwidth is comparable to traditional two-mode superconducting lumped parametric amplifiers, while the center frequency could be tuned by changing the frequencies of the pumps. Finally, we find that the calculated output S₂₂ shows reflection gain at the output port, and is consistent with the limitations of having only three parameters (the independent pumps) to optimize three of the four scattering parameters between ports 1 and 2. However, this might not be a limitation in practical use where one is more concerned with the backaction of the amplifier on a quantum circuit on the input side (port 1).

Zoom In Zoom Out Reset image size

4.1. Added noise of the three-mode amplifier

We computed the equivalent input added noise ${{\bar{n}}_{{\rm add}}}$ of the delta amplifier to check that the device can approach the standard quantum limit. We assumed that half a photon of noise is injected into every port and calculated the total output noise spectral density from the scattering parameters. The contribution from ports 2 and 3, divided by the device power gain G, constitutes the equivalent added input noise, as outlined in [30]. For a symmetric device, in the limit of large gain and at zero detuning, the result is:

$$\begin{eqnarray}&&{{\bar{n}}_{{\rm add}}}(0)=\frac{\mid {{S}_{22}}(0){{\mid }^{2}}+\mid {{S}_{23}}(0){{\mid }^{2}}}{2G}\approx \frac{\mid {{S}_{22}}(0){{\mid }^{2}}}{2G}\approx \frac{1}{2}{{\left( 1+\frac{1}{\sqrt{G}} \right)}^{2}},\end{eqnarray} \tag{ 25 }$$

where we can neglect $\mid {{S}_{23}}\mid$ as a consequence of the isolation condition. The results of the numerical calculations for nonzero detuning are shown in figure 4. The fact that the delta amplifier is quantum-limited may be surprising at first, as one might expect the third mode to contribute extra noise beyond the quantum limit. However the third mode acts only as an effective dissipative coupling between modes 1 and 2 [30]. We note that the parametric amplification scheme that is the subject of [30] constitutes a three mode device in which the added noise was determined to be at the standard quantum limit. In fact, the delta amplifier is similar to that device except for the addition of an amplification branch which serves to create a closed loop in its graph description, introducing the necessary interference for directionality.

While the previous discussion focused on using a graphical approach to synthesize a novel amplifier, it can also be used to identify aspects of other systems which enable both gain and directionality. In this section we use graphs to analyze a complicated multi-mode system, the dc-SQUID amplifier. dc-SQUIDs are a type of superconducting amplifier with noise temperatures approaching the standard quantum limit [31–33]. Since dc-SQUIDs also have inherently low power dissipation (in addition to high gain and directionality [34]), they have been the subject of several applied superconductivity efforts to produce reliable amplifiers for microwave quantum information measurements. In [19], Kamal et al showed that the directionality of the SQUID is a consequence of multiple parametric amplification and frequency conversion processes due to frequency mixing between the input microwave signal and internally-generated parametric pumps at multiples of the Josephson frequency. This approach casts the dc-SQUID as a multi-mode coupling problem and, as such, is also amenable to a graph description. Here we describe the parametric coupling in the SQUID following the description in [19] by use of graphs to map the connections between modes. In particular, we will direct our attention to the interference and dissipation that must be present to generate nonreciprocity in the presence of gain.

A dc-SQUID consists of a pair of Josephson junctions inside a superconducting loop as in figure 5. Each junction is characterized by a phase difference ${{\phi }_{1,2}}$ across its terminals, but it is more convenient to recast the dynamics of the dc-SQUID in terms of the common and differential phase difference ${{\phi }_{C,D}}={{\phi }_{1}}\pm {{\phi }_{2}}$, where the subscript C (D) indicate common (differential) mode quantities respectively. The scattering parameters of the dc-SQUID can then be obtained by a suitable rotation of the mode basis as described in detail in [19]. ${{\phi }_{C,D}}$ can be expressed as the sum of fast oscillating terms at multiples of the Josephson frequency ${{\omega }_{J}}=2\pi {{V}_{{\rm dc}}}/{{\Phi }_{0}}$, where ${{\Phi }_{0}}$ is the magnetic flux quantum, and slow oscillating terms at the signal frequency ω. The fast oscillating terms act as parametric pumps causing mode conversion between ω and the mixing product frequencies $\omega \pm k{{\omega }_{J}}$. The lowest-dimensional model that can capture nonreciprocity involves two parametric pumps and a total of ten modes at five frequencies $\omega \pm 2{{\omega }_{J}},\omega \pm {{\omega }_{J}},\omega$. Each mode is doubly-degenerate, because for each frequency there are a common and a differential mode excitation. Some scattering subgraphs to describe the SQUID are shown in figure 6, where, for example, $\beta _{\omega ,\omega -{{\omega }_{J}}}^{{\rm c}.{\rm c}.}$ indicates the coupling rate between the common modes at frequency ω and $\omega -{{\omega }_{J}}$. The input and output modes are shown at the very bottom and top respectively. Calculation of the mode coupling rates as a function of the internally generated pump amplitudes is detailed in appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The number of subgraphs to describe the SQUID is too large for analytical calculations. However if the SQUID bias current is high enough, we can look for an approximate expression for the SQUID isolation for small values of the bias parameter $\epsilon ={{I}_{c}}/{{I}_{b}}$, where I_c is the critical current of the junctions and I_b is the bias current. The pump amplitudes and the coupling coefficients can then be calculated for small and we can simplify our calculations by considering only the subgraphs of leading order in . These subgraphs are shown in figure 6. The order of the graph is given by the product of its edges and therefore depends on the length of the path connecting the input to the output modes as well as the order of the edge weights (coupling rates). The order of the coupling rates increases with the difference between the mode drive frequencies (for example $\beta _{\omega ,\omega -k{{\omega }_{J}}}^{{\rm c}.{\rm c}.}\sim {{\epsilon }^{k}}$). Moreover a small ($\sim {{\epsilon }^{3}}$) direct coupling exists between modes at the same frequency, as explained in appendix B.

The dissipation rates for the common and differential modes when the SQUID is terminated into ideal infinite transmission lines, are given by

$$\begin{eqnarray}&&{{\gamma }_{C}}=\frac{1}{RC}=\frac{{{\omega }_{c}}}{{{\beta }_{c}}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{{\gamma }_{D}}=\frac{2R}{L}=\frac{2{{\omega }_{c}}}{\pi {{\beta }_{L}}},\end{eqnarray} \tag{ 27 }$$

with ${{\beta }_{L}}=(2\;L{{I}_{c}})/{{\Phi }_{0}}$, ${{\beta }_{c}}=2\pi {{I}_{c}}{{R}^{2}}\;C/{{\Phi }_{0}}$, ${{\omega }_{c}}=2\pi {{I}_{c}}R/{{\Phi }_{0}}$, L is the inductance of the SQUID loop, and R and C are the junction shunt resistance and capacitance. Since the SQUID is a nonresonant system, the mode frequencies ${{\omega }_{i}}$ in equation (3) are zero. We can approximately compute the SQUID isolation by use of the lowest-order subgraphs shown in figure 6. The first conclusion we can draw, based on our previous discussion, is that non-zero dissipation rates in (26) are essential to obtain nonreciprocity in the SQUID. A change in the impedance terminations at the various mode frequencies $\omega \pm k{{\omega }_{J}}$ will change the amount of nonreciprocity, viz., modes that are terminated in a short or an open circuit (corresponding to zero dissipation rates) will not contribute to the SQUID isolation. Moreover, it is important that ${{\gamma }_{C}}\ne {{\gamma }_{D}}$. In fact, by taking into account the symmetries between the mode couplings, the weight of the second order subgraphs in figure 6 is given by

$$\begin{eqnarray}&&{{w}_{{{\epsilon }^{2}}}}\propto 4\omega \operatorname{Re}\left[ \beta _{\omega +{{\omega }_{J}},\omega }^{{\rm c}.{\rm c}.}\beta _{\omega ,\omega +{{\omega }_{J}}}^{cd} \right]+4({{\gamma }_{C}}-{{\gamma }_{D}}){\rm Im}\left[ \beta _{\omega +{{\omega }_{J}},\omega }^{{\rm c}.{\rm c}.}\beta _{\omega ,\omega +{{\omega }_{J}}}^{cd} \right],\end{eqnarray} \tag{ 28 }$$

apart for a common factor equal to the product of the internal modes normalized detunings. If ${{\gamma }_{C}}={{\gamma }_{D}}=\gamma$, the second term in equation (28) is zero and the SQUID becomes reciprocal at first order. The same cancellation happens in the third-order subgraphs in figure 6 and in the fourth-order subgraphs, with the exception the subgraphs in the bottom-left corner of figure 6. In a practical device isolation can however be restored for ${{\gamma }_{C}}={{\gamma }_{D}}$ if the resonant coupling rate $\beta _{\omega ,\omega }^{cd}$ is increased. In a real device we expect such a rate to be higher, due to the stray capacitive coupling between the input and output as well as other parasitics [34]. Finally, in figure 7 we show the computed isolation and power gain of the SQUID as a function of the bias points for the exact numerical solution and for the lowest-order approximations. We find that graphs up to fourth order in are enough to accurately describe the isolation properties of the SQUID. The exact power gain was also computed numerically by extracting the impedance matrix from the scattering matrix of the SQUID and then calculating the gain as discussed in [19]. In principle the power gain may also be computed from the subgraphs, but the number of necessary graphs is large, even at lowest order. In this system, graphs help identify the elements of a complex multimode parametric scattering problem that contribute to the isolation properties. Finally, we observe that at high frequency the real part of the detuning ${{\Delta }_{C,D}}$ increases and therefore we expect the dc-SQUID isolation to decrease with frequency. This is consistent with the behavior discussed in [19, 34].

Zoom In Zoom Out Reset image size

In this work we introduced a graph theoretical representation that can be used as a combinatorial accounting tool to analyze nonreciprocity in coupled-mode systems and used it to derive a minimal unilateral parametric amplifier and analyze the conditions for nonreciprocity in SQUID amplifiers. An abstract graph associated to the coupled-mode system was used to compute the input/output scattering coefficients between modes and derive general constraints that need to be satisfied in order for nonreciprocity to occur. In order for a multi-mode system to be reciprocal the sum of the phases along any loop has to be 0 or π. If this condition is violated, phase and/or amplitude nonreciprocity is present and this condition can be interpreted as interference between different permutations of connecting paths. We also find that dissipation in the remaining disjoint vertices/modes is crucial for amplitude nonreciprocity/isolation. Specifically, it breaks the symmetry between forward and backward multi-mode scattering processes such that one may obtain constructive loop interference in one direction, and destructive interference in the reverse direction. As a result we were able to design a novel three-mode parametric amplifier, characterized by forward gain, input match and reverse isolation that can be integrated on a superconducting chip.

Although the scattering matrix of a multi-mode system can always be computed numerically using more conventional approaches, graphs offer a new perspective that can be useful in the problem of synthesizing new multi-mode systems and devices. For instance, this can be leveraged as a design tool, allowing one to reduce a set of desired device characteristics to a minimal description that can be utilized to design new amplifier/converter concepts such as the delta amplifier. As an approach, it is conceptually similar to the analysis of electric networks by means of circuit representations of the network components, where graphical manipulation rules are used to simplify the circuit analysis. We also showed, with the complex example of the dc-SQUID, that graphs aid in identifying the critical elements that create isolation in the presence of forward gain. Finally, graphs may provide a useful approach to engineering nonreciprocity in recent multi-mode parametrically coupled hybrid systems such as the electrical-mechanical-optical (three-mode) transduction bridge [35] and electro-mechanical amplifiers [30]. In general, problems like these are particularly well-suited to a graph-based approach since many of the interesting properties one is interested in (gain, isolation, gain-bandwidth product, etc) actually arise from the structure of the coupling network (the topology of the graph) and not the physical particulars of a given implementation.

The authors would like to thank Prof M Devoret (Yale) and Dr Archana Kamal (MIT) for motivating discussions.

In this section we formally derive equation (1). We begin by considering a set of N_r coupled resonators described by:

$$\begin{eqnarray}&&\frac{{\rm d}A}{{\rm d}t}=f(A,{{A}^{{\rm in}}}),\end{eqnarray} \tag{ A.1 }$$

where A is a vector of N_r internal complex normal mode amplitudes and Aⁱⁿ describes the set of input signal drives for each of these modes. We assume that A and Aⁱⁿ can be expressed as a sum of periodic functions. In this case harmonic balance can be used to rewrite equation (A.1) as a system of N_m coupled-mode equations:

$$\begin{eqnarray}&&\frac{{\rm d}B}{{\rm d}t}=-{\rm i}HB+K{{B}^{{\rm in}}},\end{eqnarray} \tag{ A.2 }$$

where ${{B}^{{\rm in}}}=[b_{1}^{{\rm in}},b_{p}^{{\rm in}},\ldots ,b_{p+1}^{{\rm in}\dagger },\ldots ,b_{{{N}_{m}}}^{{\rm in}\dagger }]$ is a vector of ${{N}_{m}}=p+q$ input stimuli and $B=[{{b}_{1}},{{b}_{p}},\ldots ,b_{p+1}^{\dagger },\ldots ,b_{{{N}_{m}}}^{\dagger }]$ is a vector of N_m internal modes that can be expressed as:

$$\begin{eqnarray}&&{{b}_{j}}={{\tilde{b}}_{j}}{{{\rm e}}^{-{\rm i}\omega _{j}^{s}t}},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&b_{j}^{{\rm in}}=\tilde{b}_{j}^{{\rm in}}{{{\rm e}}^{-{\rm i}\omega _{j}^{s}t}}for1\leqslant j\leqslant p,\end{eqnarray} \tag{ A.4 }$$

and the corresponding conjugate expressions for $p+1\leqslant j\leqslant {{N}_{m}}$. Here, we allow for the number of modes to exceed the number of resonators, ${{N}_{m}}\geqslant {{N}_{r}}$ and label the drive/response frequencies with an s, while the mode itself is labeled with a subscript j so that '$\omega _{j}^{s}$' denotes the drive/response frequency of mode j which can be based within any of the N_r resonators depending on how the coupling is directed within the function $f(a,{{a}^{{\rm in}}})$. The diagonal matrix $K={\rm diag}(\sqrt{\gamma _{1}^{{\rm ext}}},\ldots ,\sqrt{\gamma _{{{N}_{m}}}^{{\rm ext}}})$ describes coupling to the environment through the external dissipation rates $\{\gamma _{j}^{{\rm ext}}\}$. The diagonal elements h_jj of the ${{N}_{m}}\times {{N}_{m}}$ coupling matrix H can then be expressed in this mode basis as [36]:

$$\begin{eqnarray}{{h}_{jj}}=\left\{ \begin{array}{ccccccccccccccc} {{\omega }_{j}}-{\rm i}{{\gamma }_{j}}/2, & \quad {\rm for} & 1\leqslant j\leqslant p \\ -{{\omega }_{j}}-{\rm i}{{\gamma }_{j}}/2, & {} & p+1\leqslant j\leqslant {{N}_{m}}, \\ \end{array} \right.\end{eqnarray} \tag{ A.5 }$$

and the off-diagonal elements h_jk are given by

$$\begin{eqnarray}{{h}_{jk}}=\frac{1}{2}\left\{ \begin{array}{ccccccccccccccc} {{g}_{jk}}{{{\rm e}}^{{\rm i}(\omega _{j}^{s}-\omega _{k}^{s})t}}+{\rm c}.{\rm c}. & {\rm for} & 1\leqslant j,k\leqslant p \\ {{g}_{jk}}{{{\rm e}}^{{\rm i}(\omega _{j}^{s}+\omega _{k}^{s})t}}+{\rm c}.{\rm c}. & {} & 1\leqslant j\leqslant p, \\ {} & {} & \;p+1\leqslant k\leqslant {{N}_{m}} \\ -{{g}_{jk}}{{{\rm e}}^{{\rm i}(\omega _{j}^{s}+\omega _{k}^{s})t}}+{\rm c}.{\rm c}. & {} & 1\leqslant k\leqslant p, \\ {} & {} & \;p+1\leqslant j\leqslant {{N}_{m}} \\ -{{g}_{jk}}{{{\rm e}}^{{\rm i}(\omega _{j}^{s}-\omega _{k}^{s})t}}+{\rm c}.{\rm c}. & {} & p+1\leqslant j,k\leqslant {{N}_{m}}. \\ \end{array} \right.\end{eqnarray} \tag{ A.6 }$$

In equations (A.5) and (A.6) ${{\omega }_{j}}$ and ${{\gamma }_{j}}$ are the natural frequency and total dissipation rate for the resonator in which mode j resides, while the '*' indicates complex conjugation, and g_jk is the complex coupling coefficient between modes j and k. We further assume that ${{g}_{jk}}=g_{kj}^{*}$ if $1\leqslant j,k\leqslant p$ or $p+1\leqslant j,k\leqslant {{N}_{m}}$ and ${{g}_{jk}}=-g_{kj}^{*}$ for $1\leqslant j\leqslant p$ and $p+1\leqslant k\leqslant {{N}_{m}}$, or for $1\leqslant k\leqslant p$ and $p+1\leqslant j\leqslant {{N}_{m}}$. For example, the coupling matrix for a two-mode frequency converter or resonantly coupled oscillator system corresponds to the case $p=2,q=0$, while a conventional parametric amplifier couples two modes, one at a positive frequency (p = 1) and the other at a negative frequency (q = 1) [21]. With this general prescription for defining our mode basis, we can perform an input/output analysis of the system (equation A.2) and calculate the vector of output (scattered) fields b^out by means of the relation [36]:

$$\begin{eqnarray}&&b_{j}^{{\rm in}}-{{{\rm e}}^{{\rm i}{{\xi }_{j}}}}b_{j}^{{\rm out}}=\sqrt{\gamma _{j}^{{\rm ext}}}b.\end{eqnarray} \tag{ A.7 }$$

In general, the phase factor ${{\xi }_{j}}$ depends on the nature of the coupling to mode j. For simplicity, we set ${{{\rm e}}^{{\rm i}{{\xi }_{j}}}}=-1$, as one would obtain with a small series capacitance or weakly reflecting mirror. Other choices of coupling will have the effect of modifying the identity matrix on the rhs of equation (1). By substituting equations (A.3) and (A.4) into equation (A.2) and by use of the input/output boundary conditions (equation (A.7)), we can obtain the expression for the scattering matrix in equation (1).

A dc-SQUID consists of a pair of Josephson junctions inside a superconducting loop as in figure 5. Each junction is characterized by a phase difference ${{\phi }_{1,2}}$ across its terminals. Defining the common and differential phase difference ${{\phi }_{C,D}}={{\phi }_{1}}\pm {{\phi }_{2}}$, the dynamical equations describing the system are [19]:

$$\begin{eqnarray}&&\frac{{\rm d}{{\phi }_{c}}}{{\rm d}t}=\frac{{{\omega }_{B}}}{2}-{{\omega }_{c}}{\rm sin} ({{\phi }_{c}}){\rm cos} ({{\phi }_{d}}),\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}&&\frac{{\rm d}{{\phi }_{d}}}{{\rm d}t}=\frac{2R}{L}(\pi {{\phi }_{e}}-{{\phi }_{d}})-{{\omega }_{c}}{\rm cos} ({{\phi }_{c}}){\rm sin} ({{\phi }_{d}}),\end{eqnarray} \tag{ B.2 }$$

where ${{\omega }_{B}}=2\pi {{I}_{b}}R/{{\Phi }_{0}}$, and ${{\omega }_{c}}=2\pi {{I}_{c}}R/{{\Phi }_{0}}$. I_c is the junction critical current and I_b is the bias current of the SQUID, ${{\Phi }_{0}}$ is the magnetic flux quantum, L is the loop geometric inductance and R is the equivalent shunt resistance across the SQUID. The system (B.1) can be linearized around the following bias point:

$$\begin{eqnarray}&&{{\phi }_{c}}={{\omega }_{J}}t+{{\delta }_{c}},\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&{{\phi }_{d}}={{\phi }_{e}}+{{\delta }_{d}},\end{eqnarray} \tag{ B.4 }$$

where ${{\delta }_{c}},{{\delta }_{d}}$ are small perturbations. In [19] it was shown that by applying a harmonic balance procedure, the phases ${{\delta }_{c,d}}$ can be expressed as

$$\begin{eqnarray}&&{{\delta }_{c,d}}={{p}_{c,d}}+{{s}_{c,d}}=\mathop{\sum }\limits_{k=-M}^{M}p_{k}^{c,d}{{{\rm e}}^{{\rm i}k{{\omega }_{J}}t}}+s_{k}^{c,d}{{{\rm e}}^{{\rm i}(\omega +k{{\omega }_{J}})t}},\end{eqnarray} \tag{ B.5 }$$

where the components $p_{k}^{c,d}$ are internally generated pumps at frequency $k{{\omega }_{J}}=2\pi k{{V}_{{\rm dc}}}/{{\Phi }_{0}}$, where V_dc is the static dc voltage across the SQUID. The pumps cause parametric frequency mixing between the various signal modes $s_{k}^{c,d}$. In [19] the harmonic balance equations were solved for the case M = 3 to show how nonreciprocity arises from the interference between the different parametric conversion paths. As shown in the main text, we can employ the graph representation of a coupled-mode system to write explicitly the interference conditions needed to maximize the amount of amplitude nonreciprocity between two modes. In the M = 3 case we can describe the SQUID as a ten-mode system. The ten modes have frequencies $\omega \pm 2{{\omega }_{J}},\omega \pm {{\omega }_{J}},\omega$ and are doubly degenerate (for each frequency there are a common and a differential excitation). In order to calculate the scattering parameters, we need to find the coupling coefficients between the modes as a function of the internal pump amplitudes. This can be obtained by linearizing equations (B.1) for small signal mode-amplitudes and equating terms at the same frequency in order to obtain a set of coupled mode equation (A.2). After taking symmetries into account, the system can be described by the following five coupling coefficients

$$\begin{eqnarray}&&g_{\omega ,\omega -{{\omega }_{J}}}^{{\rm c}.{\rm c}.}={\rm i}\epsilon ({\rm cos} {{\phi }_{e}}-{\rm i}p_{2}^{c}{\rm cos} {{\phi }_{e}}-p_{2}^{d}{\rm sin} {{\phi }_{e}}),\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&g_{\omega ,\omega -{{\omega }_{J}}}^{cd}={\rm i}\epsilon ({\rm i}\;{\rm sin} \;{{\phi }_{e}}-p_{2}^{c}{\rm sin} {{\phi }_{e}}-{\rm i}p_{2}^{d}{\rm cos} \phi ),\end{eqnarray} \tag{ B.7 }$$

$$\begin{eqnarray}&&g_{\omega ,\omega -2{{\omega }_{J}}}^{{\rm c}.{\rm c}.}={\rm i}\epsilon ({\rm i}p_{1}^{c}{\rm cos} {{\phi }_{e}}-p_{1}^{d}{\rm sin} {{\phi }_{e}}),\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}&&g_{\omega ,\omega -2{{\omega }_{J}}}^{cd}={\rm i}\epsilon ({\rm i}p_{1}^{d}{\rm cos} {{\phi }_{e}}-p_{1}^{c}{\rm sin} {{\phi }_{e}}),\end{eqnarray} \tag{ B.9 }$$

$$\begin{eqnarray}&&g_{\omega ,\omega }^{cd}={\rm i}\epsilon ({\rm Im}(p_{1}^{d}){\rm cos} {{\phi }_{e}}-\operatorname{Re}(p_{1}^{c}){\rm sin} {{\phi }_{e}}),\end{eqnarray} \tag{ B.10 }$$

where for example $g_{\omega ,\omega -{{\omega }_{J}}}^{{\rm c}.{\rm c}.}$ corresponds to the coupling between the common modes at frequency ω and $\omega -{{\omega }_{J}}.$