# Algebraic canonical quantization of lumped superconducting networks

I. L. Egusquiza<EMAIL_ADDRESS>✉ Department of Physics, University of
the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain  EHU Quantum
Centre, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao,
Spain  A. Parra-Rodriguez<EMAIL_ADDRESS>✉ Department of
Physics, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao,
Spain  Institut Quantique and Département de Physique, Université de
Sherbrooke, Sherbrooke, Québec J1K 2R1, Canada

###### Abstract

We present a systematic canonical quantization procedure for lumped-element
superconducting networks by using a redundant configuration-space description.
The algorithm is based on an original, explicit, and constructive
implementation of the symplectic diagonalization of positive semidefinite
Hamiltonian matrices, a particular instance of Williamson’s theorem. With it,
we derive canonically quantized discrete-variable descriptions of passive
causal systems. We exemplify the algorithm with representative singular
electrical networks, a nonreciprocal extension for the black-box quantization
method, as well as an archetypal Landau quantization problem.

##  I Introduction

Superconducting circuits [1] are the current leading technology for quantum
information processing [2, 3, 4]. Complex as their quantum dynamics may be,
simple effective models based on microwave circuit theory have been extremely
succesful in explaining collective quantum phenemena therein. Indeed, the
3-dimensional (3D) Maxwell-London theory together with the Josephson equations
[5] are the main phenomenological framework routinely used to describe them
[6, 7]. However, in certain cases where characteristic lengths are smaller
than the wavelengths involved, e.g., when energy is localised at very small
volumes, such 3D equations can be further simplified to Kirchhoff’s laws,
which are representable through ideal lumped-elements circuits. In fact, the
canonical quantization is that of the differential equations governing the
dynamics of the collective degrees of freedom of those elements (like
capacitors, inductors and Josephson junctions) [8, 9, 10].

Nevertheless, bigger electromagnetic environments, such as superconducting 3D
cavities, are still commonly used in experiments to achieve, for instance,
long coherence times [11]. In order to incorporate these, and more general
linear passive systems with well-defined probing ports in Kirchhoff’s
framework, a two tier process has been suggested. First, a classical
multiphysics analysis of their linear response is performed, and then, its
output is fed to classical engineering theorems for fraction expansions of
causal passive matrices [12]. In this way, it is possible to obtain lumped-
element networks of the low-lying energy spectrum to the required level of
accuracy which are fit for quantization, a routine better known as the black-
box approach [13, 14, 15]. The classical expansions can definitely include
effective breaking of time-reversal symmetry [16] through ideal gyrator
(2-port) or circulator (N-port) devices [17, 18, 19, 20, 21, 22, 23]. These
are nonreciprocal (NR) elements that can be used for noise isolation or signal
routing [24], and that implement flux-charge constraints. Necessarily, they
have to be studied in the quantum regime as well [25, 26]. Furthermore, in
order to represent general linear responses [12] the introduction of Belevitch
transformers is required [27], yet another type of multi-terminal element,
which provide direct constraints between fluxes, and between charges.

Thus, it has become of great interest to find a systematic quantum circuit
analysis [8, 28, 29, 30, 9, 31, 32, 33, 34, 13, 14, 15, 35, 36, 37, 10, 38,
39, 40, 25, 41, 42, 43, 44, 45, 43, 46, 47] of electrical networks that treats
all linear passive lumped elements (capacitors, inductors, ideal nonreciprocal
elements and Belevitch transformers) on an equal footing. This analysis must
deal with obstacles such as singular kinetic matrices that would prevent
straightforward Legendre transformations. We stress a Lagrangian formulation
since it gives us both canonical Poisson brackets, later to be promoted to
commutators, and Hamiltonian functions on phase space, later to be promoted to
operators.

In this article, we present a systematic procedure for deriving Hamiltonians
together with canonical quantizations of lossless lumped-element networks.
This procedure is based on a Lagrangian written in terms of an enlarged
configuration space with both flux and charges variables [48, 35, 44]. Such
enlargement with respect to most of previous approaches (where Lagrangians
were written in either flux or charge variables) allows us to replace the
troublesome Legendre transformation for singular kinetic matrices by a simple
nonsingular matrix inversion. Naturally, there is a cost to pay with this
approach, as redundancies will be certainly introduced in the form of free
particles and nondynamical degrees of freedom. However, making use of the
full-fledged Williamson theorem for normal forms of quadratic systems [49],
these subsectors can be systematically disentagled from the nontrivial
dynamical sector and canonical quantization can follow effortlessly.

The underlying mathematical keystone of the present work is a simple and
constructive algorithm to find systematically the symplectic normal forms of
positive semi-definite symmetric (real) matrices (with eigenvalues
$\lambda_{i}\geq 0$) inspired by Hörmander’s analysis [50] of Williamson’s
seminal article [49]. We consider and propose the application of this routine
to find canonical quantum descriptions of general passive causal systems
presenting a discrete set of frequencies. We stress that this includes the
possibility of constraints that hitherto have been considered an obstacle to
achieve that goal. In particular, we exploit it to obtain systematically exact
quantum descriptions of superconducting lumped-element networks, including
cases analogous to the Landau problem, with discrete macroscopic conjugated
charges and fluxes replacing microscopic position variables.

The results obtained here generalize systematic derivations of canonically
quantized Hamiltonians for reciprocal [9, 10, 32, 33, 35, 51], and
nonreciprocal [25, 26] superconducting lumped-element networks based on graph
theory, and will provide support for more general quasi-lumped descriptions of
circuits [40, 43, 52]. We remark that a particular and important consequence
of the application of our algorithm is the generalization of the black-box
quantization methods of non-linear degrees of freedom coupled to nonreciprocal
linear systems [13, 14, 15].

The article is structured as follows: in Sec. II we introduce an explicit
constructive algorithm of Williamson’s theorem to find symplectic
transformations that diagonalize real positive-semidefinite matrices of
dimension $2n$. For illustration in a simple understandable case, we apply it
to a Landau quantization problem with two linearly-coupled charged particles
in a homogeneous magnetic field. We make full use of the main theorem in Sec.
III to put forward a systematic process to find canonically quantized
Hamiltonians of lumped-element superconducting networks described in terms of
Kirchhoff’s equations, starting from a Lagrangian written in the enlarged
configuration space. In Sec. IV, we implement the algorithm to quantize a
singular nonreciprocal circuit containing both nondynamical degrees of freedom
and free particles using both (i) a black-box approach and (ii) keeping the
full cosine potential of the nonlinear degrees of freedom. We finally draw
some conclusions and describe possible additional applications in Sec. V.

##  II Williamson’s Theorem for Positive Semi-definite Symmetric Matrices

The study of the classification, with canonical presentations, of quadratic
forms has a long history. For instance, Weierstrass ([53] in p. 233 of [54],
and [55] published with some changes in p. 19 of [56]) looked into the
conditions for pairs of quadratic forms to be put in canonical form (Gestalt).
Part of the motivation of Weierstrass was the application of these canonical
forms to the theory of small oscillations, and this reason and others kept the
topic alive. In 1936, Williamson gave the complete classification of quadratic
forms under symplectic transformations [49], over general fields of
characteristic zero. His classical work did not provide us with explicit
constructive methods nor did it include a complete catalogue of canonical
forms. Laub and Meyer [57] did provide a listing of canonical forms, which was
refined and completed by Hörmander [50]. In fact, many other authors have
looked into the issue without being aware of each other [58, 59]. An
alternative source of information for physicists is the compilation of
symplectic normal forms for quadratic Hamiltonians presented by Arnol’d in
Appendix 6 of [60] and attributed by him to D. M. Galin, without further
reference.

In a more recent development, the study of quantum information with continuous
variables, and in particular of Gaussian states [61] has brought a special
case of Williamson’s theorem to the fore, that of positive definite quadratic
matrices. Of special relevance in this context is the symplectic
diagonalization of covariance matrices [62]. Some attention to the topic has
also been given in the context of superconducting circuits and other quantum
systems. In these applications a constructive approach is required, beyond the
descriptions in the theorems, and much has been made of this fact for the
simplest case of positive definite Hamiltonians [63], while that of positive
semidefinite Hamiltonians has been referred to but not addressed. Here we
remedy that deficiency.

Let us denote as $\mathsf{H}$ the matrix corresponding to a quadratic
Hamiltonian
$H=\operatorname{\boldsymbol{x}}^{T}\mathsf{H}\operatorname{\boldsymbol{x}}/2$,
while $\mathsf{J}$ is the canonical symplectic matrix, such that the
(classical and Heisenberg) dynamics reads

$\displaystyle\dot{\operatorname{\boldsymbol{x}}}=\operatorname{\mathsf{J}}\operatorname{\mathsf{H}}\operatorname{\boldsymbol{x}}$
(1)

for phase space coordinates
$\operatorname{\boldsymbol{x}}=(q_{1},...,q_{n},p_{1},...,p_{n})^{T}=(x_{1},...,x_{2n})^{T}$.
The content of Williamson’s theorem for positive semidefinite Hamiltonians is
that there are only three possiblities for independent degrees of freedom
present in Eq. (1) when $\operatorname{\mathsf{H}}$ belongs to that class
[50]: (i) harmonic oscillators, evolving as
$\ddot{\operatorname{\boldsymbol{y}}}_{\mathrm{ho}}=\operatorname{\mathsf{\Omega}}^{2}\operatorname{\boldsymbol{y}}_{\mathrm{ho}}$,
where $\operatorname{\mathsf{\Omega}}$ is a positive definite diagonal matrix
of frequencies, (ii) free particles with
$\ddot{\operatorname{\boldsymbol{y}}}_{\mathrm{fp}}=0$, and (iii) nondynamical
variables, i.e.,
$\operatorname{\boldsymbol{y}}_{\mathrm{nd}}(t)=\operatorname{\boldsymbol{y}}_{\mathrm{nd}}(0)$.
These variables $\operatorname{\boldsymbol{y}}$ are the desired objective,
related to the initial ones by a symplectic transformation
$\operatorname{\boldsymbol{x}}=\operatorname{\mathsf{S}}\operatorname{\boldsymbol{y}}$.

To proceed, we have to identify the linear subspaces
$K_{1}=\mathrm{ker}\left[\mathsf{JH}\right]=\mathrm{ker}\left[\mathsf{H}\right]$
and $K_{2}=\mathrm{ker}\left[\left(\mathsf{JH}\right)^{2}\right]$. If
$K_{1}=K_{2}$, there are no free particles, and all the variables
corresponding to $K_{1}$ are nondynamical. More generally,
$n_{\mathrm{f}}=\mathrm{dim}(K_{2})-\mathrm{dim}(K_{1})$ is the number of free
particles and $\mathrm{dim}(K_{1})-n_{f}=2n_{\mathrm{nd}}$ is even and twice
the number of nondynamical degrees of freedom $n_{\mathrm{nd}}$. Notice that
$\mathrm{dim}(K_{2})=2n_{\mathrm{f}}+2n_{\mathrm{nd}}$ is necessarily even.
There are $n_{\mathrm{ho}}=n-n_{\mathrm{f}}-n_{\mathrm{nd}}$ harmonic
oscillators, to complete the description, with $2n$ the dimension of phase
space.

The classification above translates into the existence of a change of basis
matrix $\mathsf{F}$ for which

$\displaystyle\mathsf{F}^{-1}\mathsf{JH}\mathsf{F}=\operatorname{\mathsf{D}}_{\operatorname{\mathsf{J}}}=\begin{pmatrix}0&0&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&0&0&\text{1}&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&i\operatorname{\mathsf{\Omega}}&0\\\
0&0&0&0&0&-i\operatorname{\mathsf{\Omega}}\end{pmatrix}.$ (2)

The order has been chosen in the sequence nondynamical/free particles/harmonic
oscillators for the degrees of freedom. The dimensions, not explicit in the
expression above, run in the sequence
$\left\\{n_{\mathrm{nd}},n_{\mathrm{nd}},n_{\mathrm{f}},n_{\mathrm{f}},n_{\mathrm{ho}},n_{\mathrm{ho}},\right\\}$
for the diagonal boxes. This matrix $\mathsf{F}$ is the output of the usual
Jordan canonical form process for $\mathsf{JH}$, up to ordering.

Generically, this change of basis matrix does not diagonalize $\mathsf{H}$ by
a symplectic transformation. For starters, it is necessarily complex. However,
an adequate transformation will be of the form

$\mathsf{S}=\operatorname{\mathsf{F}}\operatorname{\mathsf{\Lambda}}\operatorname{\mathsf{F}}_{0}^{-1}\,,$
(3)

where $\mathsf{F}_{0}$ is the matrix

$\mathsf{F}_{0}=\begin{pmatrix}\text{1}&0&0&0&0&0\\\ 0&\text{1}&0&0&0&0\\\
0&0&\text{1}&0&0&0\\\ 0&0&0&\text{1}&0&0\\\ 0&0&0&0&\text{1}&\text{1}\\\
0&0&0&0&i\text{1}&-i\text{1}\end{pmatrix}\,.$ (4)

As we show explicitly in Appendix A, the structure of $\mathsf{F}_{0}$ comes
about because the real and imaginary parts of the eigenvectors of
$\mathsf{JH}$ with purely imaginary eigenvalues can be arranged in symplectic
pairs, while the generalized eigenvectors of the zero eigenvalue are real.

To complete the description, $\mathsf{\Lambda}$ is a block diagonal matrix,
whose detailed structure depends on the existence or not of degeneracies both
in the harmonic oscillator sector and in the generalized eigenspace of
$\mathsf{JH}$ for eigenvalue 0. Let us first examine its structure in the
harmonic oscillator sector. If there are no nonzero degenerate frequencies,
the corresponding block of $\mathsf{\Lambda}$ consists of two diagonal
matrices of dimension $n_{\mathrm{ho}}$ for the harmonic oscillator sector,
with the diagonal elements being normalization factors. This is in fact the
case that has been frequently explored in the literature, since it pertains to
the positive definite case without degeneracies, the only one that has merited
full attention, as it is the easiest. Those diagonal elements in this case
simply give the normalization of each eigenvector of $\mathsf{JH}$ for purely
imaginary eigenvalue, such that its real and imaginary part are symplectic
orthonormal.

Were there a degeneracy in the nonzero spectrum of $\mathsf{JH}$, the
corresponding block for $\mathsf{\Lambda}$ will be due to symplectic
Gram–Schmidt orthonormalization of the real and imaginary parts of the
eigenvectors for that eigenvalue. Similarly, the $\mathsf{\Lambda}$ is
required, in the free particle and non-dynamical sector, to ensure that the
free particle sector is symplectic orthogonal to the non-dynamical one, and
that the resultant free particle sector is symplectic orthonormal. This can be
done explicitly, as shown in Appendix A, or through the condition that

$\mathsf{J}=\mathsf{S}^{T}\begin{pmatrix}\omega_{0}&0\\\
0&\mathsf{J}_{\mathrm{ho}}\end{pmatrix}\mathsf{S}\,,$ (5)

where $\omega_{0}$ is an antisymmetric
$(2n_{\mathrm{nd}}+2n_{\mathrm{f}})\times(2n_{\mathrm{nd}}+2n_{\mathrm{f}})$
while $\mathsf{J}_{\mathrm{ho}}$ is the canonical symplectic matrix, now for
the harmonic oscillator sector. If there are degeneracies, the harmonic
oscillator sector of $\mathsf{\Lambda}$ becomes block diagonal, with pairs of
invertible matrices for each degeneracy subspace. Complete details are
presented in Appendix A. Once $\mathsf{\Lambda}$ has been determined, we
compute $\mathsf{S}^{T}\mathsf{H}\mathsf{S}$, resulting in a block diagonal
matrix that is in fact canonically diagonal in the harmonic oscillator sector,
i.e. $\mathrm{diag}\left(\mathsf{\Omega},\mathsf{\Omega}\right)$. If one
desires full symplectic diagonalization an additional symplectic Gram–Schmidt
process will be required for $K_{2}$, as made explicit in Appendix A.

###  II.1 A familiar example: the Landau problem

The Landau problem of one charged particle moving in a homogeneous magnetic
field is very well known indeed. Simple as its statement is, its correct
solution requires a canonical transformation that involves both positions and
momenta. Furthermore, it is the starting point for our understanding of the
quantum Hall effect. We consider a slight extension of the Landau problem, and
address it with the symplectic diagonalization algorithm so as to illustrate
its applicability in a simple example. Let us thus consider a pair of
identical linearly-coupled charged particles subjected to a constant and
homogeneous external magnetic field in the $z$-direction
($\boldsymbol{B}=(0,0,B)$), see Fig.1.

Figure 1: Two coupled charged particles subjected to the action of a magnetic
field. The dynamics of this system naturally separate in two orthogonal
sectors, the in-plane ($x$–$y$) and normal ($z$) ones.

Working in the symmetric gauge, i.e.,
$\boldsymbol{A}(\operatorname{\boldsymbol{q}}_{i})=-\frac{1}{2}\left(\boldsymbol{B}\times\operatorname{\boldsymbol{q}}_{i}\right)$,
a redundant minimal-coupling Hamiltonian can be written as

$\displaystyle H=$
$\displaystyle\sum_{i}\frac{\left(\operatorname{\boldsymbol{p}}_{i}-e\boldsymbol{A}(\operatorname{\boldsymbol{q}}_{i})\right)^{2}}{2m}+\frac{k}{2}(\operatorname{\boldsymbol{q}}_{1}-\operatorname{\boldsymbol{q}}_{2})^{2},$
(6)

where the dynamics of the phase-space separates in two independent subsectors,
$2H=\operatorname{\boldsymbol{X}}_{z}^{T}\operatorname{\mathsf{H}}_{z}\operatorname{\boldsymbol{X}}_{z}+\operatorname{\boldsymbol{X}}_{xy}^{T}\operatorname{\mathsf{H}}_{xy}\operatorname{\boldsymbol{X}}_{xy}$,
i.e. vertical and in-plane motion. Williamson’s analysis of
$\operatorname{\mathsf{H}}_{z}$ yields, as symplectic diagonal degrees of
freedom, one harmonic oscillator and one free particle. In the limit of the
coupling constant going to zero ($k\rightarrow 0$), obviously, the harmonic
oscillator frequency tends to zero and in fact we have two free particles.

The analogous analysis in the $x$–$y$ plane shows the existance of three
harmonic conjugate pairs, and an extra non-dynamical one. Again, by taking the
zero-coupling limit, the dynamics in the plane split in two pairs of harmonic
and nondynamical conjugated variables. See Appendix B for further details and
an explicit expression for a symplectic diagonalizing matrix.

Observe that the procedure by itself does not prescribe a meaning for the
nondynamical sector. Thus, in the case of the mechanical system with just
these two distinguishable particles, the nondynamical sector is set to being
identically zero, and would not be detectable if the system were realized. Not
so if we take this study as part of second quantization for fermions, say, in
which case we would understand fermions in the same dynamical state to be
allowed to be in different states in the nondynamical sector, as is done in
the analysis of the quantum Hall effect. Nonetheless, our procedure allows us
the systematic identification of the different sectors.

In this particular case, one might wonder whether our analysis depends in any
way on having chosen the symmetric gauge. However, one crucial point to be
stressed in the Williamson diagonalization procedure is that the canonical
form under symplectic transformations is a symplectic invariant, by
construction. And since all the gauge transformations that preserve the
quadratic character of the Hamiltonian are actually implementable as
symplectic transformations, the physical content of the model is seen as
clearly gauge invariant, and achieved systematically. The gauge invariance of
the problem has led to some confusion in the literature as to the conserved
quantities associated with the model. Again, looking at the issue from the
point of view of symplectic invariance allows the immediate identification of
the maximum number of linear and quadratic conserved quantities for these
families of systems. Their physical interpretation is on the other hand
directly dependent on how these systems are coupled to others, and are not
intrinsic to the model. For more details on symmetries and gauge invariance in
this context, see Appendices C and C.1, respectively.

##  III Lumped-element electrical networks

In order to derive a systematic, canonically quantized theory of
superconducting lumped-element circuits, it is useful to take as an starting
point a classical Lagrangian [9, 10]. The additivity of components in the
Lagrangian is immediately obvious, while additivity for the Hamiltonian is not
necessarily as straightforwrd. Furthermore, a nonsingular Lagrangian
determines the Poisson brackets, and one does not need to inferr them to
correctly match Kirchhoff’s equations of motion [64]. One of the problems in
the analysis of circuitry, therefore, consists on the construction of
nonsingular Lagrangians in a systematic way. This objective, however, has
hitherto compelled us to restrict the type of circuit to be analyzed so that
there is a well defined mechanism for the identification of independent
variables that still does provide us with a nonsingular Lagrangian.

In contrast, we propose to construct directly a nonsingular Lagrangian and
eliminate redundancies by the systematic application of Williamson’s
symplectic diagonalization. Thus we consider a quite possibly redundant
configuration space, such that all linear energy bearing lumped elements are
assigned to a kinetic energy. In this way we are sure of the applicability of
the Legendre transform. We shall see that there are two sources of redundancy,
which, for the linear systems that we consider, correspond to the nondynamical
and free particle degrees of freedom mentioned above.

This path followed here is inspired by the one taken to describe transmission
lines coupled by ideal NR elements in [44], where a doubled configuration
space description of the lines was introduced to completely determine the
normal mode structure of the combined system. As stated, the kinetic energy
will be the addition of the contributions of all lumped inductors and
capacitors in the network. This will be the only set of linear energetic term
of the lagrangian, which will be completed by linear nonenergetic constraint
terms and relevant nonlinear energetic terms. Naturally a reduction of
variables in the configuration space is useful to simplify the algorithm, but
we are going to work on worst-case (brute force) scenarios. For the sake of
completeness, we present the Lagrangian terms describing the nonreciprocal
ideal elements [25, 44, 26, 65], the Belevitch transformer [15], the nonlinear
elements like the Josephson and the phase-slip junction [66, 67] and discuss
their effect on the Hamiltonian derivation.

###  III.1 Linear reciprocal networks

Historically, the analysis of lumped element networks has been performed using
the language of graph-theory [9, 32, 33]. Let us consider a common linear
reciprocal network of 2-terminal lumped elements consisting on a set of
inductors and capacitors connected in a tree containing $N$ number of nodes
and $M$ number of loops such that $M\leq N$, as in the example in Fig. 2.

Figure 2: LC network. Asigning loop-charge (branch-flux) variables for the
inductors (capacitors) it is possible to write systematically a nonsingular
kinetic Lagrangian term in a double space.

Currents and voltages are denoted by time derivatives of loop charges
$\dot{Q}\equiv\partial_{t}Q$ and node fluxes $\dot{\Phi}$ respectively.
However, instead of using node-flux variables, let us assign a branch-flux
variable for every capacitor ($N_{C}$), and a loop-charge variable for every
inductor ($M_{L}<M$). The equations of motion of such a generic network can be
written in the form

$\displaystyle\operatorname{\mathsf{C}}\ddot{\operatorname{\boldsymbol{\Phi}}}+\operatorname{\mathsf{D}}^{T}\dot{\operatorname{\boldsymbol{Q}}}$
$\displaystyle=0,$ (7)
$\displaystyle\operatorname{\mathsf{L}}\ddot{\operatorname{\boldsymbol{Q}}}-\operatorname{\mathsf{D}}\dot{\operatorname{\boldsymbol{\Phi}}}$
$\displaystyle=0.$ (8)

where $\operatorname{\mathsf{C}}$ and $\operatorname{\mathsf{L}}$ are positive
diagonal matrices (therefore both full rank), and $\operatorname{\mathsf{D}}$
is an adjacency matrix encoding the topology of the terminal connections. This
set of equations are a generalized discrete limit of the telegrapher’s
equations in a transmission line [68, 44] and flow from the following
Lagrangian, written in the doubled (doubled as both fluxes and charges appear)
configuration space description as

$L_{\mathrm{LC}}=\frac{1}{2}\dot{\operatorname{\boldsymbol{\Phi}}}^{T}\operatorname{\mathsf{C}}\dot{\operatorname{\boldsymbol{\Phi}}}+\frac{1}{2}\dot{\operatorname{\boldsymbol{Q}}}^{T}\operatorname{\mathsf{L}}\dot{\operatorname{\boldsymbol{Q}}}-\dot{\operatorname{\boldsymbol{Q}}}^{T}\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}.$
(9)

The total apparent number of degrees of freedom is $N_{C}+M_{L}$, but there
are constraints given by $\operatorname{\mathsf{D}}$ and
$\operatorname{\mathsf{D}}^{T}$. In fact, it is well known that the number of
node-flux and of loop-charge variables are upper bounds to the number of real
harmonic degrees of freedom of a network [35]. In contrast to previous works
on canonical quantization of electric networks, let us assume for the sake of
generality that we describe the Hamiltonian with redundancy, i.e., with a
configuration space of $N_{C}+M_{L}$ variables.

The invertibility of the kinetic matrices $\operatorname{\mathsf{C}}$ and
$\operatorname{\mathsf{L}}$ permits a Legendre transformation to obtain the
Hamiltonian

$H_{\mathrm{LC}}=\frac{1}{2}\operatorname{\boldsymbol{\Pi}}^{T}\operatorname{\mathsf{C}}^{-1}\operatorname{\boldsymbol{\Pi}}+\frac{1}{2}\left(\operatorname{\boldsymbol{P}}-\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}\right)^{T}\operatorname{\mathsf{L}}^{-1}\left(\operatorname{\boldsymbol{P}}-\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}\right),$
(10)

where the conjugate variables are $\operatorname{\boldsymbol{\Pi}}=\partial
L/\partial\dot{\operatorname{\boldsymbol{\Phi}}}$ and
$\operatorname{\boldsymbol{P}}=\partial
L/\partial\dot{\operatorname{\boldsymbol{Q}}}$. The dimension of the phase
space is thus $2\times(N_{C}+M_{L})$. This Hamiltonian is, on its face,
positive semidefinite, as $\mathsf{C}^{-1}$ and $\mathsf{L}^{-1}$ are positive
definite. The zero modes of the Hamiltonian will necessarily be associated
with the possibility that
$\operatorname{\boldsymbol{P}}-\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}$
be zero. By Williamson’s theorem, there exists a canonical transformation that
expresses the Hamiltonian as describing a set of free particles, harmonic
oscillators, and nondynamical conjugate pairs of variables, see Appendix D.
Studying the kernels of $\operatorname{\mathsf{JH}}$ and
$(\operatorname{\mathsf{JH}})^{2}$, it can be easily proven that number of
nondynamical pairs is the same as that of harmonic oscillators
$n_{nd}=n_{ho}=[N_{C}+M_{L}-n_{f}]/2$. The free particles are due to ignorable
coordinates in the Hamiltonian, i.e. due to some translation symmetry. Those
arise from the reducibility of the initial circuit by series/parallel
equivalences. The nondynamical sector, on the other hand, turns out to be
there because we have indeed doubled the number of variables, once the full
reduction in terms of those equivalences has been carried out.

Let us illustrate the idea with a trivial example of a harmonic LC-oscillator
with unnecessarily redundant descriptions.

Figure 3: (a) Equivalent representations of a harmonic oscillator with: (a) an
LC circuit; its minimal representation, (b) two capacitors in series and (c) a
gyrator and a capacitor replacing the inductor. Redundant representations (b,
c) may contain free-particle and/or nondynamical degrees of freedom.

The circuits in Fig. 3 are three different circuital representations of an LC
harmonic oscillator. It is worth realizing that even its simplest version (a)
admits irreducible representations in terms of only a flux $(\Phi)$ or charge
coordinate $(Q)$, i.e., the equations of motion can be written in terms of a
single variable, without requiring a redundant description with both variables
in the configuration space. Even more, the circuit in (b) has two capacitors
in series such that the new set of flux coordinates $(\Phi_{1},\Phi_{2})$ adds
an extra level of redudancy. Indeed, upon the use of the sum rule for series
capacitances $C_{s}=C_{1}\parallel C_{2}=(C_{1}+C_{2})/C_{1}C_{2}$, it is
possible to eliminate a flux coordinate to reach (a).

We now apply mechanically the method described above to circuit (b). That
means assigning a variable and a kinetic energy term to each capacitor and
inductor of the system, and computing the adjacency matrix to determine the
constraint terms. Thus Lagrangian (9) should now be used, with matrices

$\displaystyle\operatorname{\mathsf{C}}=\begin{pmatrix}C_{1}&0\\\
0&C_{2}\end{pmatrix},\quad\operatorname{\mathsf{L}}=\begin{pmatrix}L\end{pmatrix},\quad\operatorname{\mathsf{D}}=\begin{pmatrix}1&1\end{pmatrix},$
(11)

with vector variabless
$\operatorname{\boldsymbol{\Phi}}=(\Phi_{1},\Phi_{2})^{T}$ and
$\operatorname{\boldsymbol{Q}}=(Q)$. The Euler-Lagrange equations (7,8) reveal
the big redundancy in the description. Traditionally, one would eliminate the
charge $Q$ by integrating Eq. (8) in time while setting to zero the constant
initial charge, i.e.,
$\dot{Q}=L^{-1}\operatorname{\mathsf{D}}{\operatorname{\boldsymbol{\Phi}}}$,
in which case, eq. (7) simplifies to
$\operatorname{\mathsf{C}}\ddot{\operatorname{\boldsymbol{\Phi}}}=\operatorname{\mathsf{D}}^{T}\operatorname{\mathsf{L}}^{-1}\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}$.
Given these equations, one realizes that a simpler Lagrangian exists, and,
after a contact transformation of coordinates,
$\operatorname{\mathbf{\Psi}}=\mathsf{E}\operatorname{\boldsymbol{\Phi}}$, the
harmonic oscillator nature of the system is revealed with a new Lagrangian
$L=\frac{1}{2}(C_{s}\dot{\Psi}^{2}-L^{-1}\Psi^{2})$, where
$C_{s}=(C_{1}+C_{2})/C_{1}C_{2}$.

However, in our proposal we do not need to perform at this stage the circuital
reduction, and in fact it will come out the systematics. Namely, we write
Hamiltonian (10), and apply Williamson’s theorem. That is, we find a
symplectic-orthonormal basis such the Hamiltonian matrix in this basis is
diagonal. This computation reveals the presence of a nondynamical variable and
a free particle in addition to the harmonic oscillator. The free particle is
directly related to the reducibility to an equivalent capacitor of the two
series capacitors, while the nondynamical variable is a consequence of our
redundant description. See Appendix D for further details.

###  III.2 Nonreciprocal linear ideal devices

Having understood the analysis for the case of linear lumped energetic
elements, let us now discuss the introduction of the nonreciprocal ideal ones,
which do not store energy. In particular we address those breaking time-
reversal symmetry. These multi-terminal devices represent _constraints_
between fluxes and charges, see Fig. 4(a,b) as expressed in their constitutive
equation

$(1-\operatorname{\mathsf{S}})\dot{\bar{{\operatorname{\boldsymbol{\Phi}}}}}=R(1+\operatorname{\mathsf{S}})\dot{\bar{{\operatorname{\boldsymbol{Q}}}}},$
(12)

written in terms of a nonsymmetric scattering matrix
$\operatorname{\mathsf{S}}$($\neq\operatorname{\mathsf{S}}^{T}$), which always
exists [69]. The overbar for the fluxes and charges denotes that they have
been rescaled dimensionally, for ease of notation in the constitutive
equation. In fact, such elements are more easily introduced in the Lagrangian
formalism in terms of their impedance and/or admittance descriptions, namely

$\displaystyle
L_{\operatorname{\mathsf{Y}}}=\frac{1}{2}\dot{\bar{{\operatorname{\boldsymbol{\Phi}}}}}^{T}\operatorname{\mathsf{Y}}\bar{\operatorname{\boldsymbol{\Phi}}},\quad\text{and}\quad
L_{\operatorname{\mathsf{Z}}}=\frac{1}{2}\dot{\bar{{\operatorname{\boldsymbol{Q}}}}}^{T}\operatorname{\mathsf{Z}}\bar{\operatorname{\boldsymbol{Q}}},$
(13)

with ports described by branch-flux or loop-charge coordinates respectively.
Observe that in these descriptions only the antisymmetric part of either
admittance or impedance can have an effect on the equations of motion, since
the symmetric part gives a total time derivative. Further, notice that there
is no energy associated with these magnetic terms.

Figure 4: Multiport lossless linear system. The (a) gyrator and (b) circulator
are 2- and 3-port ideal nonreciprocal elements implementing constraints
between flux and charge variables and break of time-reversal symmetry
($\mathcal{T}$). A (c) right-hand Belevitch transformer yields flux-flux and
charge-charge constraints.

The general nonreciprocal constitutive equation described by
$\operatorname{\mathsf{S}}$ is more directly read as a constraint on the
evolution of the flux and charge variables, albeit nonholonomic. This
nonholonomic character suggests that it might best be treated with tools other
than Lagrangians on a flux coordinate space. However, Eq. (12) can always be
framed in Lagrangian terms as one of Eqs. (13) as long as one accepts the
possibility of immittance matrix being singular, in which case one further
requires an additional flux or charge constraint (which can enter as another
Lagrange multiplier) at the ports [25, 70], i.e.

$\displaystyle\operatorname{\mathsf{P}}_{-}\dot{\operatorname{\boldsymbol{\Phi}}}$
$\displaystyle=$ $\displaystyle 0,\quad\text{if
}-1\in\lambda_{\operatorname{\mathsf{S}}},$ (14)
$\displaystyle\operatorname{\mathsf{P}}_{+}\dot{\operatorname{\boldsymbol{Q}}}$
$\displaystyle=$ $\displaystyle 0,\quad\text{if
}\quad\,1\in\lambda_{\operatorname{\mathsf{S}}}.$ (15)

In either case, starting from the doubled configuration space construction, a
Legendre transformation will be available after this reduction is properly
dealt with. We remark the equivalence of the coupling terms in (13) with the
$\operatorname{\mathsf{D}}$-term in Lagrangian (9). In fact, there is a
mapping between both expressions by rescaling some of the flux (charge)
variables by $R^{-1}$ ($R$), i.e. $\tilde{Q}_{i}=\Phi_{i}/R$. This is the
reason for the enlarged configuration space approach, involving both fluxes
and charges, to be natural in dealing with circuits with nonreciprocal
elements. In other words, one can think of an inductor (capacitor) in a
circuit always as an equivalent capacitor (inductor) when seen through a
gyrator, see Fig. 2(c), in which one of the capacitors can be moved (in the
sense of equivalence of circuits) through the gyrator, becoming thus an
inductor, providing us with a read on the circuit as an LC oscillator.

###  III.3 Belevitch transformers

The final ideal passive linear element to be considered is the Belevitch
transformer [27], another even (multi) $2N$-terminal device which represents
direct constraints between voltages and currents on $N_{L}$ pairs of terminals
on the left side, and $N_{R}$ pairs on the other ($N=N_{L}+N_{R}$), see Fig.
4(c). As such, they can be written as Lagrange multipliers, and can be
directly eliminated before the Legendre transformation in the configuration,
see further in [15, 40]. They can be classified in right (R), and left (L)
transformers, such that their constraints are expressed as the Lagrangian
terms

$\displaystyle
L_{\mathrm{ET}}=\mathsf{\Lambda}_{E}\dot{\operatorname{\boldsymbol{Q}}}+\mathsf{\Xi}_{E}\dot{\operatorname{\boldsymbol{\Phi}}},$
(16)

where $E=\\{R,L\\}$, and

$\displaystyle\mathsf{\Lambda}_{R}\dot{\operatorname{\boldsymbol{Q}}}$
$\displaystyle=$
$\displaystyle\mathsf{\Lambda}(\dot{\operatorname{\boldsymbol{Q}}}_{R}+\operatorname{\mathsf{N}}\dot{\operatorname{\boldsymbol{Q}}}_{L}),$
(17) $\displaystyle\mathsf{\Xi}_{R}\dot{\operatorname{\boldsymbol{\Phi}}}$
$\displaystyle=$
$\displaystyle\mathsf{\Xi}(\dot{\operatorname{\boldsymbol{\Phi}}}_{L}-\operatorname{\mathsf{N}}^{T}\dot{\operatorname{\boldsymbol{\Phi}}}_{R}),$
(18)

and analogously for $L_{\mathrm{LT}}$, where $\operatorname{\mathsf{N}}$ is
the turn-ratios (rectangular) matrix describing the transformer. Such elements
are known to be required for the systematic analysis and synthesis of generic
linear electrical systems. For instance, in Eq. (8) we have assumed the
inductance matrix to be full rank. However, this may not be the case generally
(say when tighly-coupled mutual inductances appear, e.g.,
$M_{ij}=\sqrt{L_{i}L_{j}}$). In such situations, the system could, among other
possibilities, be written in terms of transformers and a full rank inductance
matrix. In our context, Belevitch transformers will not impinge on the
nonsingular character of the Lagrangian, and insofar as they introduce
constraints and therefore redundancies they are taken care of in the
Williamson analysis.

###  III.4 Nonlinear lumped elements

Finally, let us discuss the addition of the nonlinear elements typically used
in superconducting circuits, the Josephson junction (JJ) [5], and its dual
element, the phase-slip junction [66, 67]. The unitary dynamics of the
Josephson junction can be modeled by a linear capacitor in parallel with an
element whose current-flux relation is nonlinear,
$I_{J}=I_{c}\sin(2\pi\Phi_{J}/\Phi_{0})=\partial_{\Phi_{J}}U(\Phi_{J})$. The
complementary phase-slip element also presents a sinusoidal nonlinear
relation, in this case voltage-charge, $V_{PS}=V_{c}\sin(\pi Q_{PS}/2e)$, see
Fig. 5.

Figure 5: Conventional symbols for (a) a Josephson junction lumped element,
with its intrinsic parallel capacitance, (b) and its dual circuit, a phase-
slip junction with a series inductor.

The nonlinear nature of these elements suggests their inclusion in a
Lagrangian description in the form of potential terms

$\displaystyle L$ $\displaystyle=$ $\displaystyle
L_{\mathrm{linear}}-U_{J}(\operatorname{\boldsymbol{\Phi}}_{J})-U_{PS}(\operatorname{\boldsymbol{Q}}_{PS}).$
(19)

No issue arises in moving to a Hamiltonian description, as both flux and
charge variables belong to the configuration space _and_ because a Josephson
(phase-slip) junction will always be in parallel (series) to a capacitor
(inductor), thus ensuring that a kinetic term always exists for the
corresponding variable and mantaining the nonsingular character of the
Lagrangian,

$H=H_{\mathrm{linear}}+U_{J}(\operatorname{\boldsymbol{\Phi}}_{J})+U_{PS}(\operatorname{\boldsymbol{Q}}_{PS}).$
(20)

Once the complete classical and nonlinear Hamiltonian has been obtained,
several different strategies could be applied for its quantum analysis, and
additional considerations come into play. For definiteness, let us focus on
circuits with Josephson junctions. It is now well established that the
topology of the connection of the superconducting islands in a chip can
dramatically change in which way we are to model these nonlinear terms [71].
In the analysis of the Josephson effect one sees that a real junction island
must be modeled by a compact phase variable
$\varphi_{J}=\Phi_{J}\left(\frac{2\pi}{\Phi_{0}}\right)\in[0,2\pi]$, where
$\Phi_{J}$ is the flux variable of the junction and $\Phi_{0}$ the flux
quantum constant. This has an impact on which canonical pair to consider in
quantization, since the boundedness of the phase variable forces its conjugate
operator to have a discrete spectrum, if selfadjointness is to be mantained.
This is conveniently expressed with commutation relations of the form
$[\hat{n}_{J},e^{\pm i\varphi_{J}}]=\pm e^{\pm i\varphi_{J}}$ and the demand
that $\hat{n}_{J}$ have a discrete (integer) spectrum. The experimental data
confirming this idealization is in fact the discreteness of the spectrum of
energies.

However, if the Josephson junction is biased by an inductance, necessarily the
phase variable decompactifies and the spectrum of its quantization is the
whole real line $\varphi\in\mathbb{R}$. The conjugate operator is also a
generator of translations in this case, but its spectrum is not discrete in
this situation. The potential energy term is now a periodic function on the
real line, instead of being defined on the circle. To see the difference in
character, consider two free particle Hamiltonians of identical form, namely
the usual kinetic term, in the two situations of compact and full real line
variable. In the first case there is a discrete spectrum of eigenenergies,
while for the real line all the eigenstates of energy are scattering states,
and the spectrum is continuous.

The first question before us in this regard is whether our proposed analysis
is applicable in either case, and the second one, if indeed applicable, is
what the best strategy can be. The answer to the first question is that
_conceptually_ it is only applicable to _noncompact_ variables. In the form we
have presented the argument to this point we have assumed systematically that
the configuration space is indeed a linear space, in such a way that phase
space is also a linear space and linear symplectic transformations are global
and canonical. If we are to consider more general manifolds it is a fact that
locally we can use Williamson’s theorem to diagonalize symplectically the
quadratic part of the Hamiltonian, but the global validity and usefulness of
the construction remain to be explored, and we set that exploration aside for
future work.

Even within this limited purview, we shall apply the method in situations in
which there are compact variables present. Observe again that the presence of
these compact variables in superconducting circuits is highly dependent on the
topology of the circuit, see Fig. 6. The central physical criterion,
nonetheless, is that a proper superconducting island will be associated with a
compact variable. In the class of circuits we consider this will be ensured if
a Josephson junction is in series with a capacitor as part of a one port
element.

Figure 6: The Josephson junctions can always be modeled with compact variables
in a general network if they appear in pure islands (left). In more general
topological connections (right), it cannot be known a priori.

Now assuming that there are indeed compact variables, we examine two
situations in which we can apply the Williamson analysis. A first case takes
place when the harmonic oscillator approximation is adequate for compact
variables and, as long as the small oscillations condition is maintained, one
can effectively consider the Josephson flux variable as decompactified. One
such situation is that of the transmon qubit, with compact flux variable [72].
The kinetic energy term is due to the Josephson capacity $C_{J}$, and its
characteristic kinetic energy scale, $E_{C}$, behaves as
$E_{C}\propto\frac{1}{\tilde{C}_{J}}$. The potential energy due to the
Josephson inductance presents a characteristic potential energy scale
$E_{J}=I_{c}\left(\frac{\Phi_{0}}{2\pi}\right)$. Under the condition that
$E_{C}\ll E_{J}$ the low energy sector can be studied in the harmonic
approximation, with additional small quartic terms (Duffing oscillator)
$U(\Phi_{J})\approx\frac{\Phi_{J}^{2}}{2L_{J}}-U\frac{\Phi_{J}^{2}}{4!}$, with
$L_{J}^{-1}=E_{J}\left(\frac{2\pi}{\Phi_{0}}\right)^{2}$ and
$U=E_{J}\left(\frac{2\pi}{\Phi_{0}}\right)^{4}$. On coupling to a linear
system, therefore, one can study the quadratic Hamiltonian, with inclusion of
the quadratic part of the Josephson inductance, with the Williamson
diagonalizing procedure. Once that has been carried out one should check that
the weak nonharmonicity is mantained in the new diagonalizing phase space
variables.

An alternative way of stating that this possibility, namely of being able to
apply the symplectic diagonalization procedure to the compact case, is that if
the minima of the compact sector are very deep and essentially remain deep
after the procedure then it can be relied upon. Other than the issue of
compact/noncompact variables, the question remains as to whether the nonlinear
variables are to be kept separate or, rather, taken into account in an
application of the Williamson procedure. The criterion is the usefulness of
either choice in providing us with tools for the quantum treatment of the
system.

In the transmon regime mentioned above, for which the nonharmonicity is weak
in the fashion presented above, it would be sensible to perform Williamson’s
symplectic transformation on the whole set of coordinates in the network,
i.e., including the flux variable of the Josephson junction in the analysis as
a position variable and carrying out the harmonic approximation for its
potential. Thus, in this regime, if we have a network with capacitors,
inductors, NR elements, and JJs, one could start with a maximal set of
coordinates
$\operatorname{\boldsymbol{X}}=(\operatorname{\boldsymbol{\Phi}}_{J},\operatorname{\boldsymbol{\Phi}}_{C},\operatorname{\boldsymbol{Q}}_{L})$,
where $\operatorname{\boldsymbol{\Phi}}_{C}$ are fluxes other than those of
the junctions, perform the harmonic approximation and symplectic
diagonalization, and later check that the excursions of the coordinates inside
the nonlinear terms are indeed small.

A second interesting case in which we can use the quantization procedure
associated with Williamson diagonalization in the presence of compact
variables is when the Lagrangian coupling of those compact fluxes to the rest
of the circuit can be written as involving only their time derivatives.
Observe however that in so doing one must ensure that the compact character is
maintained properly. We will see a complete example of this idea in the next
section. Then we can carry out a Williamson analysis of the Hamiltonian, in
which the compact variables would appear as ignorable, i.e. as free particles,
and after diagonalization add the periodic potentials.

More generally, one could envisage carrying out the linear Williamson analysis
just on the non JJ subsector of the phase-space, and only then consider the
coupling to the Josephson junctions, in all their nonlinearity, to the rest of
the system. In the next section we delve more deeply on the two possibilities,
with the objective of expanding the methodology of blackbox quantization to
generic immitance descriptions.

##  IV Nonreciprocal black-box admittance quantization

As of late there has been an increasing interest in the use of black-box
linear models to describe complex 3D volumes in superconducting circuitry [13,
14, 15, 73, 42, 74]. Most of them rely on having an immitance description of
the linear system, and typically an impedance decomposition is used. This bias
towards one of the descriptions is linked to (i) the generally prevalent
preference of flux coordinates as configuration-space independent variables
for the whole network, which is in turn due to the interest in Josephson
junctions, for which it is most convenient as the nonlinearity is expressed as
a potential energy, and (ii) the fact that the derivation of a Hamiltonian for
the admittance expansion is more easily performed in terms of loop-charge
variables.

In this section we use the systematic Williamson diagonalization, as applied
to an enlarged configuration space, to quantize in a consistent and systematic
way a two-port nonreciprocal black-box generic admittance
$\operatorname{\mathsf{Y}}(\omega)$ capacitively coupled at its ports to
Josephson junctions. We show that it does not present any further obstacle
with respect to its impedance courterpart. In particular, if one were to use a
purely flux-variable based quantization, the Lagrangian of this circuit would
be necessarily written in terms of a singular kinetic (capacitance) matrix. We
sidestep this singularity by working in an enlarged configuration space [35].

The general Foster decomposition of an admittance response is composed of
inductor-gyrator series stages connected to the output ports through a
Belevitch transformer, see Fig. 7. For the sake of simplicity, and without
loss of generality, we restrict the analysis to a 1-stage resonance circuit.
From a physical perspective, the Belevitch transformer finally yields oriented
energy participation ratios between the internal (non)reciprocal modes in the
black-box and the external ports [42].

In line with the discussion above, we quantize along two different routes: (i)
by expanding the cosine potentials of the JJs and diagonalizing the quandratic
part of the Hamiltonian of the whole system, and (ii) by diagonalizing a
harmonic subsector having set aside the JJs coordinates for later (which would
permit the inclusion of the whole cosine in the posterior analysis). We show
the match of both paradigms in the full linear approximation, i.e., two-step
versus a one-step symplectic diagonalizations.

Figure 7: Josephson junctions capacitively coupled to a two-port nonreciprocal
device $\operatorname{\mathsf{Y}}(\omega)$ explicitly showing the (a)
multiport Foster representation of the admittance. (b) The Josephson junctions
are broken into a quadratic ($L_{J}$) and quartic and higher potential terms.

The Lagrangian for Josephson junctions coupled to a generic admittance can be
written by adding up energetic and constraint terms for all the different
lumped elements, i.e., $L=L_{LC}+L_{J}+L_{LT}+L_{\operatorname{\mathsf{Y}}}$.
Here $L_{LC}$, $L_{\operatorname{\mathsf{Y}}}$ and $L_{LT}$ are defined in
Eqs. (9), (13) and (16), respectively, and we collect the JJ terms in $L_{J}$.
By elimination of the trivial constraints, arising from the Belevitch right
fluxes and left charges through the Lagrange multipliers (16), it can be
rewritten as

$\displaystyle L=\,$
$\displaystyle\frac{1}{2}\left(\dot{\operatorname{\boldsymbol{\Phi}}}^{T}\operatorname{\mathsf{C}}\dot{\operatorname{\boldsymbol{\Phi}}}^{T}+\dot{\operatorname{\boldsymbol{Q}}}^{T}\operatorname{\mathsf{L}}\dot{\operatorname{\boldsymbol{Q}}}-2\dot{\operatorname{\boldsymbol{Q}}}^{T}\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}+\dot{\operatorname{\boldsymbol{Q}}}^{T}\operatorname{\mathsf{Z}}\operatorname{\boldsymbol{Q}}\right)$
$\displaystyle-U(\operatorname{\boldsymbol{\Phi}}_{J}),$ (21)

with flux coordinates describing the external variables
$\operatorname{\boldsymbol{\Phi}}^{T}=(\operatorname{\boldsymbol{\Phi}}_{J}^{T},\operatorname{\boldsymbol{\Phi}}_{c}^{T})=(\Phi_{Ja},\Phi_{Jb},\Phi_{ca},\Phi_{cb})$,
and charge coordinates for the internal variables in the admittance
representation $\operatorname{\boldsymbol{Q}}^{T}=(Q_{L1},Q_{L2})$. The
matrices read for the concrete example of Fig. 7 (a)

$\displaystyle\operatorname{\mathsf{C}}$
$\displaystyle=\begin{pmatrix}\operatorname{\mathsf{C}}_{J}&0\\\
0&\operatorname{\mathsf{C}}_{c}\end{pmatrix},\quad\operatorname{\mathsf{L}}=L_{1}\text{1}_{2},\quad\operatorname{\mathsf{D}}=\operatorname{\mathsf{D}}_{R}\operatorname{\mathsf{N}}^{T}\operatorname{\mathsf{D}}_{L}$
(22) $\displaystyle\operatorname{\mathsf{D}}_{L}$
$\displaystyle=\begin{pmatrix}1&0&1&0\\\
0&1&0&1\end{pmatrix}=-\operatorname{\mathsf{D}}_{R},\quad\operatorname{\mathsf{N}}=\begin{pmatrix}n_{11}&n_{12}&0&0\\\
0&0&n_{21}&n_{22}\end{pmatrix},$

where $\operatorname{\mathsf{C}}_{J}=C_{J}\text{1}_{2}$,
$\operatorname{\mathsf{C}}_{c}=C_{c}\text{1}_{2}$, and
$\operatorname{\mathsf{Z}}=-R\operatorname{\mathsf{J}}_{2}$. On the other
hand, the nonlinear potential of the modeling nonlinear Josephson element is
$U(\operatorname{\boldsymbol{\Phi}}_{J})=-\sum_{i}E_{Ji}\cos(\varphi_{Ji})$. A
Legendre transformation is readily available and results in the Hamiltonian

$\displaystyle H=$
$\displaystyle\,\frac{1}{2}(\operatorname{\boldsymbol{P}}-\tfrac{1}{2}\operatorname{\mathsf{Z}}\operatorname{\boldsymbol{Q}}+\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}})^{T}\operatorname{\mathsf{L}}^{-1}(\operatorname{\boldsymbol{P}}-\tfrac{1}{2}\operatorname{\mathsf{Z}}\operatorname{\boldsymbol{Q}}+\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}})$
$\displaystyle+\frac{1}{2}{\operatorname{\boldsymbol{\Pi}}}^{T}\operatorname{\mathsf{C}}^{-1}{\operatorname{\boldsymbol{\Pi}}}^{T}+U(\operatorname{\boldsymbol{\Phi}}_{J})\,,$
(23)

where the conjugate canonical pairs are defined by the Poisson brackets
$\\{\operatorname{\boldsymbol{\Phi}},\operatorname{\boldsymbol{\Pi}}^{T}\\}=\\{\operatorname{\boldsymbol{Q}},\operatorname{\boldsymbol{P}}^{T}\\}=\text{1}$.
Please, notice the similarity between the coupling of branch-fluxes
$\operatorname{\boldsymbol{\Phi}}$ and conjugated flux momenta
$\operatorname{\boldsymbol{P}}$ through the skew-symmetric matrix
$\operatorname{\mathsf{Z}}$ and the minimal-coupling Hamiltonian in the Landau
problem.

Having obtained this Hamiltonian, it is a pedagogical exercise to perform an
analytical study of e Hamiltonian (IV) in the absence of the Josephson
junctions, i.e. setting $U$ to be zero. Given that
$\mathrm{dim}\left(\mathrm{ker}\left[\operatorname{\mathsf{D}}^{T}\right]\right)=0$,
which happens whenever the transformer matrix is nontrivial, the kernels of
$K_{1}$ and $K_{2}$ can be parametrized as

$\displaystyle\operatorname{\boldsymbol{v}}_{1}=\begin{pmatrix}\operatorname{\boldsymbol{q}}\\\
\boldsymbol{\phi}\\\
\frac{\operatorname{\mathsf{Z}}}{2}\operatorname{\boldsymbol{q}}-\operatorname{\mathsf{D}}\boldsymbol{\phi}\\\
0\end{pmatrix}\in
K_{1},\,\,\operatorname{\boldsymbol{v}}_{2}=\begin{pmatrix}\operatorname{\boldsymbol{q}}\\\
\boldsymbol{\phi}\\\
\frac{\operatorname{\mathsf{Z}}}{2}\operatorname{\boldsymbol{q}}-\operatorname{\mathsf{D}}\boldsymbol{\phi}\\\
\operatorname{\mathsf{C}}\operatorname{\boldsymbol{e}}_{\operatorname{\mathsf{D}}}\end{pmatrix}\in
K_{2}\,.$ (24)

Here $\operatorname{\boldsymbol{e}}_{\operatorname{\mathsf{D}}}$ is a generic
element of $\mathrm{ker}\left[\mathsf{D}\right]$, and $\boldsymbol{\phi}$ and
$\operatorname{\boldsymbol{q}}$ span vector spaces of dimension $N_{C}$ and
$M_{L}$, respectively. Note that the coordinates are written in the order
$(\operatorname{\boldsymbol{Q}}^{T},\operatorname{\boldsymbol{\Phi}}^{T},\operatorname{\boldsymbol{P}}^{T},\operatorname{\boldsymbol{\Pi}}^{T})$.
The dimensionality can be directly read as $\mathrm{dim}(K_{1})=N_{C}+M_{L}$,
and
$\mathrm{dim}(K_{2})=N_{C}+M_{L}+\mathrm{dim}\left(\mathrm{ker}\left[\operatorname{\mathsf{D}}\right]\right)$,
such that the number of free particles is
$n_{f}=\mathrm{dim}\left(\operatorname{\boldsymbol{e}}_{\operatorname{\mathsf{D}}}\right)=2$,
and $n_{nd}=\frac{1}{2}(N_{C}+M_{L}-n_{f})=n_{ho}=2$. It can be easily proven
that removing a pair of capacitors, e.g., $C_{ci}$,
$\mathrm{dim}\left(\mathrm{ker}[\operatorname{\mathsf{D}}]\right)=\mathrm{dim}\left(\mathrm{ker}\left[\operatorname{\mathsf{D}}^{T}\right]\right)=0$,
and no free particles would appear in this description with such a choice of
coordinates. In other words, by taking together the limits
$C_{ci}\rightarrow\infty$,
$\operatorname{\boldsymbol{\Phi}}_{c}\rightarrow\operatorname{\boldsymbol{\Phi}}_{J}$
and
$\operatorname{\boldsymbol{\Pi}}_{c}\rightarrow\operatorname{\boldsymbol{\Pi}}_{J}$,
the associated free-particle dynamics vanish, i.e., in algebraic terms
$\mathrm{ker}[\operatorname{\mathsf{JH}}]=\mathrm{ker}[(\operatorname{\mathsf{JH}})^{2}]$.

However, in presence of the Josephson junctions the corresponding fluxes are
compact variables, as signalled above, and one should be more cautious in the
application of the symplectic diagonalization idea. Observe to this point that
the coupling of the Josephson fluxes to the rest of the system is through the
term
$\dot{\operatorname{\boldsymbol{Q}}}^{T}\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\Phi}}$
in Lagrangian (IV), which gives rise to the corresponding couplings in the
Hamiltonian (IV). It is therefore convenient to perform a change of variables
in configuration space, shifting the coupling to the time derivative of the
Josephson fluxes, given by

$\displaystyle\operatorname{\boldsymbol{Q}}=\tilde{\operatorname{\boldsymbol{Q}}}+2\operatorname{\mathsf{Z}}^{-1}\operatorname{\mathsf{D}}_{J}\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}\,.$
(25)

This is accompanied by no change in the flux coordinates,
$\operatorname{\boldsymbol{\Phi}}_{J}=\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}$
and
$\operatorname{\boldsymbol{\Phi}}_{c}=\tilde{\operatorname{\boldsymbol{\Phi}}}_{c}$.
Observe that this point transformation does not change the topology of
configuration space, and it induces a canonical transformation on the full
phase space (whose topology remains unchanged as well), even taking into
account that it is compact along some directions. The change of coordinates in
momentum space reads

$\displaystyle\operatorname{\boldsymbol{P}}$
$\displaystyle=\tilde{\operatorname{\boldsymbol{P}}}\,,$
$\displaystyle\operatorname{\boldsymbol{\Pi}}_{c}$
$\displaystyle=\tilde{\operatorname{\boldsymbol{\Pi}}}_{c}\,,$ (26)
$\displaystyle\operatorname{\boldsymbol{\Pi}}_{J}$
$\displaystyle=\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}+2\operatorname{\mathsf{D}}_{J}^{T}\operatorname{\mathsf{Z}}^{-1}\tilde{\operatorname{\boldsymbol{P}}}\,.$

Under this canonical transformation, the new Josephson flux variables only
appear in the potential energy terms,
$U\left(\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}\right)$, and the
couplings of the junctions to the rest of the system appear only in the
kinetic part of the Hamiltonian.

Let us now apply to this situation the two most common points of view in the
literature regarding quantization of Josephson junctions connected to an
impedance (as opposed to an admittance as we do here). Namely, either (i) a
two tier process, in which we select the part of the Hamiltonian in which
neither the flux variables for the JJs _nor_ their conjugate variables appear,
symplectically diagonalize that part, and then study the full Hamiltonian with
these new variables, or (ii) separating the cosine potentials in a quadratic
and higher order terms, setting aside the quartic and higher terms, computing
the normal modes of the quadratic part of the Hamiltonian (with the quadratic
part of the cosines included), and then expanding the quartic and higher terms
in terms of the new collective normal mode coordinates. This last approach is
known for impedance coupling as the black-box approach [13, 14, 15], see Fig.
7(b). We complete the study by showing the equivalence of both perspectives
under some specific conditions.

To facilitate further analysis of this Hamiltonian, we perform trivial
rescaling (symplectic) transformations in such a way that all the new
coordinates of phase space share the same dimensions. The inductive loop-
charge variables (and their conjugated flux momenta) transform into
$\bar{\operatorname{\boldsymbol{Q}}}=\sqrt{R}\tilde{\operatorname{\boldsymbol{Q}}}$
($\bar{\operatorname{\boldsymbol{p}}}=\tilde{\operatorname{\boldsymbol{p}}}/\sqrt{R}$),
whereas the capacitive branch-fluxes (and conjugated charge momenta) turn into
$\bar{\operatorname{\boldsymbol{\Phi}}}_{c/J}=\operatorname{\mathsf{C}}_{c/J}^{-\frac{1}{4}}\operatorname{\mathsf{L}}^{\frac{1}{4}}\tilde{\operatorname{\boldsymbol{\Phi}}}_{c/J}$
($\bar{\operatorname{\boldsymbol{\Pi}}}_{c/J}=\operatorname{\mathsf{C}}_{c/J}^{\frac{1}{4}}\operatorname{\mathsf{L}}^{-\frac{1}{4}}\tilde{\operatorname{\boldsymbol{\Pi}}}_{c/J}$).
The Hamiltonian with exclusion of the Josephson potentials is then
parameterised in the example of Fig. 7 (a), for these homogeneous dimension
coordinates, in terms of the frequencies $\Omega_{J}=1/\sqrt{C_{J}L}$,
$\Omega_{c}=1/\sqrt{C_{c}L}$, and $\Omega=R/L$. For numerical illustration
(see Appendix E) we set the transformer turn-ratios matrix to
$n_{11}=n_{12}=n_{22}=1$ and $n_{21}=0$, and we work with homogeneous
frequencies $\Omega=\Omega_{c}=\Omega_{J}=1$.

Let us now show explicitly the two tier quantization of this system.
Performing a symplectic transformation that only mixes the subset of
coordinates
$\bar{\operatorname{\boldsymbol{Y}}}=(\bar{\operatorname{\boldsymbol{Q}}},\bar{\operatorname{\boldsymbol{\Phi}}}_{c},\bar{\operatorname{\boldsymbol{P}}},\bar{\operatorname{\boldsymbol{\Pi}}}_{c})$,
and diagonalizes $\operatorname{\mathsf{H}}_{\mathrm{lin}}$, i.e.,
$(\operatorname{\mathsf{S}}_{\operatorname{\boldsymbol{Y}}}^{-1}\operatorname{\boldsymbol{Y}})^{T}=(\bar{\operatorname{\boldsymbol{y}}}_{\mathrm{nd}}^{T},\bar{\operatorname{\boldsymbol{\pi}}}_{\mathrm{nd}}^{T},\bar{\operatorname{\boldsymbol{y}}}_{\mathrm{ho}}^{T},\bar{\operatorname{\boldsymbol{\pi}}}_{\mathrm{ho}}^{T})$,we
obtain (two) harmonic oscillators linearly coupled (through momenta) to the
(dressed) nonlinear system

$\displaystyle H=$
$\displaystyle\sum_{n\in{ho}}^{2}\frac{\Omega_{n}}{2}(\bar{y}_{n}^{2}+\bar{\pi}_{n}^{2})+\left(\bar{\operatorname{\boldsymbol{y}}}^{T}\mathsf{G}+\bar{\operatorname{\boldsymbol{\pi}}}^{T}\mathsf{M}\right)\bar{\operatorname{\boldsymbol{\Pi}}}_{J}$
$\displaystyle+\sum_{n=1}^{2}\left[\Omega_{J}\bar{\Pi}_{Jn}^{2}-E_{Jn}\cos(\bar{\varphi}_{Jn})\right],$
(27)

where $\\{\bar{y}_{i},\bar{\pi}_{j}\\}=\delta_{ij}$, and $\mathsf{G}$ and
$\mathsf{M}$ are coupling matrices between the Josephson momenta and the
partially-diagonalized coordinates. A canonical quantization of this system
can now be performed given that the nondynamical sector has been singled out
and set to zero [70], promoting to quantum operators only the dynamical ones.
After reexpressing the harmonic conjugate variables in terms of annihilation
and creation operators, we obtain the Hamiltonian

$\displaystyle\hat{H}/\hbar=$
$\displaystyle\sum_{n,m}^{2}\left[\Omega_{n}a_{n}^{\dagger}a_{n}+(g_{nm}a_{n}+g_{nm}^{*}a_{n}^{\dagger})\bar{\Pi}_{Jm}\right]+H_{J},$

where $H_{J}$ stands for all the terms involving only the Josephson sector in
(27), namely, the second line. Specific values for eigenfrequencies and
coupling vectors for this set of parameters can be seen in Appendix E.

If on the other hand a black-box approach [13] is preferred, we separate out
the quartic and higher terms of the Taylor expansion of the cosine potentials
of the junctions. The remaining quadratic part (in which a
$\bar{\operatorname{\boldsymbol{\Phi}}}_{J}^{T}\operatorname{\mathsf{L}}_{J}^{-1}\bar{\operatorname{\boldsymbol{\Phi}}}_{J}$
term appears) is susceptible to symplectic diagonalization. Observe however
that the compactness of the Josephson fluxes has been discarded, and only in a
transmonlike regime would this process be valid. After this symplectic
diagonalization, quantize the resulting (four) harmonic oscillators, and
expand the nonlinear potential in terms of normal modes to obtain the
Hamiltonian

$\displaystyle\hat{H}_{bb}=\sum_{n}^{4}\hbar\bar{\Omega}_{n}b_{n}^{\dagger}b_{n}+E_{J}\sum_{l={1,2}}\sum_{m=2}^{\infty}\frac{(-1)^{m+1}\bar{\varphi}_{Jl}^{2m}}{(2m)!}.$
(28)

It is worth remarking that the flux variables of the Josephson junctions
depend in this case only on the dynamical sector after the appropiate
symplectic transformations, i.e.,
$\bar{\varphi}_{Ji}\equiv\bar{\varphi}_{Ji}(\operatorname{\boldsymbol{z}}_{ho},\operatorname{\boldsymbol{\pi}}_{ho})$,
with
$\operatorname{\boldsymbol{Z}}=\operatorname{\mathsf{S}}^{-1}\bar{\operatorname{\boldsymbol{X}}}=(\operatorname{\boldsymbol{z}}_{nd},\operatorname{\boldsymbol{z}}_{ho},\operatorname{\boldsymbol{\pi}}_{nd},\operatorname{\boldsymbol{\pi}}_{ho})$
thus setting the correct grounds for a low-lying energy investigation.
Naturally, we remind the reader that the previously described routines match,
as it can be easily proven by expanding to second order the cosine in Eq. (27)
and performing a secondary symplectic transformation, see Appendix E for
details.

As a final note, we observe that in more general scenarios the Josephson
variables might indeed depend on the nondynamical coordinates of the linear
problem, i.e.,
$\bar{\varphi}_{Ji}(\operatorname{\boldsymbol{z}}_{nd}^{T},\operatorname{\boldsymbol{\pi}}_{nd}^{T},\operatorname{\boldsymbol{z}}_{ho}^{T},\operatorname{\boldsymbol{\pi}}_{ho}^{T})$.
Thus the black-box approach as presented here would not set out an immediately
obvious path for the investigation of the low-lying energy sector, and the two
tier process, e.g., Eq. (27), should always be preferred. Furthermore, the
problem here stated would indeed increase in complexity if the Hamiltonian had
extra time-dependent classical forces, as it is the case when dealing with
pulses for quantum computation or readout [75, 76, 77], and the two tier
picture is recommended as the basis for further analysis.

##  V Conclusions and outlook

In this article we have put forward a consistent quantization procedure for
superconducting lumped-element networks modeled in an enlarged configuration
space description, with both flux and charge variables. The manifest
redundancy is eliminated in the Hamiltonian through a symplectic
transformation provided systematically by Williamson’s theorem, for which we
have given a constructive algorithm for real positive-semidefinite Hamiltonian
matrices. To the best of our knowledge, this is its first explicit application
to a simple physical scenario. The full-fledged Williamson’s theorem has been
recurrently referenced in the quantum physics community, but only positive-
definite matrices have been explicitly provided with either proof and explicit
construction [78, 79, 80]. Given its symplicity, this algorithm will find use
in a plethora of gaussian systems beyond electrical circuits.

In the context of the specific problem here tackled, we emphasize that the
enlarged configuration space idea becomes a useful starting point for the
identification and quantization of the dynamical degrees of freedom in
circuits with frequency-dependent nonreciprocal devices where a mixed charge-
flux configuration space variables may result inevitable. This reinforces the
idea we put forward for circuits with transmission lines coupled through ideal
nonreciprocal devices [70]. Naturally, our algorithm will also prove useful to
single out free particle dynamics in reciprocal circuits where a redundant
choice of either flux or charge variables is chosen as the configuration
space.

As an example of the power of our technique, we have reduced the problem of
quantizing a circuit composed of Josephson junctions capacitively coupled to a
general linear system described in terms of an admittance matrix to that of
correctly identifying the dynamical and nondynamical sectors of a redundant
Hamiltonian. Such analysis would have been much more cumbersome had we started
the analysis with a Lagrangian written in terms of just flux variables,
because the kinetic matrix would be singular. In fact, we have presented and
compared the two most standard methods to quantize canonically superconducting
circuits, i.e., performing: (i) a partial diagonalization of the purely linear
sector, setting aside at that point of the analysis the full nonlinear cosine
potentials of the Josephson junctions, and (ii) a full diagonalization of the
linear sector within a generalized black-box quantization paradigm. We point
out that in some cases the black-box approach may require further refinement
when the nonlinear potentials are written in the transformed basis. This would
arise from a (possibly) nonlinear coupling between those sectors that in the
linear approximation had been identified as dynamical and nondynamical,
respectively.

Looking to the future, further work will be required to extend these ideas for
an exact quantization prescription in the case of quasi-lumped circuits
containing frequency-dependent nonreciprocal black-boxes, as well as
transmission lines and nonlinear elements, along the lines first suggested in
[52].

###### Acknowledgements.

I. L. E. acknowledges support of the Basque Government grant IT986-16. A.
P.-R. thanks Basque Government Ph.D. grant PRE-2016-1-0284, and the Canada
First Research Excellence Fund. We thank C. Lledó-Veloso for a first
proofreading the article.

##  Appendix A Symplectic diagonalization of positive semidefinite
Hamiltonians

As introduced in the text, we investigate the symplectic diagonalization of
(even dimensional) positive semidefinite matrices. We shall denote them as
$\mathsf{H}$, and the canonical symplectic matrix $\mathsf{J}$ is

$\mathsf{J}=\begin{pmatrix}0&\text{1}\\\ -\text{1}&0\end{pmatrix}\,.$ (29)

If necessary we will distinguish canonical symplectic matrices for different
dimensions with explicit mention of the dimension.

Williamson’s theorem concerns the canonical forms under symplectic
transformations of even dimensional matrices over fields of non zero
characteristic. The crucial fact underlying the full Williamson’s theorem and
its restrictions is that the conjugacy normal forms of $\mathsf{JH}$ determine
completely the symplectic normal forms of $\mathsf{H}$. This comes about
because the space of states is a direct sum of symplectically orthogonal
subspaces, and those are associated with generalised eigenspaces of
$\mathsf{JH}$.

The detailed proof of Williamson’s theorem for the real field consists in the
analysis of each possible structure of generalised eigenspaces and the
determination of the corresponding canonical form for that sector. We refer to
[50] for the complete proof. Here we simply give a general guide to the
concepts involved in the next subsection, and later on we will concentrate on
the case of interest, that of positive semidefinite Hamiltonians, taking the
theorem as a starting point.

###  A.1 Generalities of Williamson’s theorem for real symmetric matrices

For real symmetric matrices, as is the case of interest to us, the
characteristic polynomial of $\mathsf{JH}$, $p(\lambda)$, has real
coefficients and is even, $p(-\lambda)=p(\lambda)$. This simple fact comes
about because $\mathsf{J}^{2}=-\text{1}$ and $\mathrm{det}(\mathsf{J})=1$.
Thus,

$\displaystyle p(\lambda)$
$\displaystyle=\mathrm{det}\left[\mathsf{JH}-\lambda\text{1}\right]=\mathrm{det}\left[\mathsf{H}+\lambda\mathsf{J}\right]$
$\displaystyle=\mathrm{det}\left[\left(\mathsf{H}+\lambda\mathsf{J}\right)^{T}\right]=\mathrm{det}\left[\mathsf{H}-\lambda\mathsf{J}\right]$
$\displaystyle=p(-\lambda)\,.$ (30)

Therefore, if $\lambda$ is an eigenvalue (possibly generalised), then so are
$-\lambda$, the complex conjugate $\lambda^{*}$ and $-\lambda^{*}$. Thus there
are in principle four types of eigenvalue: i) the generic $\mu+i\nu$ with real
$\mu$ and $\nu$, in a quartet
$\left\\{\mu+i\nu,\mu-i\nu,-\mu-i\nu,-\mu+i\nu\right\\}$; ii) purely imaginary
$\pm i\mu$; iii) purely real $\pm\mu$; and iv) 0.

We shall use the notation $\langle\bullet,\bullet\rangle$ for the canonical
inner product both for the real space $\mathbb{R}^{2n}$ on which $\mathsf{H}$
and $\mathsf{J}$ act and for the sesquilinear inner product of its
complexification and of $\mathbb{C}^{2n}$. In matrix terms,
$\langle\mathbf{u},\mathbf{v}\rangle=\mathbf{u}^{\dagger}\mathbf{v}$, with
${\dagger}$ denoting hermitian conjugate.

Let us now consider $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ proper eigenvectors
of $\mathsf{JH}$, with eigenvalues $\lambda_{1}$ and $\lambda_{2}$
respectively. Then

$\displaystyle\left(\lambda_{2}+\lambda_{1}^{*}\right)\langle\mathbf{v}_{1},\mathsf{J}\mathbf{v}_{2}\rangle$
$\displaystyle=\left\langle\mathbf{v}_{1},\mathsf{J}\left(\mathsf{JH}\mathbf{v}_{2}\right)\right\rangle+\left\langle\mathsf{JH}\mathbf{v}_{1},\mathsf{J}\mathbf{v}_{2},\right\rangle$
$\displaystyle=-\left\langle\mathbf{v}_{1},\mathsf{H}\mathbf{v}_{2}\right\rangle+\left\langle\mathsf{H}\mathbf{v}_{1},\mathbf{v}_{2}\right\rangle=0\,.$
(31)

Therefore, if $\lambda_{1}^{*}+\lambda_{2}\neq 0$ the corresponding two
eigenvectors are symplectic orthogonal. In fact, we can extend this property
to the generalised eigenspaces if there are nontrivial Jordan chains, by
induction. To this point, let us consider that the respective Jordan chains
are of the form

$\mathsf{JH}\mathbf{v}_{i}^{(k)}=\lambda_{i}\mathbf{v}_{i}^{(k)}+\mathbf{v}_{i}^{(k-1)}\,.$
(32)

Then, along the lines of Eq. (A.1), we have

$\displaystyle\left(\lambda_{2}+\lambda_{1}^{*}\right)\langle\mathbf{v}_{1}^{(k)},\mathsf{J}\mathbf{v}^{(l)}_{2}\rangle=\left\langle\mathbf{v}^{(k)}_{1},\mathsf{J}\left(\mathsf{JH}\mathbf{v}^{(l)}_{2}-\mathbf{v}_{2}^{(l-1)}\right)\right\rangle$
$\displaystyle+\left\langle\left(\mathsf{JH}\mathbf{v}^{(k)}_{1}-\mathbf{v}^{(k-1)}_{1}\right),\mathsf{J}\mathbf{v}^{(l)}_{2},\right\rangle$
$\displaystyle=-\left\langle\mathbf{v}^{(k)}_{1},\mathsf{H}\mathbf{v}^{(l)}_{2}\right\rangle+\left\langle\mathsf{H}\mathbf{v}^{(k)}_{1},\mathbf{v}^{(l)}_{2}\right\rangle$
$\displaystyle-\left\langle\mathbf{v}_{1}^{(k)},\mathsf{J}\mathbf{v}_{2}^{(l-1)}\right\rangle-\left\langle\mathbf{v}_{1}^{(k-1)},\mathsf{J}\mathbf{v}_{2}^{(l)}\right\rangle$
$\displaystyle=-\left\langle\mathbf{v}_{1}^{(k)},\mathsf{J}\mathbf{v}_{2}^{(l-1)}\right\rangle-\left\langle\mathbf{v}_{1}^{(k-1)},\mathsf{J}\mathbf{v}_{2}^{(l)}\right\rangle\,,$
(33)

and, by induction, the respective chains are symplectic orthogonal if
$\lambda_{2}+\lambda_{1}^{*}\neq 0$.

Please observe that $\langle\bullet,\mathsf{J}\bullet\rangle$ is _not_
antisymmetric for complex vectors, namely, given two real vectors $\mathbf{u}$
and $\mathbf{v}$ one has

$\left\langle\mathbf{u}+i\mathbf{v},\mathsf{J}\left(\mathbf{u}+i\mathbf{v}\right)\right\rangle=2i\left\langle\mathbf{u},\mathsf{J}\mathbf{v}\right\rangle\,.$
(34)

Let $V_{\lambda}$ be the generalised eigenspace of $\mathsf{JH}$ for the
eigenvalue $\lambda$. Then the result above is that the full complexified
vector space $V=\mathbb{R}^{2n}_{\mathbb{C}}$ can be written as a direct sum
of symplectic orthogonal subspaces,

$V=V_{0}\oplus\bigoplus_{\lambda\neq 0}\left(V_{\lambda}\oplus
V_{-\lambda^{*}}\right)\,,$ (35)

where the sum avoids double counting, i.e. if $\lambda$ is in the index set,
then $-\lambda^{*}$ is not.

The relevant issue at this point is the construction of a canonical symplectic
basis using these generalised eigenspaces. We remind the reader that a
symplectic basis of $\left\\{\mathbb{R}^{2n},\mathsf{J}\right\\}$ is a basis
with elements
$\left\\{\mathbf{e}_{1},\mathbf{e}_{2},\ldots,\mathbf{e}_{n},\mathbf{f}_{1},\ldots,\mathbf{f}_{n}\right\\}$
such that

$\displaystyle\left\langle\mathbf{e}_{i},\mathsf{J}\mathbf{e}_{j}\right\rangle$
$\displaystyle=0\,,$
$\displaystyle\left\langle\mathbf{f}_{i},\mathsf{J}\mathbf{f}_{j}\right\rangle$
$\displaystyle=0\,,$ (36)
$\displaystyle\left\langle\mathbf{e}_{i},\mathsf{J}\mathbf{f}_{j}\right\rangle$
$\displaystyle=\delta_{ij}\,.$

As $\mathsf{JH}$ is real, the generalised eigenvector $\mathbf{v}$
corresponding to a complex eigenvalue $\lambda$ has an associated generalised
eigenvector corresponding to the complex conjugate eigenvalue $\lambda^{*}$,
namely $\mathbf{v}^{*}$. Thus, in the construction of the _real_ symplectic
basis we desire we shall have to use combinations of $V_{\lambda}$ and
$V_{\lambda^{*}}$. Summarizing, there are four types of symplectic subspaces
to consider. First, the generic $\lambda=\mu+i\nu$ with $\mu,\nu>0$, for which
the symplectic space to be considered is the direct sum

$V_{\lambda}\oplus V_{-\lambda^{*}}\oplus V_{\lambda^{*}}\oplus
V_{-\lambda}\,.$ (37)

Secondly, the purely imaginary case for which we have the symplectic space
$V_{i\mu}\oplus V_{-i\mu}$. Thirdly, the purely real case $V_{\mu}\oplus
V_{-\mu}$, and finally $V_{0}$. As stated above, the theorem relies on the
detailed analysis of these four cases separately.

In the construction of the relevant symplectic basis it is important to notice
that the generalised eigenspaces of the different elements of the quartet (or
doublet) will have the same structure. That is to say, if $\lambda$ is a
generalised eigenvalue with a Jordan chain of length $n_{\lambda}$, then
$-\lambda$, $\lambda^{*}$ and $-\lambda^{*}$ will all have Jordan chains of
the same length, and similarly for the doublets. Similarly, if $\lambda$ is
degenerate then the degeneracy structure is shared by the other elements of
the multiplet.

For definiteness, consider now a nondegenerate quartet, with shared Jordan
chain length $n_{\lambda}$. We have a set of generalised eigenvectors of
$\mathsf{JH}$, obeying the Jordan chain structure of Eq. (32), that we will
denote as
$\left\\{\mathbf{v}^{(k)},\mathbf{u}^{(k)},\bar{\mathbf{v}}^{(k)},\bar{\mathbf{u}}^{(k)}\right\\}$,
corresponding to $\lambda,-\lambda^{*},\lambda^{*},-\lambda$ respectively. We
use here $\bar{\mathbf{u}}$ to denote the complex conjugate in order to avoid
an overburden of superscripts. According to the general result above,
$\left\langle\mathbf{v}^{(k)},\mathsf{J}\mathbf{u}^{(k)}\right\rangle$ is not
trivial, but
$\left\langle\mathbf{v}^{(k)},\mathsf{J}\bar{\mathbf{u}}^{(k)}\right\rangle=\left\langle\mathbf{v}^{(k)},\mathsf{J}\bar{\mathbf{v}}^{(k)}\right\rangle=\left\langle\mathbf{u}^{(k)},\mathsf{J}\bar{\mathbf{u}}^{(k)}\right\rangle=0$.
Let us choose to normalize $\mathbf{v}^{(k)}$ and $\mathbf{u}^{(k)}$,
including phases, such that

$\left\langle\mathbf{v}^{(k)},\mathsf{J}\mathbf{u}^{(k)}\right\rangle=2\,.$
(38)

We now define the real vectors

$\displaystyle\mathbf{e}_{2k-1}$
$\displaystyle=\frac{1}{2}\left[\mathbf{v}^{(k)}+\bar{\mathbf{v}}^{(k)}\right]\,,$
$\displaystyle\mathbf{e}_{2k}$
$\displaystyle=\frac{-i}{2}\left[\mathbf{v}^{(k)}-\bar{\mathbf{v}}^{(k)}\right]\,,$
$\displaystyle\mathbf{f}_{2k-1}$
$\displaystyle=\frac{1}{2}\left[\mathbf{u}^{(k)}+\bar{\mathbf{u}}^{(k)}\right]\,,$
(39) $\displaystyle\mathbf{f}_{2k}$
$\displaystyle=\frac{-i}{2}\left[\mathbf{u}^{(k)}-\bar{\mathbf{u}}^{(k)}\right],.$

It is easy to check that they form a symplectic basis for this level. The only
remaining task to handle is to complete Gram–Schmidt symplectic
orthogonalization (sGS, in what follows) for the whole Jordan chain. The
change of basis matrix provided by the components of the symplectic basis is a
symplectic matrix, by construction, and provides us with the canonical form of
the corresponding block.

If the eigenvalue is real, the generalised eigenvectors are real, and one
carries out sGS similarly for each sector. We set aside the case of imaginary
eigenvalues and that of eigenvalue zero, since they are relevant to the
positive semidefinite case that we investigate in the next subsection.

###  A.2 Positive semidefinite Hamiltonians

Let us now address the case of interest to us here, that of positive
semidefinite Hamiltonians. In this situation there are only three possible
types of canonical symplectic blocks, namely i) simple harmonic oscillators,
with canonical form

$\mathsf{H}_{ho}=\begin{pmatrix}\omega&0\\\ 0&\omega\end{pmatrix}\,,$ (40)

ii) free particles, with canonical form

$\mathsf{H}_{f}=\begin{pmatrix}0&0\\\ 0&1\end{pmatrix}\,,$ (41)

and iii) nondynamical variables with identically zero blocks. The Hamiltonian
will thus be symplectic diagonalizable.

The issue for an effective algorithm applicable to the cases of interest is
the disentangling of the free particle sector. It arises in the difference
between the kernel of $\mathsf{JH}$,
$K_{1}=\mathrm{ker}\left[\mathsf{JH}\right]$ and that of
$K_{2}=\mathrm{ker}\left[\left(\mathsf{JH}\right)^{2}\right]$. For this case
of positive semidefinite $\mathsf{H}$ the generalized eigenspace of eigenvalue
0, $V_{0}$, is precisely $K_{2}$. In other words, the maximum length of a
Jordan chain for this eigenvalue is two.

If $K_{2}=K_{1}$ there are no free particles, and $\mathsf{JH}$ is
diagonalizable. In this situation the speediest way of achieving our goal of
symplectic diagonalization of $\mathsf{H}$ in the relevant sector is to
compute $\mathrm{ker}\left[\mathsf{H}\right]=K_{1}$, and thus the orthogonal
decomposition $\mathbb{R}^{2n}=K_{1}\oplus K_{1}^{\perp}$. The restriction of
$\mathsf{H}$ to $K_{1}^{\perp}$, denoted here as $\mathsf{H}_{+}$, is positive
definite, and all the existing methods apply (see for instance section 8.3 of
[80]).

Let us examine the general case, which will provide us with one such method
for the simpler situation as well. Define the linear subspace

$E=\mathsf{JH}\left[K_{2}\right]\subseteq K_{1}\,,$ (42)

and identify formally

$\tilde{F}=\left(\mathsf{JH}\right)^{-1}\left[E\right]\,.$ (43)

To be definite, this means the following: obtain a basis of $E$,
$\left\\{\tilde{\mathbf{e}}_{I}\right\\}_{I=1}^{n_{f}}$, and assign to each
element of the basis its Jordan chain, i.e. a vector $\tilde{\mathbf{f}}_{I}$
such that $\mathsf{JH}\tilde{\mathbf{f}}_{I}=\tilde{\mathbf{e}}_{I}$.
$\tilde{F}$ is the linear span of these vectors,

$\tilde{F}=\mathrm{span}\left\\{\tilde{\mathbf{f}}_{I}\right\\}_{I=1}^{n_{f}}\,.$
(44)

Complete $\left\\{\tilde{\mathbf{e}}_{I}\right\\}_{I=1}^{n_{f}}$ to a basis of
$K_{1}$ with a set of vectors
$\left\\{\tilde{\mathbf{w}}_{i}\right\\}_{i=1}^{2n_{nd}}$, and denote the span
of these as $\tilde{W}$

$\tilde{W}=\mathrm{span}\left\\{\tilde{\mathbf{w}}_{i}\right\\}_{i=1}^{2n_{nd}}\,.$
(45)

It is important to notice that
$\left\langle\tilde{\mathbf{w}}_{i},\mathsf{J}\tilde{\mathbf{e}}_{I}\right\rangle=0$
($\tilde{W}$ belongs to the symplectic complement of $E$), and, similarly,
$\left\langle\tilde{\mathbf{e}}_{I},\mathsf{J}\tilde{\mathbf{e}}_{J}\right\rangle=0$
($E$ is isotropic), because

$\displaystyle\left\langle\tilde{\mathbf{w}}_{i},\mathsf{J}\tilde{\mathbf{e}}_{I}\right\rangle$
$\displaystyle=\left\langle\tilde{\mathbf{w}}_{i},\mathsf{J}\left[\mathsf{JH}\tilde{\mathbf{f}}_{I}\right]\right\rangle$
$\displaystyle=-\left\langle\tilde{\mathbf{w}}_{i},\mathsf{H}\tilde{\mathbf{f}}_{I}\right\rangle$
$\displaystyle=-\left\langle\mathsf{H}\tilde{\mathbf{w}}_{i},\tilde{\mathbf{f}}_{I}\right\rangle=0\,.$
(46)

We now carry out sGS for the free particle sector, making use of the fact that
$E$ is isotropic, defining recursively for $I=1$ to $n_{f}$

$\displaystyle\mathbf{e}_{I}$
$\displaystyle=a_{I}\left[\tilde{\mathbf{e}}_{I}-\sum_{J<I}\left\langle\tilde{\mathbf{e}}_{I},\mathsf{J}\mathbf{f}_{J}\right\rangle\mathbf{e}_{J}\right]\,,$
(47) $\displaystyle\mathbf{f}_{I}$
$\displaystyle=a_{I}\left[\tilde{\mathbf{f}}_{I}-\sum_{J<I}\left\langle\tilde{\mathbf{f}}_{I},\mathsf{J}\mathbf{f}_{J}\right\rangle\mathbf{e}_{J}+\sum_{J<I}\left\langle\tilde{\mathbf{f}}_{I},\mathsf{J}\mathbf{e}_{J}\right\rangle\mathbf{f}_{J}\right]\,.$

The normalization factors are determined to be

$a_{I}^{-2}=\left\langle\tilde{\mathbf{f}}_{I},\mathsf{H}\tilde{\mathbf{f}}_{I}\right\rangle-\sum_{J<I}\left\langle\tilde{\mathbf{e}}_{I},\mathsf{J}\mathbf{f}_{J}\right\rangle\left\langle\mathbf{e}_{J},\mathsf{J}\tilde{\mathbf{f}}_{I}\right\rangle\,.$
(48)

By induction one can prove that for all $I$ from 1 to $n_{f}$
$\mathsf{JH}\mathbf{f}_{I}=\mathbf{e}_{I}$, so the canonical Jordan form is
preserved in the new basis for this sector. Observe that $E$ is also the span
of $\left\\{\mathbf{e}_{I}\right\\}_{I=1}^{n_{f}}$, and with the sGS we have
constructed a linear space
$F=\mathrm{span}\left\\{\mathbf{f}_{I}\right\\}_{I=1}^{n_{f}}$ that satisfies
both the condition $\mathsf{JH}F=E$ and that of forming the symplectic
subspace $W_{F}=E\oplus F$. Now we want to obtain the symplectic complement of
$W_{F}$ in $K_{2}$. To that purpose, we have to substract from $\tilde{W}$ the
symplectic projection of $W_{F}$ on it, defining

$\mathbf{w}_{i}=\tilde{\mathbf{w}}_{i}-\sum_{J=1}^{n_{f}}\left\langle\tilde{\mathbf{w}}_{i},\mathsf{J}\mathbf{f}_{J}\right\rangle\mathbf{e}_{J}$
(49)

and their linear span
$W=\mathrm{span}\left\\{\mathbf{w}_{i}\right\\}_{i=1}^{2n_{nd}}$. This is a
symplectic subspace, and it is symplectic orthogonal to $W_{f}$ and to the
harmonic oscillator sector.

The only complication in the harmonic oscillator sector might come from the
presence of degenerate frequencies, i.e. from the presence of degenerate
$i\omega_{a}$ eigenvalues of $\mathsf{JH}$. Let us denote with
$\mathbf{v}^{(\alpha)}_{a}$ an eigenvector of $\mathsf{JH}$ for eigenvalue
$i\omega_{a}$, where the superscript $\alpha$ is a degeneracy index. Then
$\bar{\mathbf{v}}^{(\alpha)}_{a}$ is an eigenvector of $\mathsf{JH}$ for
eigenvalue $-i\omega_{a}$. Furthermore, let $V_{a}$ be the (proper) eigenspace
of $\mathsf{JH}$ for eigenvalue $i\omega_{a}$, and $\bar{V}_{a}$ its complex
conjugate. Then $V_{a}\oplus\bar{V}_{a}$ is a symplectic subspace.
Alternatively, we can write this as
$\mathrm{Re}\left(V_{a}\right)\otimes\mathrm{Im}\left(V_{a}\right)$ as a real
symplectic subspace. Thus the only remaining task is to construct a symplectic
basis of this subspace, again by sGS. Namely, define

$\displaystyle\tilde{\mathbf{e}}_{a}^{(\alpha)}$
$\displaystyle=\mathrm{Re}\left(\mathbf{v}_{a}^{(\alpha)}\right)\,,$
$\displaystyle\tilde{\mathbf{f}}_{a}^{(\alpha)}$
$\displaystyle=\mathrm{Im}\left(\mathbf{v}_{a}^{(\alpha)}\right)\,.$ (50)

We perform this assignment because
$\left\langle\mathrm{Re}\left(\mathbf{v}_{a}^{(\alpha)}\right),\mathsf{J}\,\mathrm{Im}\left(\mathbf{v}_{a}^{(\alpha)}\right)\right\rangle$
is positive. This statement might seem surprising, since one can change the
phase of $\mathbf{v}_{a}^{(\alpha)}$ ad libitum. However, as
$\mathsf{JH}\mathbf{v}_{a}^{(\alpha)}=i\omega_{a}\mathbf{v}_{a}^{(\alpha)}$,
it follows that

$\displaystyle\mathsf{JH}\,\mathrm{Re}\left(\mathbf{v}_{a}^{(\alpha)}\right)$
$\displaystyle=-\omega_{a}\,\mathrm{Im}\left(\mathbf{v}_{a}^{(\alpha)}\right)\,,$
$\displaystyle\mathsf{JH}\,\mathrm{Im}\left(\mathbf{v}_{a}^{(\alpha)}\right)$
$\displaystyle=\omega_{a}\,\mathrm{Re}\left(\mathbf{v}_{a}^{(\alpha)}\right)\,,$
(51)

whence

$\displaystyle\omega_{a}\left\langle\tilde{\mathbf{e}}_{a}^{(\alpha)},\mathsf{J}\,\tilde{\mathbf{f}}_{a}^{(\alpha)}\right\rangle$
$\displaystyle=-\left\langle\mathrm{Re}\left(\mathbf{v}_{a}^{(\alpha)}\right),\tilde{\mathbf{e}}_{a}^{(\alpha)}\mathsf{J}\,\mathsf{JH}\left[\tilde{\mathbf{e}}_{a}^{(\alpha)}\right]\right\rangle$
$\displaystyle=\left\langle\tilde{\mathbf{e}}_{a}^{(\alpha)},\mathsf{H}\,\tilde{\mathbf{e}}_{a}^{(\alpha)}\right\rangle>0\,.$
(52)

Now, from this assignment, and as usual, define recursively

$\displaystyle\mathbf{e}_{a}^{(a)}$
$\displaystyle=A_{a}^{(\alpha)}\left[\tilde{\mathbf{e}}_{a}^{(a)}-\sum_{\beta<\alpha}\left\langle\tilde{\mathbf{e}}_{a}^{(\alpha)},\mathsf{J}\mathbf{f}_{a}^{(\beta)}\right\rangle\mathbf{e}_{a}^{(\beta)}\right.$
$\displaystyle\quad\left.+\sum_{\beta<\alpha}\left\langle\tilde{\mathbf{e}}_{a}^{(\alpha)},\mathsf{J}\mathbf{e}_{a}^{(\beta)}\right\rangle\mathbf{f}_{a}^{(\beta)}\right]\,,$
(53)

and similarly for $\mathbf{f}_{a}^{(\alpha)}$, with the normalization selected
to ensure that we obtain a symplectic basis. We choose the same factor
$A_{a}^{(\alpha)}$ for both $\mathbf{e}_{a}^{(\alpha)}$ and
$\mathbf{f}_{a}^{(\alpha)}$, to ensure that their combination
$\mathbf{e}_{a}^{(\alpha)}+i\mathbf{f}_{a}^{(\alpha)}$ is also an eigenvector
of $\mathsf{JH}$ and therefore that the relations (A.2) hold as

$\displaystyle\mathsf{JH}\,\mathbf{e}_{a}^{(\beta)}$
$\displaystyle=-\omega_{a}\mathbf{f}_{a}^{(\beta)}\,,$
$\displaystyle\mathsf{JH}\,\mathbf{f}_{a}^{(\beta)}$
$\displaystyle=\omega_{a}\mathbf{e}_{a}^{(\beta)}\,.$ (54)

We now prove this statement inductively. To do so, define

$\mathbf{u}_{a}^{(\alpha)}=\mathbf{e}_{a}^{(\alpha)}+i\mathbf{f}_{a}^{(\alpha)}\,.$
(55)

If the inductive hypothesis Eq. (A.2) holds for $\beta<\alpha$, then for those
$\beta<\alpha$

$\mathsf{JH}\,\mathbf{u}_{a}^{(\beta)}=i\omega_{a}\mathbf{u}_{a}^{(\beta)}\,,$
(56)

and therefore

$\mathsf{JH}\,\bar{\mathbf{u}}_{a}^{(\beta)}=-i\omega_{a}\bar{\mathbf{u}}_{a}^{(\beta)}\,.$
(57)

Consequently,

$\left\langle\mathbf{v}_{a}^{(\alpha)},\mathsf{J}\,\bar{\mathbf{u}}_{a}^{(\beta)}\right\rangle=0\,,$
(58)

as follows from Eq. (A.1), from which one can conclude that

$\displaystyle\left\langle\tilde{\mathbf{e}}_{a}^{(\alpha)},\mathsf{J}\,{\mathbf{e}}_{a}^{(\beta)}\right\rangle-\left\langle\tilde{\mathbf{f}}_{a}^{(\alpha)},\mathsf{J}\,{\mathbf{f}}_{a}^{(\beta)}\right\rangle$
$\displaystyle=0\,,$
$\displaystyle\left\langle\tilde{\mathbf{e}}_{a}^{(\alpha)},\mathsf{J}\,{\mathbf{f}}_{a}^{(\beta)}\right\rangle+\left\langle\tilde{\mathbf{f}}_{a}^{(\alpha)},\mathsf{J}\,{\mathbf{e}}_{a}^{(\beta)}\right\rangle$
$\displaystyle=0$ (59)

for $\beta<\alpha$ if the inductive hypothesis holds. Clearly Eq. (A.2) holds
for $\beta=1$, the first inductive step. It is now easy to check that if it
holds for $\beta<\alpha$, then Eq. (A.2) and its $\mathbf{f}$ companion allow
us to conclude that (A.2) also holds for $\beta=\alpha$.

Now we complete our task. The matrix constructed as

$\mathsf{S}=\begin{pmatrix}\left\\{\mathbf{w}_{i}\right\\}_{i=1}^{2n_{nd}}&\left\\{\mathbf{e}_{I}\right\\}_{I=1}^{n_{f}}&\left\\{\mathbf{f}_{I}\right\\}_{I=1}^{n_{f}}&\left\\{\mathbf{e}_{a}^{(\alpha)}\right\\}&\left\\{\mathbf{f}_{a}^{(\alpha)}\right\\}\end{pmatrix}$
(60)

is not necessarily symplectic, because the first set,
$\left\\{\mathbf{w}_{i}\right\\}_{i=1}^{2n_{nd}}$, need not be a symplectic
basis. If insisted upon, sGS can also be applied to the nondynamical sector,
and after the procedure $\mathsf{S}$ will indeed be symplectic (up to
reordering, as we shall see). But even without going through that step, it is
indeed the case that

$\mathsf{S}^{T}\mathsf{H}\mathsf{S}=\mathsf{H}_{D}=\mathrm{diag}\begin{pmatrix}0&0&\text{1}&\operatorname{\mathsf{\Omega}}&\operatorname{\mathsf{\Omega}}\end{pmatrix}\,,$
(61)

where $\operatorname{\mathsf{\Omega}}$ is the diagonal matrix constructed with
$\omega_{a}$ with its possible degeneracies.

Let us prove this statement, making use of all the previous analysis. To do
so, let us use a more compact notation introducing the rectangular matrices
$\mathsf{W}$, $\mathsf{E}$, $\mathsf{F}$, $\mathcal{E}$ and $\mathcal{F}$
whose columns are, respectively, the vectors $\mathbf{w}_{i}$,
$\mathbf{e}_{I}$, $\mathbf{f}_{I}$, $\mathbf{e}_{a}^{(\alpha)}$ and
$\mathbf{f}_{a}^{(\alpha)}$. Bearing in mind that $\mathsf{J}$ is invertible
and considering the action of $\mathsf{JH}$ on each of this rectangular
matrices, we have

$\displaystyle\mathsf{H}\,\mathsf{W}$
$\displaystyle=0\,\mathsf{H}\,\mathsf{E}=0\,,\mathsf{H}\,\mathsf{F}=-\mathsf{J}\,\mathsf{E}\,,$
(62) $\displaystyle\mathsf{H}\,\mathcal{E}$
$\displaystyle=\mathsf{J}\,\mathcal{F}\mathsf{\Omega}\,,\mathsf{H}\,\,\mathcal{\mathcal{F}}=-\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\,.$

The first two are simply the statement that the $\mathbf{w}_{i}$ and
$\mathbf{e}_{I}$ vectors belong to the kernel of $\mathsf{H}$. The third line
that $\mathsf{F}$ and $\mathsf{E}$ are related as elements of a Jordan chain.
The last two lines are actually a restatement of Eq. (A.2).

We now see that

$\displaystyle\mathsf{H}\,\mathsf{S}$
$\displaystyle=\mathsf{H}\begin{pmatrix}\mathsf{W}&\mathsf{E}&\mathsf{F}&\mathcal{E}&\mathcal{F}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}0&0&-\mathsf{J}\,\mathsf{E}&\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\end{pmatrix}\,.$
(63)

From here we compute Eq. (61) explicitly:

$\displaystyle\mathsf{S}^{T}\mathsf{HS}$
$\displaystyle=\begin{pmatrix}\mathsf{W}^{T}\\\ \mathsf{E}^{T}\\\
\mathsf{F}^{T}\\\ \mathcal{E}^{T}\\\
\mathcal{F}^{T}\end{pmatrix}\begin{pmatrix}0&0&-\mathsf{J}\,\mathsf{E}&\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}0&0&-\mathsf{W}^{T}\mathsf{J}\,\mathsf{E}&\mathsf{W}^{T}\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathsf{W}^{T}\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\\\
0&0&-\mathsf{E}^{T}\mathsf{J}\,\mathsf{E}&\mathsf{E}^{T}\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathsf{E}^{T}\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\\\
0&0&-\mathsf{F}^{T}\mathsf{J}\,\mathsf{E}&\mathsf{F}^{T}\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathsf{F}^{T}\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\\\
0&0&-\mathcal{E}^{T}\mathsf{J}\,\mathsf{E}&\mathcal{E}^{T}\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathcal{E}^{T}\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\\\
0&0&-\mathcal{F}^{T}\mathsf{J}\,\mathsf{E}&\mathcal{F}^{T}\mathsf{J}\mathcal{F}\mathsf{\Omega}&-\mathcal{F}^{T}\mathsf{J}\,\mathcal{E}\mathsf{\Omega}\end{pmatrix}$
(64) $\displaystyle=\begin{pmatrix}0&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&\text{1}&0&0\\\ 0&0&0&\mathsf{\Omega}&0\\\ 0&0&0&0&\mathsf{\Omega}\\\
\end{pmatrix}=\mathsf{H}_{D}\,.$ (65)

In the first row the third zero follows from Eq. (A.2). The fourth and the
fifth from that row, as well as from the second and third one, are a
consequence of Eq. (A.1). Similarly for the third zero of the fourth and fifth
rows. The third zero of the second row is from the analogue of Eq. (A.2) in
the $E$ subspace. The identity in the third element of the third row follows
from the fact that
$\left\\{\mathbf{e}_{I},\mathbf{f}_{I}\right\\}_{I=1}^{n_{f}}$ is a symplectic
basis. And analogously for the harmonic oscillator block.

Regarding whether this matrix $\mathsf{S}$ is symplectic or not, we can
compute $\mathsf{S}^{T}\mathsf{JS}$ explicitly, using the same arguments, and
it results in

$\mathsf{S}^{T}\mathsf{JS}=\begin{pmatrix}\mathsf{W}^{T}\mathsf{JW}&0&0&0&0\\\
0&0&\text{1}&0&0\\\ 0&-\text{1}&0&0&0\\\ 0&0&0&0&\text{1}\\\
0&0&0&-\text{1}&0\end{pmatrix}\,.$ (66)

We are guaranteed that $\mathsf{W}^{T}\mathsf{JF}=0$ because we have
reorganized $\tilde{W}$ into $W$ in Eq. (49). The other entries are due to sGS
in the different sectors and Eq. (A.2). We are not overly concerned with the
presence of the antisymmetric entry $\mathsf{W}^{T}\mathsf{JW}$, because it
pertains to the nondynamical sector which will be eliminated from any
posterior analysis. To make this statement abundantly clear, consider the
classical equations of motion for the initial variables, Eq. (1), that we
restate here:

$\displaystyle\dot{\operatorname{\boldsymbol{x}}}=\operatorname{\mathsf{J}}\operatorname{\mathsf{H}}\operatorname{\boldsymbol{x}}\,.$

Let us define now $\operatorname{\boldsymbol{y}}$ by
$\operatorname{\boldsymbol{x}}=\operatorname{\mathsf{S}}\operatorname{\boldsymbol{y}}$,
with $\operatorname{\mathsf{S}}$ defined in Eq. (60), i.e.

$\mathsf{S}=\begin{pmatrix}\mathsf{W}&\mathsf{E}&\mathsf{F}&\mathcal{E}&\mathcal{F}\end{pmatrix}\,.$
(67)

From Eq. (66) one can determine its inverse,

$\mathsf{S}^{-1}=\begin{pmatrix}\left(\mathsf{W}^{T}\mathsf{JW}\right)^{-1}&0&0&0&0\\\
0&0&-\text{1}&0&0\\\ 0&\text{1}&0&0&0\\\ 0&0&0&0&-\text{1}\\\
0&0&0&\text{1}&0\end{pmatrix}\mathsf{S}^{T}\mathsf{J}\,.$ (68)

Thus

$\displaystyle\dot{\operatorname{\boldsymbol{y}}}$
$\displaystyle=\mathsf{S}^{-1}\mathsf{JH}\,\mathsf{S}\operatorname{\boldsymbol{y}}$
$\displaystyle=\begin{pmatrix}0&0&0&0&0\\\ 0&0&\text{1}&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&\mathsf{\Omega}\\\
0&0&0&-\mathsf{\Omega}&0\end{pmatrix}\operatorname{\boldsymbol{y}}\,,$ (69)

as promised throughout the analysis. The possibly nontrivial
$\mathsf{W}^{T}\mathsf{JW}$ plays no role, as it pertains to the nondynamical
sector.

##  Appendix B Landau quantization

Working on the symmetric gauge, we can write the Hamiltonian for the two
linearly-coupled charged particles in a magnetic field ($z$ direction) as
$H=\frac{1}{2}\operatorname{\boldsymbol{X}}_{z}^{T}\operatorname{\mathsf{H}}_{z}\operatorname{\boldsymbol{X}}_{z}+\frac{1}{2}\operatorname{\boldsymbol{X}}_{xy}^{T}\operatorname{\mathsf{H}}_{xy}\operatorname{\boldsymbol{X}}_{xy}$
where

$\displaystyle\operatorname{\mathsf{H}}_{z}$
$\displaystyle=\begin{pmatrix}k&-k&0&0\\\ -k&k&0&0\\\ 0&0&\frac{1}{m}&0\\\
0&0&0&\frac{1}{m}\end{pmatrix},$ (70)
$\displaystyle\operatorname{\mathsf{H}}_{xy}$
$\displaystyle=\begin{pmatrix}\tilde{k}&-k&0&0&0&0&\frac{-g}{m}&0\\\
-k&\tilde{k}&0&0&0&0&0&\frac{-g}{m}\\\ 0&0&\tilde{k}&-k&\frac{g}{m}&0&0&0\\\
0&0&-k&\tilde{k}&0&\frac{g}{m}&0&0\\\ 0&0&\frac{g}{m}&0&\frac{1}{m}&0&0&0\\\
0&0&0&\frac{g}{m}&0&\frac{1}{m}&0&0\\\ \frac{-g}{m}&0&0&0&0&0&\frac{1}{m}&0\\\
0&\frac{-g}{m}&0&0&0&0&0&\frac{1}{m}\end{pmatrix}.$ (71)

with $g=eB/2$, $\tilde{k}=k+\frac{g^{2}}{m}$,
$\operatorname{\boldsymbol{X}}_{z}=(z_{1},z_{2},p_{z,1},p_{z,2})^{T}$ and
$\operatorname{\boldsymbol{X}}_{xy}=(x_{1},x_{2},y_{1},y_{2},p_{x,1},p_{x,2},p_{y,1},p_{y,2})^{T}$.
The harmonic oscillator in the $z$ direction has frequency
$\Omega_{z}=\sqrt{2k/m}$. Now, let us show how to perform a blindfold
diagonalization of the dynamics in the $x$–$y$. By first making a rescaling of
the variables to simplify the parameters, the configuration space variables in
the plane are modified as $q_{i}\rightarrow\frac{q_{i}}{\sqrt{m\Omega_{c}}}$
and the momenta accordingly $p_{i}\rightarrow\sqrt{m\Omega_{c}}p_{i}$, where
$\Omega=g/2m=\Omega_{c}/4$, with $\Omega_{c}=eB/m$ the cyclotron frequency.
This Hamiltonian can thus be written as

$\displaystyle\operatorname{\mathsf{H}}_{xy}$
$\displaystyle=\Omega_{c}\begin{pmatrix}\chi^{2}+\frac{1}{4}&-\chi^{2}&0&0&0&0&-\frac{1}{2}&0\\\
-\chi^{2}&\chi^{2}+\frac{1}{4}&0&0&0&0&0&-\frac{1}{2}\\\
0&0&\chi^{2}+\frac{1}{4}&-\chi^{2}&\frac{1}{2}&0&0&0\\\
0&0&-\chi^{2}&\chi^{2}+\frac{1}{4}&0&\frac{1}{2}&0&0\\\
0&0&\frac{1}{2}&0&1&0&0&0\\\ 0&0&0&\frac{1}{2}&0&1&0&0\\\
-\frac{1}{2}&0&0&0&0&0&1&0\\\ 0&-\frac{1}{2}&0&0&0&0&0&1\end{pmatrix},$ (72)

where we have defined the unitless parameter $\chi=\Omega_{z}/\Omega\sqrt{2}$.
The normal frequencies of this sector can be found by diagonalizing
$\operatorname{\mathsf{J}}_{xy}\operatorname{\mathsf{H}}_{xy}$, and correspond
to the cyclotron motion ($\Omega$), and one symmetric and one antisymmetric
modes $\left(\Omega\sqrt{\sqrt{1+8\chi^{2}}\pm(1+4\chi^{2})/2}\right)$.
Following the algorithm, we can obtain systematically symplectic
transformations that diagonalize Hamiltonian (72). For example, at the
degeneracy point $\chi=1$,

$\displaystyle\operatorname{\mathsf{S}}|_{\chi=1}=\begin{pmatrix}\operatorname{\mathsf{W}}&\mathcal{E}&\mathcal{F}\end{pmatrix}=\left(\begin{array}[]{cccccccc}\frac{1}{\sqrt{2}}&0&-\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}&0&0&0\\\
\frac{1}{\sqrt{2}}&0&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{6}}&-\frac{1}{\sqrt{6}}&0&0&0\\\
0&-\frac{1}{\sqrt{2}}&0&0&0&-\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}\\\
0&-\frac{1}{\sqrt{2}}&0&0&0&-\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{6}}&-\frac{1}{\sqrt{6}}\\\
0&\frac{1}{2\sqrt{2}}&0&0&0&-\frac{1}{2\sqrt{2}}&-\sqrt{\frac{3}{8}}&\sqrt{\frac{3}{8}}\\\
0&\frac{1}{2\sqrt{2}}&0&0&0&-\frac{1}{2\sqrt{2}}&\sqrt{\frac{3}{8}}&-\sqrt{\frac{3}{8}}\\\
\frac{1}{2\sqrt{2}}&0&\frac{1}{2\sqrt{2}}&-\sqrt{\frac{3}{8}}&-\sqrt{\frac{3}{8}}&0&0&0\\\
\frac{1}{2\sqrt{2}}&0&\frac{1}{2\sqrt{2}}&\sqrt{\frac{3}{8}}&\sqrt{\frac{3}{8}}&0&0&0\\\
\end{array}\right),$ (81)

such that the diagonalized Hamiltonian becomes

$\displaystyle\operatorname{\mathsf{H}}_{D}|_{\chi=1}=\Omega\left(\begin{array}[]{cccccccc}0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&2&0&0&0\\\
0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&2\\\ \end{array}\right).$
(90)

##  Appendix C Conserved linear and quadratic quantities

In any study of the Landau problem one will immediately encounter a
description of the linear and quadratic conserved quantities. For the one
particle Landau problem, the latter is the angular momentum, which reduces to
the energy minus a constant associated with the nondynamical sector when on-
shell. The characterization and identification of all linear and quadratic
conserved quantities for a linear system is, as that example shows, useful and
relevant, most importantly with a view towards the nonlinear case. In this
section we observe that the application of Williamson’s theorem to positive
semidefinite Hamiltonians allows us to give a complete set of independent
linear and quadratic conserved quantities.

The idea is the following: as there is a canonical symplectically diagonalized
form of the Hamiltonian matrix, characterize the invariants in that canonical
form, and undo the symplectic transformation to write the full set in the
initial variables. We shall carry out the analysis using the classical
equations of motion, because for operators linear and quadratic in the initial
variables the classical equations match Heisenberg’s equations of motion.

The classical equations of motion for functions on phase space without
explicit time dependence are

$\dot{f}=\left\\{f,H\right\\}\,,$ (91)

which for the case under consideration become

$\dot{f}=\mathsf{J}_{\alpha\beta}\partial_{\alpha}f\,\mathsf{H}_{\beta\gamma}x^{\gamma}\,.$
(92)

Alternatively,

$\dot{f}=-\operatorname{\boldsymbol{x}}^{T}\mathsf{H}\mathsf{J}\nabla f\,.$
(93)

If $f$ is a linear function on phase space,
$f(\operatorname{\boldsymbol{x}})=\operatorname{\boldsymbol{x}}^{T}\mathbf{f}+d$,
with $\mathbf{f}$ a column vector, and thus the condition of being a conserved
quantity is that for all points of phase space $\operatorname{\boldsymbol{x}}$
one has

$\operatorname{\boldsymbol{x}}^{T}\mathsf{H}\mathsf{J}\mathbf{f}=0\,.$ (94)

It follows that $\mathsf{J}\mathbf{f}$ must belong to the kernel of
$\mathsf{H}$. Thus we identify the set of linear conserved quantities with
$\mathsf{J}\mathrm{ker}\left[\mathsf{H}\right]$. If the reader is wondering
about the need for the symplectic matrix, please observe, first, that it is
invertible, so we indeed have that a linear space of linear conserved
quantities that is isomorphic to the kernel of $\mathsf{H}$ and, second, that
we have seen that the $\mathrm{ker}\left[\mathsf{H}\right]$ consists of the
nondynamical sector of phase space and the position like part of the phase
space of free particles, since it is the $E$ space. The conserved quantity for
a free particle is its momentum, i.e. $F$ space in the notation above, and
this the result from applying $\mathsf{J}$ to $E$, essentially.

Coming now to quadratic conserved quantities, we write the symmetric matrix
associated with one such $f$ as $\mathsf{B}=\mathsf{B}^{T}$ (to avoid
confusion with a previous use of $\mathsf{F}$), as

$f(\operatorname{\boldsymbol{x}})=\frac{1}{2}\operatorname{\boldsymbol{x}}^{T}\mathsf{B}\operatorname{\boldsymbol{x}}\,.$
(95)

Now Eq. (93) becomes

$\dot{f}=-\operatorname{\boldsymbol{x}}^{T}\mathsf{H}\mathsf{J}\mathsf{B}\operatorname{\boldsymbol{x}}\,.$
(96)

For $f$ to be a conserved quantity we need the right hand side to vanish for
all points on phase space. Therefore the symmetric part of the matrix must be
zero, that is

$\mathsf{HJB}-\mathsf{BJH}=0\,,$ (97)

which is the condition we must now study, given a Hamiltonian matrix
$\mathsf{H}$. In fact, premultiplying with the symplectic matrix $\mathsf{J}$
one obtains the satisfactory physical interpretation that the evolution in
time controlled by $\mathsf{H}$ commutes with that generated by
$\operatorname{\mathsf{B}}$. Quantum mechanically, the associated Hamiltonians
must conmute, and in terms of the matrices we have

$\left[\mathsf{JH},\mathsf{JB}\right]=0\,.$ (98)

Thus the problem is reduced to the study of the commutant of the normal form
of $\mathsf{JH}$. If $\mathsf{H}$ is positive, the normal form is
$\mathrm{diag}\left(i\Omega,-i\Omega\right)$, and is obtained from
$\mathsf{JH}$ by conjugation. If there are no degeneracies, the commutant is
given by diagonal matrices. A further condition, however, must be met, namely
that $\mathsf{B}$ be real and symmetric. Given $\mathbf{v}_{a}$ eigenvector of
$\mathsf{JH}$ with eigenvalue $-i\omega_{a}$, it must be the case that it be
an eigenvector of $\mathsf{JB}$, with eigenvalue $\beta_{a}$. It follows that
$\mathbf{v}_{a}^{*}$ must also be an eigenvector of $\mathsf{JB}$, now with
eigenvalue $\beta^{*}_{a}$. In terms of the real and imaginary parts of
$\mathbf{v}_{a}$, $\mathbf{v}_{a}=\mathbf{e}_{a}+i\mathbf{f}_{a}$, and
denoting the real and imaginary parts of $\beta_{a}$ as
$\beta_{a}=\mu_{a}+i\nu_{a}$, we have that in the real symplectic basis
$\left\\{\mathbf{e}_{a},\mathbf{f}_{a}\right\\}$ (or proportional to
symplectic) $\mathsf{JB}$ has the matrix

$\mathsf{JB}\Big{|}_{a}=\begin{pmatrix}\mu_{a}&\nu_{a}\\\
-\nu_{a}&\mu_{a}\end{pmatrix}\,,$ (99)

and thus, in the same basis

$\operatorname{\mathsf{B}}_{a}=\begin{pmatrix}\nu_{a}&-\mu_{a}\\\
\mu_{a}&\nu_{a}\end{pmatrix}\,.$ (100)

Since $\operatorname{\mathsf{B}}$ is symmetric, this means that (as expected)
$\mu_{a}=0$.

Let us now examine the positive definite degenerate case, and use the notation
of Eq. (55) in Section A.2. The matrix $\operatorname{\mathsf{B}}$ restricted
to the $\omega_{a}$ eigenspace in the basis
$\left\\{e_{a}^{(\alpha)}\right\\}\cup\left\\{f_{a}^{(\alpha)}\right\\}$ will
be of the form

$\operatorname{\mathsf{B}}_{a}=\begin{pmatrix}\operatorname{\mathsf{N}}_{a}&-\operatorname{\mathsf{M}}_{a}\\\
\operatorname{\mathsf{M}}_{a}&\operatorname{\mathsf{N}}_{a}\end{pmatrix}$
(101)

with $\mathsf{M}_{a}$ and $\mathsf{N}_{a}$ matrices. As $\mathsf{B}$ is
symmetric, this means that $\mathsf{M}_{a}$ must be antisymmetric and
$\mathsf{N}_{a}$ symmetric.

For example, two harmonic oscillators with different frequencies have no
quadratic invariants other than their separate energy functions and linear
superpositions thereof. On the other hand the isotropic harmonic oscillator on
the plane presents a basis of four quadratic invariants: the two energy
functions $\left(x_{i}^{2}+p_{i}^{2}\right)$, the canonical angular momentum
$x_{1}p_{2}-x_{2}p_{1}$, and a final $x_{1}x_{2}+p_{1}p_{2}$ (or alternatively
$\left(x_{1}+x_{2}\right)^{2}+\left(p_{1}+p_{2}\right)^{2}$).

Before examining the general positive semidefinite case, let us illustrate the
difference with the one particle Landau problem. After some symplectic
transformations, the Hamiltonian matrix is proportional to

$\mathsf{H}=\begin{pmatrix}1&0&0&1\\\ 0&1&-1&0\\\ 0&-1&1&0\\\
1&0&0&1\end{pmatrix}\,,$ (102)

with one non dynamical degree of freedom and one harmonic oscillator. Thus the
canonical form is proportional to

$\tilde{\mathsf{H}}_{D}=\mathrm{diag}\begin{pmatrix}0&1&0&1\end{pmatrix}\,.$
(103)

On inserting this form in Eq. (98) one obtains that in the canonical basis the
quadratic invariants are the trivial ones involving the nondynamical sector
(three independent ones, $\xi^{2}$, $p_{\xi}^{2}$, and $\xi p_{\xi}$, with
$\xi,p_{\xi}$ canonical coordinates of the nondynamical sector), and the
energy invariant of the oscillator, $x^{2}+p_{x}^{2}$.

In this example one sees a feature of the general case: the harmonic
oscillator sector will not couple with the nondynamical and free particle
sectors in invariant quadratic forms. On the other hand, as can readily be
seen from the trivial example of one free particle and one nondynamical, these
do couple into invariants. These invariants, however, are in that sector not
independent from the linear invariants, that generate all possible invariants.
In summary, the only independent quadratic invariants are to be found in the
harmonic oscillator sector.

We can rephrase the presentation above in more abstract terms as follows: (i)
$\mathsf{JH}$ belongs to the Lie algebra for the symplectic group; (2) By
using symplectic transformations on $\mathsf{H}$, $\mathsf{JH}$ can be written
as being upper triangular; (3) the quadratic conserved quantities are
determined by the commutant of $\mathsf{JH}$ in the Lie algebra.

###  C.1 Gauge invariance and symplectic diagonalization

The physical content of Hamiltonians such as the minimal coupling example of
Eq. (6) is independent of the gauge in which it is written. However, there is
a limited set of gauges in which the Hamiltonian is a quadratic function on
phase space. Namely, the symmetric gauge and those obtained from the symmetric
gauge vector potential and adding the gradient of a quadratic function,

$\displaystyle\boldsymbol{A}=\boldsymbol{A}_{S}+\nabla\xi$ (104)

with

$\displaystyle\xi=\frac{1}{2}\operatorname{\boldsymbol{x}}^{T}\mathsf{\Xi}\operatorname{\boldsymbol{x}}\,.$
(105)

In fact, this induces a symplectic transformation.

As the symplectic diagonalization approach is precisely a tool to identify the
evolution of a system in terms of its symplectic invariants, changes of gauge
such as these have no impact on the process and the result.

##  Appendix D Linear reciprocal networks: counting of degrees of freedom

The Hamiltonian in Eq. (10) is evidently positive semidefinite. The
corresponding Hamiltonian matrix is

$\displaystyle\mathsf{H}=\begin{pmatrix}0&0&0&0\\\
0&\mathsf{D}^{T}\mathsf{L}^{-1}\mathsf{D}&-\mathsf{D}^{T}\mathsf{L}^{-1}&0\\\
0&-\mathsf{L}^{-1}\mathsf{D}&\operatorname{\mathsf{L}}^{-1}&0\\\
0&0&0&\operatorname{\mathsf{C}}^{-1}\end{pmatrix}\,.$ (106)

$\mathsf{C}$ and $\mathsf{L}$ are positive definite matrices defined on
vectors with $N_{C}$ and $M_{L}$ components respectively. The adjacency matrix
$\mathsf{D}$ has dimension $N_{C}\times M_{L}$ and is possibly not of maximum
rank. Denote with $\tilde{\operatorname{\boldsymbol{e}}}$ a generic element of
the kernel of $\mathsf{D}$ and with
$\tilde{\operatorname{\boldsymbol{e}}}_{T}$ a generic element of
$\mathrm{ker}\left[\mathsf{D}^{T}\right]$. Similarly represent a generic
element of $\mathbb{R}^{N_{C}}$ as $\operatorname{\boldsymbol{\phi}}$, and a
generic element of $\mathbb{R}^{M_{L}}$ as $\boldsymbol{q}$.

In order to assess the number of nondynamical, free particle and harmonic
oscillator degrees of freedom we have to compute the dimensions of
$K_{1}=\mathrm{ker}\left[\mathsf{JH}\right]$ and
$K_{2}=\mathrm{ker}\left[\left(\mathsf{JH}\right)^{2}\right]$. In this case it
is easy to compute an explicit form for the generic elements of $K_{1}$ and
$K_{2}$, namely

$\displaystyle\begin{pmatrix}\operatorname{\boldsymbol{q}}\\\
\operatorname{\boldsymbol{\phi}}\\\
\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\phi}}\\\ 0\end{pmatrix}\in
K_{1},\qquad\text{and}\qquad\begin{pmatrix}\operatorname{\boldsymbol{q}}\\\
\operatorname{\boldsymbol{\phi}}\\\
\operatorname{\mathsf{D}}\operatorname{\boldsymbol{\phi}}+\operatorname{\mathsf{L}}\tilde{\operatorname{\boldsymbol{e}}}_{T}\\\
\operatorname{\mathsf{C}}\tilde{\operatorname{\boldsymbol{e}}}\end{pmatrix}\in
K_{2}\,.$ (107)

We read directly the dimensions and the number of free particles,

$\displaystyle\mathrm{dim}\left(K_{1}\right)$ $\displaystyle=N_{C}+M_{L}\,,$
$\displaystyle\mathrm{dim}\left(K_{2}\right)$
$\displaystyle=N_{C}+M_{L}+n_{f}\,,$ $\displaystyle n_{f}$
$\displaystyle=\mathrm{dim}\left(\mathrm{ker}\left[\mathsf{D}\right]\right)+\mathrm{dim}\left(\mathrm{ker}\left[\mathsf{D}^{T}\right]\right)\,.$
(108)

It is easy to build $E$ and $\tilde{F}$ spaces, in the notation of Appendix A,
as being generated, respectively by vectors of the form

$\displaystyle\begin{pmatrix}\tilde{\operatorname{\boldsymbol{e}}}_{T}\\\
\tilde{\operatorname{\boldsymbol{e}}}\\\ 0\\\ 0\end{pmatrix}\in
E\qquad\mathrm{and}\quad\begin{pmatrix}0\\\ 0\\\
\mathsf{L}\tilde{\operatorname{\boldsymbol{e}}}_{T}\\\
\mathsf{C}\tilde{\operatorname{\boldsymbol{e}}}\end{pmatrix}\in\tilde{F}\,.$
(109)

Since the dimension of $K_{1}$ is the number of free particles plus twice the
number of nondynamical variables, we express the number of nondynamical
variables for this case of linear reciprocal networks as

$\displaystyle n_{nd}$
$\displaystyle=\frac{1}{2}\left(\mathrm{dim}\left(K_{1}\right)-n_{f}\right)$
$\displaystyle=\frac{1}{2}\left[N_{C}+M_{L}-n_{f}\right]\,.$ (110)

Passing now to the harmonic oscillator sector, notice that
$n_{ho}+n_{f}+n_{nd}$ is the total number of configuration space variables, in
our case $N_{C}+M_{L}$, whence we compute

$\displaystyle n_{ho}$ $\displaystyle=N_{C}+M_{L}-n_{f}-n_{nd}$
$\displaystyle=\frac{1}{2}\left[N_{C}+M_{L}-n_{f}\right]=n_{nd}\,,$ (111)

so the number of harmonic oscillators is the same as that of nondynamical
variables in this case of linear reciprocal networks, as stated in the main
text.

###  D.1 Example

In the text the example of Fig. 3 (b) has been put forward. In this situation,
with $\mathsf{D}=\begin{pmatrix}1&1\end{pmatrix}$, we see that
$\mathrm{ker}\left[\mathsf{D}^{T}\right]=\emptyset$, while
$\mathrm{ker}\left[\mathsf{D}\right]$ is one-dimensional, so $n_{f}=1$. We
have $N_{C}=2$, $M_{L}=1$, whence $n_{nd}=n_{ho}=1$. Applying Eqs. (47) and
(48), we obtain for the free particle sector

$\displaystyle\mathbf{e}_{f}=\frac{1}{\sqrt{C_{1}+C_{2}}}\begin{pmatrix}0\\\
1\\\ -1\\\ 0\\\ 0\\\
0\end{pmatrix}\quad\mathrm{and}\quad\mathbf{f}_{f}=\frac{1}{\sqrt{C_{1}+C_{2}}}\begin{pmatrix}0\\\
0\\\ 0\\\ 0\\\ C_{1}\\\ -C_{2}\end{pmatrix}\,.$ (112)

From inspection of $K_{1}$ in (107) and comparison (112) we construct
$\tilde{W}$, as defined in Eq. (45), and by symplectic orthogonalization, Eq.
(49) we obtain

$\displaystyle\mathbf{w}_{1}=\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\\\ 0\\\
0\end{pmatrix}\,,\qquad\mathbf{w}_{2}=\frac{1}{C_{1}+C_{2}}\begin{pmatrix}0\\\
C_{2}\\\ C_{1}\\\ C_{1}+C_{2}\\\ 0\\\ 0\end{pmatrix}\,.$ (113)

The full symplectic basis is completed with the harmonic oscillator sector,

$\displaystyle\mathbf{e}=\frac{1}{\sqrt{\Omega L}}\begin{pmatrix}1\\\ 0\\\
0\\\ 1\\\
1\end{pmatrix}\,,\quad\mathbf{f}=\frac{-1}{\Omega^{3/4}L^{1/2}}\begin{pmatrix}0\\\
1/C_{1}\\\ 1/C_{2}\\\ 0\\\ 0\\\ 0\end{pmatrix},\,$ (114)

where $\Omega=1/\sqrt{L(C_{1}+C_{2})}$.

##  Appendix E Example circuit: JJs coupled to a NR black-box admittance

Here we provide details for the quantization of the nonreciprocal two-port
admittance capacitively coupled to Josephson junctions in Fig. 7. After
writing systematically the Lagrangian (IV) and Hamiltonian (IV), let us
perform a triangular transformation (shift) of the internal charge degrees of
freedom with the Josephson fluxes,

$\displaystyle\tilde{\operatorname{\boldsymbol{Q}}}$
$\displaystyle=\operatorname{\boldsymbol{Q}}-2\operatorname{\mathsf{Z}}^{-1}\operatorname{\mathsf{D}}_{J}\operatorname{\boldsymbol{\Phi}}_{J},\quad\tilde{\operatorname{\boldsymbol{P}}}=\operatorname{\boldsymbol{P}},$
(115) $\displaystyle\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}$
$\displaystyle=\operatorname{\boldsymbol{\Phi}}_{J},\quad\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}=\operatorname{\boldsymbol{\Pi}}_{J}-2\operatorname{\mathsf{D}}_{J}^{T}\operatorname{\mathsf{Z}}^{-1}\operatorname{\boldsymbol{P}},$
(116)

while not modiying the coupling capacitor variables
$\tilde{\operatorname{\boldsymbol{\Phi}}}_{c}=\operatorname{\boldsymbol{\Phi}}_{c}$
and
$\tilde{\operatorname{\boldsymbol{\Pi}}}_{c}=\operatorname{\boldsymbol{\Pi}}_{c}$.
The Hamiltonian can be rewritten then

$\displaystyle H=$
$\displaystyle\,\frac{1}{2}(\tilde{\operatorname{\boldsymbol{P}}}-\tfrac{1}{2}\operatorname{\mathsf{Z}}\tilde{\operatorname{\boldsymbol{Q}}}+\operatorname{\mathsf{D}}_{c}\tilde{\operatorname{\boldsymbol{\Phi}}}_{c})^{T}\operatorname{\mathsf{L}}^{-1}(\tilde{\operatorname{\boldsymbol{P}}}-\tfrac{1}{2}\operatorname{\mathsf{Z}}\tilde{\operatorname{\boldsymbol{Q}}}+\operatorname{\mathsf{D}}_{c}\tilde{\operatorname{\boldsymbol{\Phi}}}_{c})$
$\displaystyle+\frac{1}{2}(\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}+2\operatorname{\mathsf{D}}_{J}^{T}\operatorname{\mathsf{Z}}^{-1}\tilde{\operatorname{\boldsymbol{P}}})^{T}\operatorname{\mathsf{C}}_{J}^{-1}(\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}+2\operatorname{\mathsf{D}}_{J}^{T}\operatorname{\mathsf{Z}}^{-1}\tilde{\operatorname{\boldsymbol{P}}})+\frac{1}{2}\tilde{\operatorname{\boldsymbol{\Pi}}}_{c}\operatorname{\mathsf{C}}_{c}^{-1}\tilde{\operatorname{\boldsymbol{\Pi}}}_{c}+U(\tilde{\operatorname{\boldsymbol{\Phi}}}_{J})$
$\displaystyle=$
$\displaystyle\tilde{\operatorname{\boldsymbol{X}}}^{T}\operatorname{\mathsf{H}}\tilde{\operatorname{\boldsymbol{X}}}+U(\tilde{\operatorname{\boldsymbol{\Phi}}}_{J})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\tilde{\operatorname{\boldsymbol{Y}}}^{T}\operatorname{\mathsf{H}}_{\text{lin}}\tilde{\operatorname{\boldsymbol{Y}}}+\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}^{T}\operatorname{\mathsf{C}}_{J}^{-1}\left(\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}+4\operatorname{\mathsf{D}}_{J}^{T}\operatorname{\mathsf{Z}}^{-1}\tilde{\operatorname{\boldsymbol{P}}}\right)\right]+U(\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}),$
(117)

where we have defined vectors
$\tilde{\operatorname{\boldsymbol{Y}}}^{T}=(\tilde{\operatorname{\boldsymbol{Q}}}^{T},\tilde{\operatorname{\boldsymbol{\Phi}}}_{c}^{T},\tilde{\operatorname{\boldsymbol{P}}}^{T},\tilde{\operatorname{\boldsymbol{\Pi}}}_{c}^{T})$,
and
$\tilde{\operatorname{\boldsymbol{X}}}^{T}=(\tilde{\operatorname{\boldsymbol{Q}}}^{T},\tilde{\operatorname{\boldsymbol{\Phi}}}_{c}^{T},\tilde{\operatorname{\boldsymbol{P}}}^{T},\tilde{\operatorname{\boldsymbol{\Pi}}}_{c}^{T},\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}^{T},\tilde{\operatorname{\boldsymbol{\Pi}}}_{J}^{T})$,
the nonlinear potential
$U(\operatorname{\boldsymbol{\Phi}}_{J})=-\sum_{n}E_{Jn}\cos(\varphi_{Jn})$,
with the matrices

$\displaystyle\operatorname{\mathsf{H}}=$
$\displaystyle\begin{pmatrix}\operatorname{\mathsf{H}}_{\text{lin}}&\operatorname{\mathsf{H}}_{\text{int}}\\\
\operatorname{\mathsf{H}}_{\text{int}}^{T}&\operatorname{\mathsf{H}}_{J}\end{pmatrix},$
(118) $\displaystyle\operatorname{\mathsf{H}}_{\text{lin}}=$
$\displaystyle\begin{pmatrix}\frac{\operatorname{\mathsf{Z}}^{T}\operatorname{\mathsf{L}}^{-1}\operatorname{\mathsf{Z}}}{4}&\frac{-\operatorname{\mathsf{Z}}^{T}\operatorname{\mathsf{L}}^{-1}\operatorname{\mathsf{D}}_{c}}{2}&\frac{-\operatorname{\mathsf{Z}}^{T}\operatorname{\mathsf{L}}^{-1}}{2}&0\\\
\frac{-\operatorname{\mathsf{D}}_{c}^{T}\operatorname{\mathsf{L}}^{-1}\tilde{\operatorname{\mathsf{Z}}}}{2}&\operatorname{\mathsf{D}}_{c}^{T}\operatorname{\mathsf{L}}^{-1}\operatorname{\mathsf{D}}_{c}&\operatorname{\mathsf{D}}_{c}^{T}\operatorname{\mathsf{L}}^{-1}&0\\\
\frac{-\operatorname{\mathsf{L}}^{-1}\operatorname{\mathsf{Z}}}{2}&\operatorname{\mathsf{L}}^{-1}\operatorname{\mathsf{D}}_{c}&\tilde{\operatorname{\mathsf{L}}}^{-1}&0\\\
0&0&0&\operatorname{\mathsf{C}}_{c}^{-1}\end{pmatrix},$ (119)

where
$\tilde{\operatorname{\mathsf{L}}}^{-1}=\operatorname{\mathsf{L}}^{-1}-4\operatorname{\mathsf{Z}}^{-1}\operatorname{\mathsf{D}}_{J}\operatorname{\mathsf{C}}_{J}^{-1}\operatorname{\mathsf{D}}_{J}^{T}\operatorname{\mathsf{Z}}^{-1}$.
Note that the Josephson variables
$\tilde{\operatorname{\boldsymbol{\Phi}}}_{J}$ only appear in the nonlinear
potential, and as previously mentioned, they should be modeled with compact
variables $\varphi_{Jn}=(2\pi)\frac{\Phi_{Jn}}{\Phi_{0}}\in\mathcal{S}^{1}$.
We then perform a rescaling transformations (another trivial transformation)
so that all classical coordinates have the same units,

$\displaystyle\bar{\operatorname{\boldsymbol{Q}}}$
$\displaystyle=\sqrt{R}\tilde{\operatorname{\boldsymbol{Q}}},\qquad\bar{\operatorname{\boldsymbol{P}}}=\tilde{\operatorname{\boldsymbol{P}}}/\sqrt{R},$
(120) $\displaystyle\bar{\operatorname{\boldsymbol{\Phi}}}_{c/J}$
$\displaystyle=\operatorname{\mathsf{C}}_{c/J}^{-\frac{1}{4}}\operatorname{\mathsf{L}}^{\frac{1}{4}}\tilde{\operatorname{\boldsymbol{\Phi}}}_{c/J},\qquad\bar{\operatorname{\boldsymbol{\Pi}}}_{c/J}=\operatorname{\mathsf{C}}_{c/J}^{\frac{1}{4}}\operatorname{\mathsf{L}}^{-\frac{1}{4}}\tilde{\operatorname{\boldsymbol{\Pi}}}_{c/J}.$

The matrix of Hamiltonian (IV) can now be written, for the specific example at
hand, Eq. (22), in terms of three frequencies, $\Omega_{J}=1/\sqrt{C_{J}L}$,
$\Omega_{c}=1/\sqrt{C_{c}L}$, and $\Omega=R/L$. Furthermore, without loss of
generality, we set the transformer turn-ratios matrix to
$n_{11}=n_{12}=n_{22}=1$ and $n_{21}=0$.

###  E.1 No Josephson junctions

Let us begin the analysis of the circuit by setting to zero the nonlinear
potential, i.e., $E_{Jn}\rightarrow 0$ such that
$U(\operatorname{\boldsymbol{\Phi}}_{J})=0$ in Hamiltonian (IV). In the
example, Eq. (22), after the transformations of Eq. (120) and using the
transformer turn-ratios mentioned above, the Hamiltonian in this case becomes
$2H_{U=0}=\bar{\operatorname{\boldsymbol{X}}}^{T}\operatorname{\mathsf{H}}\bar{\operatorname{\boldsymbol{X}}}$,
with

$\displaystyle\operatorname{\mathsf{H}}=\left(\begin{array}[]{cccccccccccc}\frac{\Omega}{4}&0&\frac{\sqrt{\Omega\Omega_{c}}}{2}&\frac{\sqrt{\Omega\Omega_{c}}}{2}&0&-\frac{\Omega}{2}&0&0&0&0&0&0\\\
0&\frac{\Omega}{4}&-\frac{\sqrt{\Omega\Omega_{c}}}{2}&0&\frac{\Omega}{2}&0&0&0&0&0&0&0\\\
\frac{\sqrt{\Omega\Omega_{c}}}{2}&-\frac{\sqrt{\Omega\Omega_{c}}}{2}&2\Omega_{c}&\Omega_{c}&-\sqrt{\Omega\Omega_{c}}&-\sqrt{\Omega\Omega_{c}}&0&0&0&0&0&0\\\
\frac{\sqrt{\Omega\Omega_{c}}}{2}&0&\Omega_{c}&\Omega_{c}&0&-\sqrt{\Omega\Omega_{c}}&0&0&0&0&0&0\\\
0&\frac{\Omega}{2}&-\sqrt{\Omega\Omega_{c}}&0&\frac{8\Omega_{J}^{2}}{\Omega}+\Omega&-\frac{4\Omega_{J}^{2}}{\Omega}&0&0&0&0&\frac{2\Omega_{J}^{3/2}}{\sqrt{\Omega}}&\frac{2\Omega_{J}^{3/2}}{\sqrt{\Omega}}\\\
-\frac{\Omega}{2}&0&-\sqrt{\Omega\Omega_{c}}&-\sqrt{\Omega\Omega_{c}}&-\frac{4\Omega_{J}^{2}}{\Omega}&\frac{4\Omega_{J}^{2}}{\Omega}+\Omega&0&0&0&0&-\frac{2\Omega_{J}^{3/2}}{\sqrt{\Omega}}&0\\\
0&0&0&0&0&0&\Omega_{c}&0&0&0&0&0\\\ 0&0&0&0&0&0&0&\Omega_{c}&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&\frac{2\Omega_{J}^{3/2}}{\sqrt{\Omega}}&-\frac{2\Omega_{J}^{3/2}}{\sqrt{\Omega}}&0&0&0&0&\Omega_{J}&0\\\
0&0&0&0&\frac{2\Omega_{J}^{3/2}}{\sqrt{\Omega}}&0&0&0&0&0&0&\Omega_{J}\\\
\end{array}\right)$ (133)

and
$\bar{\operatorname{\boldsymbol{X}}}^{T}=(\bar{\operatorname{\boldsymbol{Q}}}^{T},\bar{\operatorname{\boldsymbol{\Phi}}}_{c}^{T},\bar{\operatorname{\boldsymbol{P}}}^{T},\bar{\operatorname{\boldsymbol{\Pi}}}_{c}^{T},\bar{\operatorname{\boldsymbol{\Phi}}}_{J}^{T},\bar{\operatorname{\boldsymbol{\Pi}}}_{J}^{T})$.
Solving the eigenvalue problem for $\operatorname{\mathsf{JH}}$ we can find
the symplectic transformation. For the homogeneous case of frequencies
$\Omega=\Omega_{c}=\Omega_{J}=1$, say, we have

$\displaystyle\operatorname{\mathsf{S}}|_{\Omega_{\alpha}=1}\stackrel{{\scriptstyle\text{num.}}}{{=}}$
$\displaystyle\left(\begin{array}[]{cccccccccccc}-1&1&-1&-1&\sqrt{2}&\frac{1}{\sqrt{2}}&0&0&1.38612&0.435045&0.&0.\\\
1&0&0&1&-\sqrt{2}&0&0&0&-0.634047&0.634047&0.800274&-0.251173\\\
0&0&0&\frac{1}{2}&-\frac{1}{\sqrt{2}}&0&0&0&0.317023&-0.317023&-0.125587&0.400137\\\
0&0&\frac{1}{2}&0&0&-\frac{1}{2\sqrt{2}}&0&0&0.217523&0.693058&0.&0.\\\
-\frac{1}{2}&0&0&0&0&0&0&0&0.&0.&0.27455&0.27455\\\
-\frac{1}{2}&\frac{1}{2}&0&0&0&0&0&0&-0.158512&0.158512&-0.125587&0.400137\\\
0&0&0&0&0&0&-\frac{1}{\sqrt{2}}&0&0.317023&-0.317023&0.800274&-0.251173\\\
0&0&0&0&0&0&0&-\frac{1}{\sqrt{2}}&0.&0.&0.5491&0.5491\\\
0&0&0&\frac{1}{2}&\frac{1}{\sqrt{2}}&0&0&0&0.317023&-0.317023&-0.125587&0.400137\\\
0&0&\frac{1}{2}&0&0&\frac{1}{2\sqrt{2}}&0&0&0.217523&0.693058&0.&0.\\\
0&1&0&0&0&0&\frac{1}{\sqrt{2}}&0&0.&0.&0.&0.\\\
1&0&0&0&0&0&0&\frac{1}{\sqrt{2}}&0.&0.&0.&0.\\\ \end{array}\right)$ (146)

such that the diagonalized Hamiltonian has the form (A.2) where the matrix of
frequencies is

$\displaystyle\operatorname{\mathsf{\Omega}}=\begin{pmatrix}2.52434&0\\\
0&0.792287\\\ \end{pmatrix}.$ (147)

As previously advanced in the analysis of the main text, this system has two
free-particles (which correspond to the fact that we have two pairs of
capacitors in series and we have used a redudant description), two
nondynamical pairs of conjugated variables (associated with the intrinsically
redundant gyrator description), and two harmonic oscillators.

###  E.2 Keeping the full cosine

Let us now analyse the full problem with the Josephson junctions. One of the
routines suggested for the quantization of this problem is to perform
Williamson’s analysis only on the fully linear subsector, in this case,
$\operatorname{\mathsf{H}}_{\text{lin}}$ of Eq. (119), i.e., taking into
consideration only the inductor and coupling capacitor variables. This partial
classical diagonalization will be followed by a quantum one on the coupled
quantized system, thus the name two tier method. The partial symplectic
diagonalization of the Hamiltonian yields

$\displaystyle\operatorname{\mathsf{H}}_{pD}=\operatorname{\mathsf{S}}_{\operatorname{\boldsymbol{Y}}}^{T}\operatorname{\mathsf{H}}\operatorname{\mathsf{S}}_{\operatorname{\boldsymbol{Y}}}$
$\displaystyle=\begin{pmatrix}0&0&0&0&0\\\
0&\operatorname{\mathsf{\Omega}}&0&0&\operatorname{\mathsf{G}}\\\
0&0&\operatorname{\mathsf{\Omega}}&0&\operatorname{\mathsf{M}}\\\ 0&0&0&0&0\\\
0&\operatorname{\mathsf{G}}^{T}&\operatorname{\mathsf{M}}^{T}&0&\operatorname{\mathsf{\Omega}}_{J}\end{pmatrix},$
(148)

with particular numbers for the set of parameters chosen

$\displaystyle\operatorname{\mathsf{H}}_{pD}\stackrel{{\scriptstyle
num}}{{=}}\left(\begin{array}[]{cccccccccccc}0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&2.52434&0&0&0&0&0&0.800274&0.5491\\\
0&0&0&0&0&0.792287&0&0&0&0&-0.251173&0.5491\\\
0&0&0&0&0&0&2.52434&0&0&0&-0.317023&0\\\
0&0&0&0&0&0&0&0.792287&0&0&0.317023&0\\\ 0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0.800274&-0.251173&-0.317023&0.317023&0&0&1.&0\\\
0&0&0&0&0.5491&0.5491&0&0&0&0&0&1.\\\ \end{array}\right),$ (161)

where the coordinates have been transformed accordingly to the partially-
diagonalizing symplectic transformation such that
$(\operatorname{\mathsf{S}}_{\operatorname{\boldsymbol{Y}}}^{-1}\bar{\operatorname{\boldsymbol{X}}})^{T}=(\bar{\operatorname{\boldsymbol{y}}}^{T},\bar{\operatorname{\boldsymbol{\pi}}}^{T},\bar{\operatorname{\boldsymbol{\Phi}}}_{J}^{T},\bar{\operatorname{\boldsymbol{\Pi}}}_{J}^{T})$,
and
$\bar{\operatorname{\boldsymbol{y}}}^{T}=(\bar{\operatorname{\boldsymbol{y}}}_{\mathrm{nd}}^{T},\bar{\operatorname{\boldsymbol{y}}}_{\mathrm{ho}}^{T})$
are diagonalized position coordinates conjugated to the canonical momenta
$\bar{\operatorname{\boldsymbol{\pi}}}^{T}=(\bar{\operatorname{\boldsymbol{\pi}}}_{\mathrm{nd}}^{T},\bar{\operatorname{\boldsymbol{\pi}}}_{\mathrm{ho}}^{T})$.
After the canonical quantization of purely oscillator conjugated pairs and
their transformation to annihilation and creation operators, i.e.,
$\bar{y}_{i}=\sqrt{\frac{\hbar}{2}}(a_{i}+a_{i}^{\dagger})$, and
$\bar{\pi}_{i}=i\sqrt{\frac{\hbar}{2}}(a_{i}^{\dagger}-a_{i})$, we obtain

$\displaystyle\hat{H}/\hbar\stackrel{{\scriptstyle\text{q.}}}{{=}}$
$\displaystyle\sum_{n,m}^{2}\left[\Omega_{n}a_{n}^{\dagger}a_{n}+(g_{nm}a_{n}+g_{nm}^{*}a_{n}^{\dagger})\hat{\bar{\Pi}}_{Jm}\right]+\hat{H}_{J}/\hbar,$
(162)

where
$g_{nm}=\sqrt{\frac{1}{2\hbar}}\left[\operatorname{\mathsf{G}}_{nm}-i\operatorname{\mathsf{M}}_{nm}\right]$.
The purely quadratic part of the Josephson Hamiltonian can be read from the
lower right corner $4\times 4$ submatrix of $\operatorname{\mathsf{H}}_{pD}$
(or $\operatorname{\mathsf{H}}$), and thus the Josephson part of the total
Hamiltonian becomes

$\displaystyle\hat{H}_{J}/\hbar=$
$\displaystyle\sum_{n=1}^{2}\hat{\bar{\Pi}}_{Jn}^{2}-E_{Jn}\cos(\hat{\bar{\varphi}}_{Jn}),$
(163)

where the renormalized phase variable is
$\bar{\varphi}_{Jl}=2\pi\frac{\bar{\Phi}_{Jl}}{\Phi_{0}}$.

###  E.3 Black-box quantization approach

Let us now present a black-box quantization approach. Under the condition that
the Josephson flux variables have small fluctuations, we carry out an extended
harmonic analysis. To this point, extract the quadratic term of the cosine
potentials, and include its contribution,
$\bar{\operatorname{\boldsymbol{\Phi}}}_{J}^{T}\operatorname{\mathsf{L}}_{J}^{-1}\bar{\operatorname{\boldsymbol{\Phi}}}_{J}$,
to the Hamiltonian matrix (118), with

$\displaystyle\qquad\qquad\operatorname{\mathsf{L}}_{J}^{-1}=$
$\displaystyle\sqrt{\frac{L}{L_{J}}}\Omega_{J}\text{1}_{2}.$ (164)

Setting $L_{J}$, as way of illustration, to be equal to $L$, one would obtain