# Supplemental Material for
"Competition between two-photon driving, dissipation and interactions in
bosonic lattice models: an exact solution"

David Roberts1,2, A. A. Clerk  Pritzker School of Molecular Engineering,
University of Chicago, Chicago, IL, USA
2Department of Physics, University of Chicago, Chicago, IL, USA

(August 27, 2024)

###### Contents

  1. I Augmenting the model with hopping terms
    1. I.1 Imaginary hopping arrangements
    2. I.2 Chirally-symmetric hopping arrangements
  2. II Hidden TRS conditions
    1. II.1 Exact solution via representation theory
    2. II.2 Diagonal form of the representation
  3. III Exact solution: collective moments
    1. III.1 Unitary case
    2. III.2 Nonunitary case
  4. IV Exact solution: local moments
    1. IV.1 Addition of angular momentum for $SU(1,1)$
    2. IV.2 Unitary case
    3. IV.3 Nonunitary case
  5. V Phenomenology of the exact solution
    1. V.1 $SU(1,1)$ coherent states
    2. V.2 DMFT analysis of the $D=0$ model
    3. V.3 Semiclassical limit
  6. VI Mathematical background
    1. VI.1 Proof of the $SU(1,1)$ decomposition theorem
    2. VI.2 Exact solution for the steady-state Wigner function
    3. VI.3 Using the Wigner function to verify the nonthermal character of the steady state
  7. VII Experimental realization using superconducting circuits ($D=0$)
    1. VII.1 Adding a two-photon drive
    2. VII.2 Effect of junction capacitances

##  I Augmenting the model with hopping terms

Figure S1: Exactly-solvable Bose-Hubbard models with hopping and pair-driving.
(a) A bipartite lattice is subjected to chirally-symmetric hopping (depicted
using blue and red arrows) as well as pair-driving of the form
$\delta\hat{H}=\sum_{j}(\hat{a}_{j}^{\dagger}\hat{b}_{j}^{\dagger}+h.c.)$ (not
shown in this subpanel), and an infinite-range Bose-Hubbard interaction (not
shown in this subpanel). (b) At late times, the steady state becomes
stationary with respect to the hopping process depicted in panel (a), and so
the steady state can be analyzed by ignoring the hopping processes and
considering only the pair-driving process (depicted by blue ovals) and the
Hubbard interaction (not shown in this subpanel) in isolation.

The driven-dissipative bosonic lattice model that we solve in this work is
given by the following Lindblad master equation:

$\displaystyle\partial_{t}\hat{\rho}$
$\displaystyle=-i[\hat{H},\hat{\rho}]+\sum_{j}\mathcal{D}[\hat{L}_{j}],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\hat{H}=u\hat{N}^{2}-\Delta\hat{N}+\sum_{ij}\big{(}M_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}+h.c.\big{)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{L}_{j}=\sqrt{\kappa}\hat{a}_{j}$ (S1)

where here, $\hat{N}\equiv\sum_{j=1}^{N}\hat{a}_{j}^{\dagger}\hat{a}_{j}$
denotes total boson number, and $u\equiv U/N$. In addition, we have an
arbitrary complex-valued two-mode squeezing array $M_{ij}$.

We will now describe how single-photon tunneling terms of the form
$J\hat{a}_{i}^{\dagger}\hat{a}_{j}+h.c.$ can be added to the Hamiltonian while
preserving its solvability. A proper explanation of this necessitates a
discussion of the symmetries of the dissipative evolution generated by the
Lindbladian $\mathcal{L}$. We begin by discussing the two-mode squeezing array
$M$. We define the symmetry group $S$ of the array $M$ to be the group formed
by unitary beam-splitter transformations that leave the pairing array
invariant. Under such a transformation, represented by a unitary matrix $W$,
the pairing array $M$ transforms as $M\to WMW^{T}$. Therefore

$\displaystyle S$ $\displaystyle=\bigg{\\{}W\in U(N):\leavevmode\nobreak\
WMW^{T}=M\bigg{\\}}=V\prod_{\lambda\in\Sigma}O(s_{\lambda},\mathbb{R})\,V^{\dagger},$
(S2)

where $M=V\Sigma V^{T}$ is the Autonne-Takagi factorization of the matrix $M$
with $\Sigma$ the diagonal matrix of singular values, and where $s_{\lambda}$
denotes the degeneracy of the singular value $\lambda\in\Sigma$. In arriving
at the result (S2) we have used the fact that the symmetry group of $\Sigma$
is the product group $\prod_{\lambda\in\Sigma}O(s_{\lambda},\mathbb{R})$.
Since $S$ by definition conserves total particle number, any beam-splitter
transformation that lies in $S$ is a weak symmetry of the Lindbladian (S1).

Let us now assume that $\mathcal{L}$ has a unique steady state. Then this
steady state must respect the weak symmetries of the Lindbladian $\mathcal{L}$
Buča and Prosen (2012). In particular, the steady state is invariant under the
symmetry group $S$ of the pairing array. In particular, we can add any element
in the Lie algebra $\mathfrak{s}$ of $S$ to the Hamiltonian in (S1) while
still preserving the steady state. The Lie algebra of $S$ can be explicitly
computed:

$\displaystyle\mathfrak{s}=\mathcal{V}\bigg{\\{}i\sum_{\lambda\in\Sigma}\sum_{i,j\in\underline{s}_{\lambda}}t_{ij}(\hat{a}_{i}^{\dagger}\hat{a}_{j}-\hat{a}_{j}^{\dagger}\hat{a}_{i}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
t_{ij}\in\mathbb{R}\bigg{\\}}\mathcal{V}^{\dagger}.$ (S3)

where $\underline{s}_{\lambda}$ denotes the subset of mode indices
corresponding to the singular value $s_{\lambda}$, and $\mathcal{V}$ is some
beam-splitter transformation that implements the unitary transformation $V$,
i.e.
$\mathcal{V}^{\dagger}\hat{a}_{i}\mathcal{V}\equiv\sum_{j}V_{ij}\hat{a}_{j}$.

###  I.1 Imaginary hopping arrangements

When the singular values of the pairing array are fully degenerate, the
pairing array is maximally symmetric: in this case, the symmetry group $S$ is,
up to the unitary similarity transform $V$, just the full rotation group as a
subgroup of the unitary group. The simplest example is the case where the
pairing array is diagonal and uniform: $M_{ij}=G\delta_{ij}$, which lets us
set $V=1$. Therefore, in this case, the symmetry group of the pairing array is
generated by terms of the form

$\displaystyle\delta\hat{H}\equiv
i\sum_{i,j=1}^{N}t_{ij}(\hat{a}_{i}^{\dagger}\hat{a}_{j}-\hat{a}_{j}^{\dagger}\hat{a}_{i}).$
(S4)

Recall that, since $\delta\hat{H}$ generates a weak symmetry of $\mathcal{L}$,
one can add $\delta\hat{H}$ to the Hamiltonian without changing the steady
state of $\mathcal{L}$. That is, one can add any set of imaginary hopping
terms to the Hamiltonian without changing the steady state of $\mathcal{L}$.

###  I.2 Chirally-symmetric hopping arrangements

Thanks to the similarity transform $V$ appearing in the symmetry group $S$,
there are many situations where our solvable model can be augmented with
conventional hopping terms. As a particularly striking example of this,
consider a global Bose-Hubbard model defined on a bipartite lattice with $2N$
sites, with corresponding modes
$\hat{a}_{1},\dots\hat{a}_{N},\hat{b}_{1},\dots,\hat{b}_{N}$. We can consider
a special case of the Lindbladian (S1) where the pairing array is fully
"dimerized", that is, where the Hamiltonian takes the form

$\displaystyle\hat{H}=u\hat{N}^{2}-\Delta\hat{N}+\Lambda\sum_{j}(\hat{a}_{j}^{\dagger}\hat{b}_{j}^{\dagger}+h.c.).$
(S5)

In this case, one can augment the above Hamiltonian with the following hopping
terms:

$\displaystyle\delta\hat{H}\equiv\sum_{ij}J_{ij}\bigg{(}(\hat{a}_{i}^{\dagger}\hat{a}_{j}+h.c.)-(\hat{b}_{i}^{\dagger}\hat{b}_{j}+h.c.)\bigg{)}.$
(S6)

Note that such a model describes arbitrary hopping within the $A$-sublattice,
which is mirrored in the $B$-sublattice. This situation is schematically
depicted in Fig. S1a. To see explicitly that the hopping terms generate
rotational symmetries of $\mathcal{L}$, we have to diagonalize the pairing
array. In particular, one can do this by defining modes

$\displaystyle\hat{c}_{j,+}^{\dagger}\equiv\frac{\hat{a}_{j}^{\dagger}+\hat{b}_{j}^{\dagger}}{\sqrt{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{c}_{j,-}^{\dagger}\equiv\frac{\hat{a}_{j}^{\dagger}-\hat{b}_{j}^{\dagger}}{i\sqrt{2}}.$
(S7)

In this mode basis, the Hamiltonian, with the hopping terms included, reads as

$\displaystyle\hat{H}=u\hat{N}^{2}-\Delta\hat{N}+\Lambda\sum_{j,\pm}(\hat{c}_{j,\pm}^{\dagger
2}+h.c.)+i\sum_{i,j=1}^{N}J_{ij}(c_{i,-}^{\dagger}\hat{c}_{j,+}-h.c.)+i\sum_{i,j=1}^{N}J_{ij}(c_{i,+}^{\dagger}\hat{c}_{j,-}-h.c.).$
(S8)

In this diagonal mode basis, it is extremely clear that the hopping terms
generate rotations of the modes $\hat{c}_{j,\pm},\hat{c}_{j,\pm}^{\dagger}$
that are weak symmetries of the Lindbladian. Therefore, the results of this
work can be used to analytically describe the steady states of the model (S5),
including the chirally-symmetric hopping arrangement depicted in Fig. S1a.

##  II Hidden TRS conditions

Given a steady-state $\hat{\rho}_{\text{ss}}$ and corresponding time-reversal
operation $\hat{T}$, one can construct the thermofield double state

$\displaystyle|\Psi_{\hat{T}}\rangle=\sum_{n}\sqrt{p_{n}}|n\rangle_{L}\hat{T}|n\rangle_{R},$
(S9)

where $|n\rangle$ are the eigenvectors of $\hat{\rho}_{\text{ss}}$, and
$p_{n}$ the corresponding eigenvalues. Sufficient conditions for $\hat{T}$ to
be a hidden time-reversal symmetry are Fagnola and Umanità (2010)

$\displaystyle\Bigg{[}\bigg{(}u(\hat{N}_{L}+\hat{N}_{R})-\Delta_{\text{eff}}\bigg{)}(\hat{N}_{L}-\hat{N}_{R})+2\sum_{ij}M_{ij}\hat{\alpha}_{i,+}^{\dagger}\hat{\alpha}_{j,-}^{\dagger}\Bigg{]}|\Psi_{\hat{T}}\rangle$
$\displaystyle=0,$ $\displaystyle\hat{\alpha}_{j,-}|\Psi_{\hat{T}}\rangle$
$\displaystyle=0,$ (S10)

where $\hat{\alpha}_{j,\pm}\equiv 2^{-1/2}(\hat{a}_{j,L}\pm\hat{a}_{j,R})$.
Since the Bogoliubov transformation
$\hat{a}_{j,L},\hat{a}_{j,R}\to\hat{\alpha}_{j,+},\hat{\alpha}_{j,-}$ is
boson-number preserving, we have
$\hat{N}_{L}+\hat{N}_{R}=\hat{N}_{+}+\hat{N}_{-}$, where
$\hat{N}_{\pm}:=\sum_{j}\hat{\alpha}_{j,\pm}^{\dagger}\hat{\alpha}_{j,\pm}$.
Therefore, we can write the above conditions equivalently as

$\displaystyle\sum_{j}\hat{\alpha}_{j,-}^{\dagger}\Bigg{[}\bigg{(}u(\hat{N}_{+}+\hat{N}_{-}+1)-\Delta_{\text{eff}}\bigg{)}\hat{\alpha}_{j,+}+2\sum_{i}M_{ij}\hat{\alpha}_{i,+}^{\dagger}\Bigg{]}|\Psi_{\hat{T}}\rangle$
$\displaystyle=0,$ (S11)
$\displaystyle\hat{\alpha}_{j,-}|\Psi_{\hat{T}}\rangle$ $\displaystyle=0.$
(S12)

Each term in the sum in Eq. (S11) has exactly one boson in the antisymmetric
mode $\hat{\alpha}_{j,-}$ and is thus mutually orthogonal with the remaining
terms. Therefore, each term in the sum must separately vanish. Factoring the
thermofield double into symmetric and antisymmetric components
$|\Psi_{\hat{T}}\rangle=|\Psi_{+}\rangle|\Psi_{-}\rangle$, we obtain
$|\Psi_{-}\rangle=|0_{-}\rangle$, i.e. the antisymmetric component is in
vacuum, and

$\displaystyle\bigg{(}(u(\hat{N}_{+}+1)-\Delta_{\text{eff}})\hat{\alpha}_{j,+}+2\sum_{i}M_{ij}\hat{\alpha}_{i,+}^{\dagger}\bigg{)}|\Psi_{+}\rangle=0.$
(S13)

###  II.1 Exact solution via representation theory

We now proceed to solve the system of equations Eq. (S13). To begin, we
multiply each constraint on the left by $\hat{\alpha}_{j,+}^{\dagger}$,
yielding

$\displaystyle\bigg{(}u(\hat{N}_{+}-\Delta_{\text{eff}})\hat{n}_{j,+}+2\sum_{i}M_{ij}\hat{\alpha}_{i,+}^{\dagger}\hat{\alpha}_{j,+}^{\dagger}\bigg{)}|\Psi_{+}\rangle=0,$
(S14)

where $\hat{n}_{j,+}:=\hat{\alpha}_{j,+}^{\dagger}\hat{\alpha}_{j,+}$. We then
sum over $j$:

$\displaystyle\bigg{(}u(\hat{N}_{+}-\Delta_{\text{eff}})\hat{N}_{+}+2\sum_{ij}M_{ij}\hat{\alpha}_{i,+}^{\dagger}\hat{\alpha}_{j,+}^{\dagger}\bigg{)}|\Psi_{+}\rangle=0.$
(S15)

Note how the above operator annihilating $|\Psi_{\hat{T}}\rangle$ is almost
like the effective nonhermitian Hamiltonian
$\hat{H}_{\text{eff}}:=\hat{H}-i\kappa\hat{N}/2$, except that here the
lowering operators associated with the drive are not present. We now identify
a "hidden" nonunitary representation of $SU(1,1)$ via

$\displaystyle\hat{K}_{+}$
$\displaystyle=\frac{1}{2}\sum_{ij}\frac{M_{ij}}{u}\hat{\alpha}_{i,+}^{\dagger}\hat{\alpha}_{j,+}^{\dagger},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{-}=\frac{1}{2}\sum_{ij}(uM^{-1})_{ij}\hat{\alpha}_{i,+}\hat{\alpha}_{j,+},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{z}=\frac{1}{2}\sum_{j}\bigg{(}\hat{\alpha}_{j,+}^{\dagger}\hat{\alpha}_{j,+}+\frac{1}{2}\bigg{)}.$
(S16)

Note that the representation presented here is unitary if and only if the
pairing array is unitary, i.e. $(M/u)^{-1}=(M/u)^{\dagger}$. In terms of the
Lie algebra generators, the hTRS condition simplifies to
$\mathcal{H}_{\text{eff},+}|\Psi_{+}\rangle=0$, where

$\displaystyle\mathcal{H}_{\text{eff},+}:=\big{(}\hat{K}_{z}-N/4-\Delta_{\text{eff}}/2u\big{)}\big{(}\hat{K}_{z}-N/4\big{)}+\hat{K}_{+}.$
(S17)

To solve this condition, we make the ansatz that the dark state inhabits an
irreducible representation of $SU(1,1)$, with the vacuum state constituting
the vector of lowest weight (with respect to the subalgebra generated by
$\hat{K}_{z}$). This irreducible representation is spanned by vectors of the
form

$\displaystyle|\Psi\rangle=\sum_{m}c_{m}\hat{K}_{+}^{m}|h_{0}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
|h_{0}\rangle:=|\Omega\rangle$ (S18)

The use of the notation $h_{0}$ to denote the bosonic vacuum
$|\Omega\rangle:=|0\rangle_{L}|0\rangle_{R}$, as well as the fact that the
solution Eq. (S18) to the constraints Eq. (S13) is unique will become apparent
in later sections when we review the general representation theory of
$SU(1,1)$. For now, however, we will treat Eq. (S18) as an ansatz and proceed
to compute the exact coefficients $c_{m}$ occuring in the expansion:

$\displaystyle\mathcal{H}_{\text{eff},+}|\Psi_{+}\rangle=\sum_{m=1}^{\infty}\Bigg{(}m(m-\Delta_{\text{eff}}/2u)c_{m}+c_{m-1}\Bigg{)}\hat{K}_{+}^{m}|h_{0}\rangle=0,$
(S19)

which yields the solution $c_{m}\propto(-1)^{m}/m!(\delta)_{m}$, where
$\delta:=1-\Delta_{\text{eff}}/2u$, and
$(\delta)_{m}:=\delta(\delta+1)\cdots(\delta+m-1)$ denotes the Pochhammer
symbol (rising factorial).

###  II.2 Diagonal form of the representation

By treating the raising operators $\hat{\alpha}_{j,+}^{\dagger}$ as
indeterminates, one can interpret the $SU(1,1)$ raising operator $\hat{K}_{+}$
presented in Eq. (S16) as a quadratic form over the complex numbers. Its
matrix representation $M/u$ is invariant under transposition and thus admits a
kind of singular-value decomposition $M/u=V\Sigma V^{T}$ (called the Autonne-
Takagi factorization Autonne (1915)), with $V$ an $N\times N$ unitary matrix
and $\Sigma=\text{diag}(\lambda_{1},\dots,\lambda_{N})$ the matrix of singular
values. Writing
$\hat{\beta}_{i,+}^{\dagger}:=\sum_{i}V_{ij}\hat{\alpha}_{j,+}^{\dagger}$ for
the corresponding basis of singular vectors, we can re-express the $SU(1,1)$
representation (S16) as follows:

$\displaystyle\hat{K}_{+}$
$\displaystyle=\frac{1}{2}\sum_{j}\lambda_{j}\hat{\beta}_{j,+}^{\dagger
2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\hat{K}_{-}=\frac{1}{2}\sum_{i}\lambda_{j}^{-1}\hat{\beta}_{j,+}^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{z}=\frac{1}{2}\sum_{j}\bigg{(}\hat{\beta}_{j,+}^{\dagger}\hat{\beta}_{j,+}+\frac{1}{2}\bigg{)}.$
(S20)

##  III Exact solution: collective moments

The simplest quantities to compute in our model are collective observables,
that is, moments of normally-ordered monomials in the $SU(1,1)$ algebra.

###  III.1 Unitary case

Our job is easiest when the representation Eq. (S16) satisfies
$\hat{K}_{+}^{\dagger}\propto\hat{K}_{-}$. This happens if and only if the
singular values of the matrix $M/u$ are all the same. In this case, the
$SU(1,1)$ representation can be made unitary by rescaling $\hat{K}_{-}$, and
the diagonal form of the representation is as follows:

$\displaystyle\hat{K}_{+}=\frac{1}{2}\sum_{j}\lambda\hat{\beta}_{j,+}^{\dagger
2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\hat{K}_{-}=\frac{1}{2}\sum_{j}\lambda^{-1}\hat{\beta}_{j,+}^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{z}=\frac{1}{2}\sum_{j}\bigg{(}\hat{\beta}_{j,+}^{\dagger}\hat{\beta}_{j,+}+\frac{1}{2}\bigg{)},$
(S21)

with $\lambda>0$ the unique singular value of the matrix $M/u$. With this, we
now turn to compute moments of monomials of $SU(1,1)$ within the
representation:

$\displaystyle\langle\Psi_{\hat{T}}|\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}|\Psi_{\hat{T}}\rangle=\sum_{p,q}\lambda^{2p}c_{p}^{*}c_{q}\langle
h_{0}|\hat{K}_{-}^{p}\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}\hat{K}_{+}^{q}|h_{0}\rangle$
(S22)

We reparametrize the sum by defining an auxiliary non-negative integer
$l=p-n=q-m$. In terms of $l$:

$\displaystyle\langle\Psi_{\hat{T}}|\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}|\Psi_{\hat{T}}\rangle$
$\displaystyle=\sum_{l=0}^{\infty}\lambda^{2(l+n)}c_{l+n}^{*}c_{l+m}\langle
h_{0}|\hat{K}_{-}^{l}\hat{K}_{-}^{n}\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}\hat{K}_{+}^{m}\hat{K}_{+}^{l}|h_{0}\rangle$
$\displaystyle=\lambda^{2n}c_{n}^{*}c_{m}\sum_{l=0}^{\infty}(2l)^{k}|c_{l}|^{2}\lambda^{2l}\frac{c_{n+l}^{*}}{c_{l}^{*}c_{n}^{*}}\frac{c_{m+l}}{c_{l}c_{m}}f(n,m,l).$
(S23)

Using the $SU(1,1)$ commutation relations, we evaluate

$\displaystyle f(n,m,l)$ $\displaystyle:=\langle
h_{0}|\hat{K}_{-}^{l}\hat{K}_{-}^{n}\hat{K}_{+}^{n}\hat{K}_{+}^{m}\hat{K}_{-}^{m}\hat{K}_{+}^{l}|h_{0}\rangle$
$\displaystyle=l!(N/2)_{l}\frac{(n+l)!}{l!}\frac{(m+l)!}{l!}\frac{(N/2)_{n+l}}{(N/2)_{l}}\frac{(N/2)_{m+l}}{(N/2)_{l}}=\frac{(N/2)_{n}(N/2)_{m}}{l!}\frac{(n+l)!(m+l)!(N/2+n)_{l}(N/2+m)_{l}}{(N/2)_{l}}.$
(S24)

We then compute

$\displaystyle(2l)^{k}|c_{l}|^{2}\lambda^{2l}\frac{c_{n+l}^{*}}{c_{l}^{*}c_{n}^{*}}\frac{c_{m+l}}{c_{l}c_{m}}=|c_{l}|^{2}\lambda^{2l}\frac{l!n!l!m!}{(n+l)!(m+l)!}\frac{(\delta^{*})_{l}(\delta^{*})_{n}(\delta)_{l}(\delta)_{m}}{(\delta^{*})_{l+n}(\delta^{*})_{l+m}}=\frac{\lambda^{2l}n!m!}{(n+l)!(m+l)!(\delta^{*}+n)_{l}(\delta+m)_{l}}.$
(S25)

Putting everything together,

$\displaystyle\langle\Psi_{\hat{T}}|\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}|\Psi_{\hat{T}}\rangle$
$\displaystyle=\lambda^{2n}c_{n}^{*}c_{m}(N/2)_{n}(N/2)_{m}n!m!\sum_{l=0}^{\infty}(2l)^{k}\frac{\lambda^{2l}}{l!}\frac{(N/2+n)_{l}(N/2+m)_{l}}{(N/2)_{l}(\delta^{*}+n)_{l}(\delta+m)_{l}}$
$\displaystyle=\frac{\lambda^{2n}(-1)^{n+m}(N/2)_{n}(N/2)_{m}}{(\delta^{*})_{n}(\delta)_{m}}\bigg{(}2z\frac{d}{dz}\bigg{)}^{k}\,_{2}F_{3}(N/2+n,N/2+m;N/2,\delta^{*}+n,\delta+m;z)\Bigg{|}_{z=\lambda^{2}},$
(S26)

where here, $\,{}_{p}F_{q}$ denotes the generalized hypergeometric function,
defined as the analytic extension of the following infinite series:

$\,{}_{p}F_{q}(a_{1},\dots a_{p};b_{1},\dots
b_{q};z)=\sum_{l=0}^{\infty}\frac{\prod_{m=1}^{p}(a_{m})_{l}}{\prod_{m=1}^{q}(b_{m})_{l}}\frac{z^{l}}{l!}.$
(S27)

Up to this point, we have computed $|\Psi_{\hat{T}}\rangle$ up to a
normalization prefactor. This normalization can also be computed in closed
form:
$\langle\Psi_{\hat{T}}|\Psi_{\hat{T}}\rangle=\,_{1}F_{2}(N/2;\delta,\delta^{*};\lambda^{2})$,
thus completing the characterization of collective moments in the case that
the representation (S16) is unitary.

###  III.2 Nonunitary case

We now accomplish the same task, but in the more general regime where the
singular values $\lambda_{j}$ are possibly non-degenerate. In this situation,
there is less structure to the calculation, but moments of $SU(1,1)$ monomials
may nonetheless be extracted efficiently:

$\displaystyle\langle\Psi_{\hat{T}}|\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}|\Psi_{\hat{T}}\rangle$
$\displaystyle=\sum_{p,q}^{\infty}c_{p}^{*}c_{q}\langle
h_{0}|\hat{K}_{+}^{\dagger
p}\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}\hat{K}_{+}^{q}|h_{0}\rangle$
(S28)

We reparametrize the sum by defining an auxiliary non-negative integer
$l=p-n=q-m$. In terms of $l$:

$\displaystyle\langle\Psi_{\hat{T}}|\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}|\Psi_{\hat{T}}\rangle$
$\displaystyle=\sum_{l=0}^{\infty}c_{l+n}^{*}c_{l+m}(2l)^{k}\langle
h_{0}|\hat{K}_{+}^{\dagger l}\hat{K}_{+}^{\dagger
n}\hat{K}_{+}^{n}\hat{K}_{-}^{m}\hat{K}_{+}^{m}\hat{K}_{+}^{l}|h_{0}\rangle$
$\displaystyle=c_{n}^{*}c_{m}\sum_{l=0}^{\infty}|c_{l}|^{2}\frac{c_{l+n}^{*}}{c_{l}^{*}c_{n}^{*}}\frac{c_{l+m}}{c_{l}c_{m}}(2l)^{k}g(n,m,l).$
(S29)

We evaluate:

$\displaystyle g(n,m,l)$
$\displaystyle=\frac{(N/2)_{m+l}(m+l)!}{(N/2)_{l}l!}\langle
h_{0}|\hat{K}_{+}^{\dagger(n+l)}\hat{K}_{+}^{n+l}|h_{0}\rangle:=\frac{(N/2)_{m+l}(m+l)!(n+l)!(n+l)!}{(N/2)_{l}l!}\Phi_{n+l}(\vec{\lambda};N),$
(S30)

where $\Phi_{n}(\vec{\lambda};N):=(n!)^{2}\langle h_{0}|\hat{K}_{+}^{\dagger
n}\hat{K}_{+}^{n}|h_{0}\rangle$ is a form factor to be evaluated later on.
Therefore, in total we have

$\displaystyle\langle\Psi_{\hat{T}}|\hat{K}_{+}^{n}\hat{N}_{+}^{k}\hat{K}_{-}^{m}|\Psi_{\hat{T}}\rangle=\frac{(-1)^{n+m}(N/2)_{m}}{(\delta^{*})_{n}(\delta)_{m}}\sum_{l=0}^{\infty}\frac{(N/2+m)_{l}(n+l)!(2l)^{k}}{(\delta^{*}+n)_{l}(\delta+m)_{l}(N/2)_{l}l!}\Phi_{n+l}(\vec{\lambda};N).$
(S31)

All that is left is to compute the normalization, which comes out to
$\langle\Psi_{\hat{T}}|\Psi_{\hat{T}}\rangle=\sum_{l=0}^{\infty}\Phi_{l}(\vec{\lambda};N)/(\delta)_{l}(\delta^{*})_{l}$.
Up until this point, we have given formulae in terms of the form factors
$\Phi_{n}(\vec{\lambda};N):=\langle h_{0}|\hat{K}_{+}^{\dagger
n}\hat{K}_{+}^{n}|h_{0}\rangle/(n!)^{2}$, which, in the nonunitary case,
cannot be evaluated purely using the internal relations of the $SU(1,1)$
algebra. Instead, via the multinomial expansion, we are left with a sum over
an exponential number of terms:

$\displaystyle\Phi_{n}(\vec{\lambda};N)=\frac{1}{n!^{2}}\langle
h_{0}|\Bigg{(}\frac{1}{2}\sum_{j}\lambda_{j}\hat{\beta}_{j,+}^{2}\Bigg{)}^{n}\Bigg{(}\frac{1}{2}\sum_{j}\lambda_{j}\hat{\beta}_{j,+}^{\dagger
2}\Bigg{)}^{n}|h_{0}\rangle=\sum_{\sum_{j}k_{j}=n}\prod_{j}(1/2)_{k_{j}}\lambda_{j}^{2k_{j}}$
(S32)

Our work would not be complete if we did not give an efficient prescription
for evaluating these sums. The key to evaluating these sums efficiently is the
following observation:

$\displaystyle\Phi_{n}(\vec{\lambda};N)=\sum_{p=0}^{n}(1/2)_{p}\lambda_{N}^{2p}\sum_{\sum_{j=1}^{N-1}k_{j}=n-p}\prod_{j}(1/2)_{k_{j}}\lambda_{j}^{2k_{j}}=\sum_{p=0}^{n}(1/2)_{p}\lambda_{N}^{2p}\Phi_{n-p}(\vec{\lambda};N-1)$
(S33)

We note that the above constitutes a recursion relation for
$\Phi_{n}(\vec{\lambda};N)$ with boundary condition
$\Phi_{p}(\vec{\lambda};1)=(1/2)_{p}\lambda_{1}^{2p}$. Note that via this
recursion relation, and fixing some boson-number cutoff $k$, the collection
$\\{\Phi_{p}(\vec{\lambda};N)\\}_{p\leq k}$ of form factors may be evaluated
using only $O(k^{2}N)$ floating-point operations.

##  IV Exact solution: local moments

We now turn to a more difficult task: that of computing (equal-time)
correlation functions of local observables, that is, observables that are not
collective in nature. To accomplish this efficiently in the unitary case, we
must solve the "addition of angular momentum" problem for $SU(1,1)$, i.e. we
must understand how to decompose a tensor product of $SU(1,1)$ representations
into irreducible components. Luckily, this task is easily solved in terms of
the theory of harmonic functions on $\mathbb{R}^{N}$. In the nonunitary case,
however, the $SU(1,1)$ structure is not so useful, and instead some
generalization of the combinatorics in (S33) is needed to compute observables.

To make our task easier, we will compute observables in the basis of modes
$\hat{b}_{i}^{\dagger}:=\sum_{j}V_{ij}\hat{a}_{j}^{\dagger}$ that diagonalizes
the pair-driving matrix $M/u$.

###  IV.1 Addition of angular momentum for $SU(1,1)$

We can write our global $SU(1,1)$ representation (S20) as a tensor product of
local representations:

$\displaystyle\hat{K}_{+}^{(j)}\equiv\frac{1}{2}\lambda_{j}\hat{\beta}_{j,+}^{\dagger
2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\hat{K}_{-}^{(j)}\equiv\frac{1}{2}\lambda^{-1}_{j}\hat{\beta}_{j,+}^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{z}^{(j)}\equiv\frac{1}{2}\bigg{(}\hat{\beta}_{j,+}^{\dagger}\hat{\beta}_{j,+}+\frac{1}{2}\bigg{)}$
(S34)

It is easy to see that each local representation is reducible, and has the
decomposition $V^{(j)}=V^{(j)}_{+}\oplus V^{(j)}_{-}$, where $V^{(j)}_{\pm}$
denotes the subspace consisting of all states with a fixed boson-number parity
$\pm 1$. Our global $SU(1,1)$ representation
$\hat{K}_{+},\hat{K}_{-},\hat{K}_{z}$ is the $N$-fold tensor product
$\otimes_{j}V^{(j)}$, and has the following decomposition into irreducible
subrepresentations:

$\displaystyle\otimes_{j}V^{(j)}\simeq\bigoplus_{l=0}^{\infty}\Bigg{(}\bigoplus_{p=1}^{d_{l}}V_{l}^{(p)}\Bigg{)},$
(S35)

where each irreducible subrepresentation $V_{l}^{(p)}$ takes the form

$\displaystyle|\Psi_{l}^{(p)}\rangle=\sum_{m=0}^{\infty}c_{m}\hat{K}_{+}^{m}|h_{l}^{(p)}\rangle,$
(S36)

and $|h_{l}^{(1)}\rangle,\dots,|h_{l}^{(d_{l})}\rangle$ is some orthonormal
basis of the subspace of the kernel of $\hat{K}_{-}$ consisting of those
states with fixed total photon number equal to $l$. Later on we sketch a proof
of the decomposition (S35) using the Segal-Bargmann representation Bargmann
(1961, 1967), which represents bosonic states as multivariate analytic
functions, and bosonic creation- and annihilation operators as partial
differential operators. In the Segal-Bargmann representation, the $SU(1,1)$
vacua $|h_{l}^{(p)}\rangle$ are represented as homogeneous harmonic
polynomials of degree $l$, hence the notation "$h_{l}$".

###  IV.2 Unitary case

In the unitary case, the decomposition (S35) becomes an orthogonal direct-sum
decomposition, with the superselection rules $\langle
h_{l}^{(m)}|\hat{K}_{+}^{\dagger
p}\hat{K}_{+}^{q}|h_{l^{\prime}}^{(m^{\prime})}\rangle=\delta_{l,l^{\prime}}\delta_{m,m^{\prime}}\delta_{p,q}p!\lambda^{2p}(N/2+l)_{p}$.
These simple rules allow us to compute any local quantity of interest in our
model. For the case of quadratic correlation functions, however, such rules
are not necessary. To see this, it suffices to use the weak permutation
symmetry $\hat{b}_{i}\leftrightarrow\hat{b}_{j}$, as well as the weak parity
symmetry $\hat{b}_{i}\to-\hat{b}_{i}$ of the Lindbladian $\mathcal{L}$:

$\displaystyle\langle\hat{b}_{i}^{\dagger}\hat{b}_{i+r}\rangle$
$\displaystyle=\delta_{r,0}\frac{1}{N}\sum_{j}\langle\hat{b}_{j}^{\dagger}\hat{b}_{j}\rangle=\frac{\delta_{r,0}}{N}\langle\Psi_{\hat{T}}|\hat{N}_{+}|\Psi_{\hat{T}}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\langle\hat{b}_{i}\hat{b}_{i+r}\rangle=\delta_{r,0}\frac{1}{N}\sum_{j}\langle\hat{b}_{j}^{2}\rangle=\frac{\delta_{r,0}}{N}\langle\Psi_{\hat{T}}|\hat{K}_{-}|\Psi_{\hat{T}}\rangle$
(S37)

Higher-order correlations, are not expressible in terms of collective moments,
and so in general the decomposition (S35), along with the associated
superselection rules, serve as a useful guide. As an example, we compute the
pair-correlation function, as well as the density-density correlation
function, and leave the general case to the reader. We compute the pair
correlation function first:

$\displaystyle\langle\hat{b}_{i}^{\dagger 2}\hat{b}_{i+r}^{2}\rangle$
$\displaystyle=\frac{1}{4}\langle\Psi_{\hat{T}}|\hat{\beta}_{i+r,+}^{2}\hat{\beta}_{i,+}^{\dagger
2}|\Psi_{\hat{T}}\rangle-\frac{2\delta_{r,0}}{N}\langle\Psi_{\hat{T}}|\hat{K}_{z}|\Psi_{\hat{T}}\rangle=\frac{1}{4}\sum_{p,q}c_{p}^{*}c_{q}\langle
h_{0}|\hat{\beta}_{i,+}^{2}\hat{K}_{+}^{\dagger
p}\hat{K}_{+}^{q}\hat{\beta}_{i+r,+}^{\dagger
2}|h_{0}\rangle-\frac{2\delta_{r,0}}{N}\langle\Psi_{\hat{T}}|\hat{K}_{z}|\Psi_{\hat{T}}\rangle$
(S38)

For the calculation to proceed, it is necessary to decompose the states
$\hat{\beta}_{i,+}^{\dagger 2}|h_{0}\rangle$ into $SU(1,1)$ vacua. Although
this can be done by hand in this simple case, in general it is useful to
automate this process (especially for higher-order observables). For this
purpose, the HFT Mathematica package Axler (2020) is especially relevant, in
particular the function $\text{harmonicDecomposition}[]$, which automatically
extracts the decomposition

$\displaystyle\hat{\beta}_{i,+}^{\dagger
2}|h_{0}\rangle=\sqrt{\frac{2N-2}{N}}|h_{2}^{(i)}\rangle+\frac{2\hat{K}_{+}}{\lambda
N}|h_{0}\rangle,$ (S39)

where
$|h_{2}^{(i)}\rangle:=\frac{1}{\lambda}\sqrt{\frac{2N}{N-1}}(\hat{K}_{+}^{(i)}-\frac{\hat{K}_{+}}{N})|h_{0}\rangle$
generates an irreducible subrepresentation with weight $N/4+1$. We thus
obtain:

$\displaystyle\frac{1}{4}\sum_{p,q}c_{p}^{*}c_{q}\langle
h_{0}|\hat{\beta}_{i,+}^{2}\hat{K}_{+}^{\dagger
p}\hat{K}_{+}^{q}\hat{\beta}_{i+r,+}^{\dagger 2}|h_{0}\rangle$
$\displaystyle=\frac{2N-2}{N^{3}\lambda^{2}}\sum_{l}|c_{l}|^{2}\langle
h_{0}|\hat{K}_{+}^{\dagger
l+1}\hat{K}_{+}^{l+1}|h_{0}\rangle+\frac{N-1}{2N}\sum_{l}|c_{l}|^{2}\langle
h_{2}^{(i)}|\hat{K}_{+}^{\dagger l}\hat{K}_{+}^{l}|h_{2}^{(i+r)}\rangle$
$\displaystyle=\Bigg{(}\frac{2N-2}{N^{3}}\bigg{[}z\frac{\partial}{\partial
z}+1\bigg{]}\,_{1}F_{2}(N/2;\delta,\delta^{*};z)+\frac{N\delta_{r,0}-1}{2N}\,_{1}F_{2}(N/2+2;\delta,\delta^{*};z)\Bigg{)}\Bigg{|}_{z=\lambda^{2}}.$
(S40)

We now compute the density-density correlation function
$\langle\hat{b}_{i}^{\dagger}\hat{b}_{i}\hat{b}_{i+r}^{\dagger}\hat{b}_{i+r}\rangle$.
Since the case $r=0$ is subsumed by the previous calculation, without loss of
generality we can assume $r\neq 0$. Again, we split the calculation into a
collective part and a noncollective part (that populates a higher weight
representation):

$\displaystyle\langle\hat{b}_{i}^{\dagger}\hat{b}_{i}\hat{b}_{i+r}^{\dagger}\hat{b}_{i+r}\rangle=\frac{1}{4}\sum_{p,q}c_{p}^{*}c_{q}\langle
h_{0}|\hat{\beta}_{i}\hat{\beta}_{i+r}\hat{K}_{+}^{\dagger
p}\hat{K}_{+}^{q}\hat{\beta}_{i}^{\dagger}\hat{\beta}_{i+r}^{\dagger}|h_{0}\rangle-\frac{1}{N}\langle\Psi_{\hat{T}}|\hat{K}_{z}|\Psi_{\hat{T}}\rangle.$
(S41)

We evaluate the noncollective part, by noticing that
$|h_{2}^{(i,j)}\rangle:=\hat{\beta}_{i,+}^{\dagger}\hat{\beta}_{j,+}^{\dagger}|h_{0}\rangle$
generates an irreducible subrepresentation with weight $N/4+1$:

$\displaystyle\frac{1}{4}\sum_{p,q}c_{p}^{*}c_{q}\langle
h_{0}|\hat{\beta}_{i}\hat{\beta}_{i+r}\hat{K}_{+}^{\dagger
p}\hat{K}_{+}^{q}\hat{\beta}_{i}^{\dagger}\hat{\beta}_{i+r}^{\dagger}|h_{0}\rangle$
$\displaystyle=\frac{1}{4}\sum_{l}|c_{l}|^{2}\langle
h_{2}^{(i,i+r)}|\hat{K}_{+}^{\dagger
l}\hat{K}_{+}^{l}|h_{2}^{(i,i+r)}\rangle=\frac{1}{4}\,_{1}F_{2}(N/2+2;\delta,\delta^{*};\lambda^{2}).$
(S42)

#### Higher-order moments

The general procedure for calculating higher order moments is no different
from what we have done in the previous subsection: given a higher-order moment
to be evaluated, one writes the expectation value in terms of an antinormally-
ordered correlation function involving the mode operators
$\hat{\beta}_{j,+},\hat{\beta}_{j,+}^{\dagger}$, and then uses standard
harmonic analysis software Axler (2020) to decompose the resulting states into
components lying in higher-weight subrepresentations. Iterating this process
yields formulae in terms of generalized hypergeometric functions and their
derivatives.

###  IV.3 Nonunitary case

We now turn to the most general task of computing local correlation functions
in our global Bose-Hubbard model, in the nonunitary regime. In this case the
calculation has the least amount of structure, and so here we just present the
most general result. To better organize the calculation, we will write the
purification $|\Psi_{\hat{T}}\rangle$ as a power series

$\displaystyle|\Psi_{\hat{T}}\rangle=\sum_{\vec{m}\in\mathbb{N}^{N}}\frac{c_{\vec{m}}}{\vec{m}!}\prod_{j}\hat{\beta}_{j,+}^{\dagger
m_{j}}|h_{0}\rangle,$ (S43)

with
$c_{\vec{m}}=\prod_{j}(4\lambda_{j})^{m_{j}}(1/2)_{m_{j}}/(\delta)_{\sum_{j}m_{j}}$.
From the above form of the purification, we can calculate the parametric form
of any normally-ordered correlation function:

$\displaystyle\langle\hat{b}_{1}^{\dagger n_{1}}$
$\displaystyle\cdots\hat{b}_{N}^{\dagger
n_{N}}\hat{b}_{1}^{m_{1}}\cdots\hat{b}_{N}^{m_{N}}\rangle=\frac{1}{\sqrt{2^{\sum_{j}n_{j}+m_{j}}}}\langle\Psi_{\hat{T}}|\hat{\beta}_{1,+}^{\dagger
n_{1}}\cdots\hat{\beta}_{N,+}^{\dagger
n_{N}}\hat{\beta}_{1,+}^{m_{1}}\cdots\hat{\beta}_{N,+}^{m_{N}}|\Psi_{\hat{T}}\rangle=\frac{1}{\sqrt{2^{\sum_{j}n_{j}+m_{j}}}}\sum_{\vec{k}\in\mathbb{N}^{N}}\frac{c^{*}_{\vec{n}+\vec{k}}c_{\vec{m}+\vec{k}}}{\vec{k}!}$
(S44)

There are a number of issues, however, with the series expression given above:
firstly, due to the weak symmetry $\hat{b}_{j}\to-\hat{b}_{j}$ of
$\mathcal{L}$, only correlation functions with $n_{j}\equiv m_{j}$ modulo two
are nonzero. Secondly, the series expression naively seems to be useless, as,
for a fixed total boson number cutoff, the series contains a number terms that
is exponentially growing with $N$. Therefore, the naive way of evaluating
(S44) scales no better than a direct simulation of the original master
equation.

We will resolve both of these issues now: first of all, we can efficiently
parametrize all of the nonzero correlators by replacing $\vec{n}\to
2\vec{n}+\vec{b},\vec{m}\to 2\vec{m}+\vec{b}$, where
$\vec{b}\in\mathbb{F}_{2}^{N}$ is a fixed vector of booleans. To exponentially
reduce the complexity of summing the series (S44), we define generalized
combinatorial form factors

$\displaystyle\Phi_{l}(\vec{\lambda},\vec{n},\vec{m},\vec{b};N):=\sum_{\sum_{j}k_{j}=l}\prod_{j}\frac{(1/2+n_{j}+b_{j})_{k_{j}}(1/2+m_{j}+b_{j})_{k_{j}}}{(2k_{j}+b_{j})!}(4\lambda_{j})^{2k_{j}}.$
(S45)

Indeed, in terms of the above form factors, the series expressions for
normally-ordered moments simplify quite considerably:

$\displaystyle\langle$ $\displaystyle\hat{b}^{\dagger
2\vec{n}+\vec{b}}\hat{b}^{2\vec{m}+\vec{b}}\rangle=\frac{1}{\sqrt{2^{\sum_{j}n_{j}+m_{j}+b_{j}}}}\sum_{\vec{k}\in\mathbb{N}^{N}}\frac{c^{*}_{2\vec{n}+\vec{b}+\vec{k}}c_{2\vec{m}+\vec{b}+\vec{k}}}{\vec{k}!}$
$\displaystyle=\frac{1}{\sqrt{2^{\sum_{j}n_{j}+m_{j}+b_{j}}}}\sum_{\vec{k}\in\mathbb{N}^{N}}\frac{c^{*}_{2(\vec{n}+\vec{b}+\vec{k})}c_{2(\vec{m}+\vec{b}+\vec{k})}}{(2\vec{k}+\vec{b})!}=\frac{c_{2(\vec{n}+\vec{b})}^{*}c_{2(\vec{m}+\vec{b})}}{\sqrt{2^{\sum_{j}n_{j}+m_{j}+b_{j}}}}\sum_{l=0}^{\infty}\frac{\Phi_{l}(\vec{\lambda},\vec{n},\vec{m},\vec{b})}{(\delta^{*}+\sum_{j}n_{j}+\sum_{j}b_{j})_{l}(\delta+\sum_{j}m_{j}+\sum_{j}b_{j})_{l}}.$
(S46)

Our job would not be finished if we did not give an efficient prescription for
evaluating the form factors $\Phi_{l}$. Our task is made easier, however, by
observing that, when $\vec{n}=\vec{m}=\vec{b}=\vec{0}$, the form factors (S45)
reduce to the form factor $\Phi_{n}(\vec{\lambda};N)$ used previously to
compute collective moments. In fact, these more general form factors satisfy
an analogous recursion relation:

$\displaystyle\Phi_{l}(\vec{\lambda},\vec{n},\vec{m},\vec{b};N)$
$\displaystyle=\sum_{p=0}^{n}\frac{(1/2+n_{N}+b_{N})_{p}(1/2+m_{N}+b_{N})_{p}}{(2p+b_{N})!}(4\lambda_{N})^{2p}\Phi_{n-p}(\vec{\lambda},\vec{n},\vec{m},;N-1).$
(S47)

We note that the above recursion relation has the boundary condition
$\Phi_{p}(\vec{\lambda},\vec{n},\vec{m},\vec{b};1)=(1/2+n_{1}+b_{1})_{p}(1/2+m_{1}+b_{1})_{p}(4\lambda_{1})^{2p}/(2p+b_{1})!$.
Therefore, the following $k$th-order approximant for each normally-ordered
moment,

$\displaystyle\frac{c_{2(\vec{n}+\vec{b})}^{*}c_{2(\vec{m}+\vec{b})}}{\sqrt{2^{\sum_{j}n_{j}+m_{j}+b_{j}}}}\sum_{l=0}^{k}\frac{\Phi_{l}(\vec{\lambda},\vec{n},\vec{m},\vec{b})}{(\delta^{*}+\sum_{j}n_{j}+\sum_{j}b_{j})_{l}(\delta+\sum_{j}m_{j}+\sum_{j}b_{j})_{l}},$
(S48)

may be evaluated using only $O(k^{2}N)$ operations. We now factor in
considerations as to the scaling of the total-particle-number cutoff $k$.
Assuming the onsite density $\bar{n}$ converges as $N\to\infty$, the total
number of particles is $O(N)$, and so the cutoff typically scales with the
system size: $k\sim O(N)$. Therefore, factoring in all considerations, time-
complexity of evaluating a single normally-ordered moment is roughly
$O(N^{3})$.

#### Quadratic correlation functions

Up until this point, we have given expressions for normally-ordered moments in
the basis of singular modes
$\hat{b}_{i}^{\dagger}:=V_{ij}\hat{a}_{j}^{\dagger}$. To obtain correlation
functions in the spatial mode basis, we must expand each $\hat{a}$-mode in
terms of $\hat{b}$-modes:

$\displaystyle\langle\hat{a}_{i}^{\dagger}\hat{a}_{i+r}\rangle$
$\displaystyle=\sum_{j}V_{i,j}V_{i+r,j}^{*}\langle\hat{b}_{j}^{\dagger}\hat{b}_{j}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\langle\hat{a}_{i}\hat{a}_{i+r}\rangle=\sum_{j}V_{i,j}^{*}V_{i+r,j}^{*}\langle\hat{b}_{j}^{2}\rangle$
(S49)

In both cases, one has to evaluate $N$ normally-ordered moments in the
$\hat{b}$-basis. By the previous estimates, each normally-ordered moment takes
$O(N^{3})$ floating-point operations to evaluate, so that in general it takes
$O(N^{4})$ floating-point operations to evaluate a quadratic correlation
function.

#### Quartic correlation functions

We now turn to evaluate quartic correlators using the same re-expansion
technique. We first compute the pair-correlation function:

$\displaystyle\langle\hat{a}_{i}^{\dagger
2}\hat{a}_{i+r}^{2}\rangle=\sum_{j_{1},j_{2},j_{3},j_{4}}V_{i,j_{1}}V_{i,j_{2}}V_{i+r,j_{3}}^{*}V_{i+r,j_{4}}^{*}\langle\hat{b}_{j_{1}}^{\dagger}\hat{b}_{j_{2}}^{\dagger}\hat{b}_{j_{3}}\hat{b}_{j_{4}}\rangle=2\sum_{j\neq
j^{\prime}}V_{i,j}V_{i,j^{\prime}}V_{i+r,j}^{*}V_{i+r,j^{\prime}}^{*}\langle\hat{b}_{j}^{\dagger}\hat{b}_{j^{\prime}}^{\dagger}\hat{b}_{j}\hat{b}_{j^{\prime}}\rangle+\sum_{j,j^{\prime}}V_{i,j}^{2}V_{i+r,j^{\prime}}^{*2}\langle\hat{b}_{j}^{\dagger
2}\hat{b}_{j^{\prime}}^{2}\rangle.$ (S50)

The calculation for the density-density correlation function (assuming $r\neq
0$) is very similar:

$\displaystyle\langle\hat{n}_{i}\hat{n}_{i+r}\rangle$
$\displaystyle=\sum_{j_{1},j_{2},j_{3},j_{4}}V_{i,j_{1}}V_{i+r,j_{2}}V_{i,j_{3}}^{*}V_{i+r,j_{4}}^{*}\langle\hat{b}_{j_{1}}^{\dagger}\hat{b}_{j_{2}}^{\dagger}\hat{b}_{j_{3}}\hat{b}_{j_{4}}\rangle$
$\displaystyle=\frac{1}{2}\sum_{j\neq
j^{\prime}}|V_{i,j}V_{i+r,j^{\prime}}+V_{i+r,j}V_{i,j^{\prime}}|^{2}\langle\hat{b}_{j}^{\dagger}\hat{b}_{j^{\prime}}^{\dagger}\hat{b}_{j}\hat{b}_{j^{\prime}}\rangle+\sum_{j,j^{\prime}}V_{i,j}V_{i+r,j}V_{i,j^{\prime}}^{*}V_{i+r,j^{\prime}}^{*}\langle\hat{b}_{j}^{\dagger
2}\hat{b}_{j^{\prime}}^{2}\rangle.$ (S51)

In both cases, one has to evaluate $O(N^{2})$ normally-ordered moments in the
$\hat{b}$-basis. By the previous estimates, each normally-ordered moment takes
$O(N^{3})$ floating-point operations to evaluate, so that in general it takes
$O(N^{5})$ floating-point operations to evaluate a quartic correlation
function.

##  V Phenomenology of the exact solution

###  V.1 $SU(1,1)$ coherent states

Figure S2: Bosonic pairing fluctuations near the PCS regime. (a) Here, we plot
the normalized fluctuations
$g^{(2)}_{K}:=(\langle\hat{k}_{-}^{\dagger}\hat{k}_{-}\rangle-|K|^{2})/|K|^{2}$
in the nonlocal pairing observable $K:=\langle\hat{k}_{-}\rangle$, for
different values of loss $\kappa\in\\{0.01U,0.1U,U\\}$ (red curves;
transparency increases with increasing loss). Note the sharp dip exactly at
$\Delta_{\text{PCS}}=U(2-N)/N$. Here, $G=U,\Lambda=0,N=500$. Inset: We plot
the normalized fluctuations $g^{(2)}_{\phi}:=(\langle\hat{a}^{\dagger
2}\hat{a}^{2}\rangle-|\phi|^{2})/|\phi|^{2}$ in the local onsite pairing
$\phi:=\langle\hat{a}_{j}^{2}\rangle$. Parameters same as before. (b) Same as
panel (a), but for $D=2$ with $N=8\times 8$ and periodic boundary conditions.
Here, $\Lambda=2U=4G$.

When $\delta=N/2$, the purification $|\Psi_{\hat{T}}\rangle$ is an
$SU(1,1)$-coherent state in the sense of Barut and Girardello (1971), that is,
an eigenstate of the lowering operator $\hat{K}_{-}$. This holds regardless of
whether the representation is unitary or not:

$\displaystyle\hat{K}_{-}|\Psi_{\hat{T}}\rangle=\sum_{m=0}^{\infty}\frac{(-1)^{m}}{(N/2)_{m}}\frac{\hat{K}_{-}\hat{K}_{+}^{m}}{m!}|h_{0}\rangle=\sum_{m=0}^{\infty}\frac{(-1)^{m+1}}{(N/2)_{m+1}}(m+1)(N/2+m)\frac{\hat{K}_{+}^{m}}{m!}|h_{0}\rangle=-|\Psi_{\hat{T}}\rangle.$
(S52)

The above identity has a dramatic effect on the physical steady state
$\hat{\rho}_{\text{ss}}$. In particular, fluctuations in
$K:=\langle\hat{k}_{-}\rangle$ exactly vanish when $\delta=N/2$, whereas
fluctuations in the onsite pairing $\phi:=\langle\hat{a}_{j}^{2}\rangle$
remain nonzero and show no special behavior. Here,
$\hat{k}_{-}:=\frac{1}{2}\sum_{ij}(uM^{-1})_{ij}\hat{a}_{i}\hat{a}_{j}$ is a
pair-lowering operator involving only the physical lattice modes
$\hat{a}_{j}$.

###  V.2 DMFT analysis of the $D=0$ model

#### DMFT expansion

We now use dynamical mean field theory (DMFT) to derive a large-$N$ expansion
for our $D=0$ model. To derive this expansion, we first fix a site (labelled
$j=0$) in our lattice model, and attempt to integrate-out the remaining
degrees of freedom. Typically this integration procedure is carried out within
the Schwinger-Keldysh formalism, which states that the action for the full
$D=0$ lattice model is

$\displaystyle S$ $\displaystyle=\int dt\sum_{\sigma=\pm
1}\sigma\bigg{[}\sum_{j}(\alpha_{j,\sigma}\partial_{t}\alpha_{j,\sigma}^{*}-\Delta
n_{j,\sigma}+G(\alpha_{j,\sigma}^{2}+c.c.))+\frac{U}{N}\bigg{(}\sum_{j}n_{j,\sigma}\bigg{)}^{2}\bigg{]}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ +i\kappa\int
dt\sum_{j}\bigg{(}\alpha_{j,+}\alpha_{j,-}-\frac{1}{2}\sum_{\sigma=\pm
1}|\alpha_{j,\sigma}|^{2}\bigg{)},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
n_{j,\sigma}:=|\alpha_{j,\sigma}|^{2}$ (S53)

where $\alpha_{j,\pm}=\alpha_{j,\pm}(t)$ are complex fields associated to the
bosonic degrees of freedom $\hat{a}_{j},\hat{a}_{j}^{\dagger}$, as is standard
in the Keldysh formalism. Within this formalism, the goal is now to compute
the effective action for a fixed site $j=0$. One can compute a large-$N$
asymptotic expansion for the effective action for the fields
$\alpha_{0,+},\alpha_{0,-}$, which takes the following self-consistent form at
leading-order:

$\displaystyle S_{\text{eff}}=S_{\text{free}}+\int dt\sum_{\sigma=\pm
1}2\sigma Un_{0,\sigma}\langle n_{0,\sigma}\rangle_{\text{eff}}+O(N^{-1}),$
(S54)

where $S_{\text{free}}$ is the Keldysh action for a single site, without the
Bose-Hubbard interaction, and $\langle\cdot\rangle_{\text{eff}}$ denotes an
average taken with respect to the effective action. Note that this leading-
order theory is Markovian. Therefore, as $N\to\infty$, one can evolve
observables for a fixed site using the self-consistent Lindbladian

$\displaystyle\mathcal{L}_{\text{eff}}\hat{\rho}_{0}=-i[\Delta_{\text{eff}}(\hat{\rho}_{0})\hat{n}_{0}+G\hat{a}_{0}^{\dagger
2}+h.c.,\hat{\rho}_{0}]+\kappa\mathcal{D}[\hat{a}_{0}]\hat{\rho}_{0},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Delta_{\text{eff}}(\hat{\rho}_{0})=\Delta+2U\text{Tr}[\hat{\rho}_{0}\hat{n}_{0}].$
(S55)

One can then solve for the steady-state density
$\bar{n}\equiv\bar{n}^{\text{MF}}$ within this leading-order mean-field
description, leading to a cubic equation for $\bar{n}^{\text{MF}}$. This leads
to a large-$N$ phase diagram where regions exist with up to three self-
consistent solutions for the density. We obtain good agreement between this
leading-order DMFT description and the exact solution, and find that

Figure S3: Benchmarking DMFT using the exact solution (in $D=0$). (a) Here, we
plot the mean onsite occupation $\bar{n}$ as a function of detuning, for
$G=U,\kappa=0.01U$. Note that we obtain asymptotic agreement with DMFT in the
limit that $N\to\infty$. To see work where a similar kind of mean field theory
was benchmarked by exact diagonalization results in a permutation-symmetric
spin model, see Wang and Fazio (2021). (b) Here, we plot the rms fluctuations
in $\hat{N}_{+}$ using the exact solution. We observe the empirical scaling
$\Delta
N_{+}:=\sqrt{\langle\hat{N}_{+}^{2}\rangle-\langle\hat{N}_{+}\rangle^{2}}=O(N^{1/2})$
so that $\Delta N_{+}/N$ vanishes as $N\to\infty$, ensuring the asymptotic
convergence of the wavefunction $|\Psi_{\hat{T}}\rangle$ of the paired boson
gas to the form predicted by DMFT. Here, $G=U/10,\kappa=U/100$, and
$\Delta=0$.

  * •

The location $\Delta_{c}$ of the 1st-order phase transition obtained from the
exact solution approaches the location of the bifurcation in the MFT cubic
self-consistency condition as $N\to\infty$. However, this convergence is
extremely slow (finite-size effects are still noticeable for $N\sim 10^{3}$).

  * •

The tristable region of the mean-field phase diagram has a unique point with
maximal loss rate $\kappa_{*}$. Let $(\Delta_{*},\kappa_{*})$ denote this
point. We call this point the mean-field critical-point. We find that the
mean-field critical point is precisely the same as the critical point in the
phase diagram obtained from the exact solution, that is,

$\kappa_{*}=\kappa_{c},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \Delta_{c}(\kappa_{c})=\Delta_{*}.$

#### Establishing the validity of DMFT

Normally DMFT expansions are hard to rigorously justify directly, even in the
limit of large coordination number $z\to\infty$. We will nonetheless give the
standard justification here, and then see how the exact hTRS solution yields a
much more direct perspective, at least with respect to the steady state
problem. The usual justification for the asymptotic result (S54) is as
follows: note that we can view our Hubbard interaction as a general (extended)
Bose-Hubbard interaction for a general graph $\mathcal{G}$ with
$\mathcal{G}=K_{N}$, i.e. a complete graph,

$\displaystyle\frac{U}{N}\bigg{(}\sum_{j}\hat{n}_{j}\bigg{)}^{2}=\frac{U}{z}\sum_{\langle
i,j\rangle\in\mathcal{G}}\hat{n}_{i}\hat{n}_{j}\Bigg{|}_{\mathcal{G}=K_{N}}.$
(S56)

The large-$N$ asymptotic result (S54) can be established, by copying exactly
the calculation in Scarlatella  _et al._ (2021); Strand  _et al._ (2015), but
by performing the cumulant expansion therein with respect to the density
fields $n_{j,\sigma}$ instead of the ordinary complex fields
$\alpha_{j,\sigma}$ (the calculation is relatively unilluminating, and for the
sake of brevity, we omit it here). This has the effect of establishing the
desired result for $\mathcal{G}$ a regular tree graph with coordination number
$N$, instead of a complete graph. After performing such a calculation, one
then waves one’s hands and claims that the same asymptotics (S54) holds on a
complete graph.

The exact solution yields a more direct way to check the validity of (S54), at
least from the perspective of the steady-state. The steady-state(s) predicted
by the leading-order DMFT dynamics $\mathcal{L}_{\text{eff}}$ admit the
following purification:

$\displaystyle|\Psi_{\text{DMFT}}\rangle=\sum_{m=0}^{\infty}\frac{1}{m!}\bigg{(}\frac{2\hat{K}_{+}}{2U\bar{n}^{\text{MF}}-\Delta_{\text{eff}}}\bigg{)}^{m}|h_{0}\rangle,$
(S57)

where $\bar{n}^{\text{MF}}$ is any solution to the cubic self-consistency
condition mentioned in the preceeding subsection. We can directly see that the
exact solution indeed converges to this form as $N\to\infty$. We begin by
noticing the asymptotics
$(N/2U)^{m}\Gamma(\delta)/\Gamma(m+\delta)\underset{N\to\infty}{\sim}(2Um/N-\Delta_{\text{eff}})^{-m}$
for the Gamma function, which yields the asymptotic estimate

$\displaystyle|\Psi_{\hat{T}}\rangle\underset{N\to\infty}{\sim}\sum_{m=0}^{\infty}\frac{1}{m!}\bigg{(}\frac{2\hat{K}_{+}}{2Um/N-\Delta_{\text{eff}}}\bigg{)}^{m}|h_{0}\rangle.$
(S58)

Notice that the above is identical to (S57), provided that we replace $m\sim
2\hat{N}_{+}$ with its average value. Indeed, fluctuations in $\hat{N}_{+}$
with respect to the exact solution $|\Psi_{\hat{T}}\rangle$ can be evaluated
directly via (S26) and shown to be of order $O(N^{1/2})$, so that $m/N$
becomes a deterministic variable in the large-$N$ limit. The series above thus
becomes well-concentrated about $m/N\sim\bar{n}^{\text{MF}}$, leading to the
self-consistent form (S57) predicted by leading-order DMFT.

#### Establishing the location of the critical point

We will now test the claim made previously, namely that the critical point
obtained from the exact solution lies exactly at the mean-field critical
point, i.e. we will compute the maximum magnitude of the susceptibility

$\displaystyle\chi_{\text{max}}(\kappa,N):=\sup_{\Delta}\Bigg{|}\frac{\partial\bar{n}}{\partial\Delta}\Bigg{|}$
(S59)

where all other parameters are held fixed. We confirm that $\chi_{\text{max}}$
diverges as $N\to\infty$ whenever $\kappa<\kappa_{*}$, and converges
otherwise. Repeatedly testing for the convergence of $\chi_{\text{max}}$ in
the thermodynamic limit for different values of $\kappa$ allows us to
approximately compute the rate at which $\chi_{\text{max}}$ diverges, thus
obtaining the critical exponent $\gamma$:

$\displaystyle\lim_{N\to\infty}\chi_{\text{max}}(\kappa,N)\underset{\kappa\to\kappa_{*}^{+}}{\sim}\tau^{-\gamma},$
(S60)

with $\tau:=(\kappa-\kappa_{*})/\kappa_{*}$. Using a very crude polynomial
fitting algorithm, we estimate $\gamma\approx-1$ (see Figure S4(b)).

Figure S4: Confirming the location of the critical point using the exact
solution. (a) Average density as a function of detuning $\Delta$ and loss
$\kappa$, with $N=500$, $\Lambda=0$, and $G=U$. Phase boundaries can be seen,
provided that $\kappa<\kappa_{*}\approx 4U$. (b) Maximum absolute value of the
susceptibility as a function of $\kappa$, for $\kappa\to\kappa_{*}^{+}$. Here,
$\Lambda=0$, and $G=U/10$. The polynomial fit used to estimate $\gamma$ is
depicted (dashed blue line).

###  V.3 Semiclassical limit

We now investigate our model without the $D\equiv 0$ restriction, but in the
semiclassical limit $\bar{n}\gg 1$. When the onsite photon occupation is
large, the dynamics of the field amplitudes
$\alpha_{j}(t):=\text{Tr}[\hat{\rho}(t)\hat{a}_{j}]$ is well captured by the
semiclassical equation of motion

$\displaystyle iu^{-1}\partial_{t}\beta_{j}$
$\displaystyle=2\lambda_{j}\beta^{*}_{j}+\beta_{j}\bigg{[}\bigg{(}2\sum_{k}|\beta_{k}|^{2}+1\bigg{)}-\Delta_{\text{eff}}/u\bigg{]},$
(S61)

where

$\displaystyle\beta_{j}^{*}\equiv\sum_{k=1}^{N}V_{j,k}\alpha_{k}^{*}$ (S62)

is the change-of-coordinates on the classical phase space induced by the
unitary $V$ in the Autonne-Takagi factorization $M/u=V\Sigma V^{T}$. We now
demonstrate the claim made in the main text, namely that the stable stationary
states of the above dynamical system are spheres in phase space formed by the
max-pairing modes. In particular, the semiclassical fixed-points satisfy the
equations

$\displaystyle
0=2\lambda_{j}\beta^{*}_{j}+\beta_{j}\bigg{[}\bigg{(}2\sum_{k}|\beta_{k}|^{2}+1\bigg{)}-\Delta_{\text{eff}}/u\bigg{]}$
(S63)

From now on, for clarity, we will use the symbol $X$ to denote the set of
fixed points. Note that, trivially, $\vec{0}\in X$. However, we are interested
in the nonzero fixed points. Therefore, let $\vec{\beta}_{\text{ss}}\in X$ be
a nonzero fixed point, i.e. a nonzero solution to (S63). Whenever
$\beta_{j,\text{ss}}\neq 0$, we can divide through by $\beta_{j}$ in (S63),
yielding a constraint on the phase
$e^{i\theta_{j}}\equiv\beta_{j,\text{ss}}/|\beta_{j,\text{ss}}|$:

$\displaystyle
e^{-2i\theta_{j}}=-\frac{2R_{\text{ss}}^{2}+1-\Delta_{\text{eff}}/u}{2\lambda_{j}},$
(S64)

where $R^{2}_{\text{ss}}\equiv\sum_{j}|\beta_{j,\text{ss}}|^{2}$. In
particular, if $\beta_{i,\text{ss}},\beta_{j,\text{ss}}\neq 0$ is any pair of
nonzero components of $\vec{\beta}_{\text{ss}}$, then, by taking the absolute
value of the above equation,

$\displaystyle\frac{1}{2\lambda_{i}}|2R^{2}_{\text{ss}}+1-\Delta_{\text{eff}}/u|=\frac{1}{2\lambda_{j}}|2R_{\text{ss}}^{2}+1-\Delta_{\text{eff}}/u|.$
(S65)

In particular, $\kappa\neq 0$ and so we can divide both sides by
$|2R^{2}_{\text{ss}}+1-\Delta_{\text{eff}}/u|$. Therefore,
$\lambda_{i}=\lambda_{j}$. We can go even further: (S64) also means that
$e^{2i\theta_{i}}=e^{2i\theta_{j}}$. In summary, for each
$\vec{\beta}_{\text{ss}}\in X$, there exists a $\lambda,\theta$ such that, for
all nonzero components of $\vec{\beta}_{\text{ss}}$,

$\displaystyle\lambda_{j}=\lambda,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\beta_{j,\text{ss}}=e^{i\theta}x_{j},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
x_{j}\in\mathbb{R},$ (S66)

where $\theta$ is independent of $j$ and uniquely determined by $\lambda$ in
the following way:

$\displaystyle
e^{-2i\theta}=-\frac{2R^{2}_{\text{ss}}+1-\Delta_{\text{eff}}/u}{2\lambda}.$
(S67)

By taking the real and imaginary parts of the above equation, we also have

$\displaystyle 2\lambda\sin 2\theta$ $\displaystyle=-\kappa/2u,$ (S68)
$\displaystyle 2\lambda\cos 2\theta$
$\displaystyle=-(2R_{\text{ss}}^{2}+1-\Delta/u).$ (S69)

#### Instability of solutions with $\lambda<\lambda_{*}$

Let $\vec{\beta}_{\text{ss}}\in X$ be a nonzero fixed point, and $\lambda$ the
corresponding singular value. Also, let $\lambda_{*}\equiv\sup_{j}\lambda_{j}$
denote the maximum singular value. We will now show that if
$\lambda\neq\lambda_{*}$, then $\vec{\beta}_{\text{ss}}$ is unstable. For this
purpose, it will be useful to rewrite the equations of motion in a coordinate
system $(\beta_{1}^{\prime},\dots,\beta_{N}^{\prime})$ that is adapted to
$\vec{\beta}_{\text{ss}}$, namely such that

$\displaystyle(\beta_{1,\text{ss}}^{\prime},\dots,\beta_{N,\text{ss}}^{\prime})=(e^{i\theta}R,0,\dots,0).$
(S70)

Crucially, by the arguments in the preceeding subsection, we can achieve this
via a rotation of the mode amplitudes $\beta_{j}$ with
$\lambda_{j}\equiv\lambda$:

$\displaystyle\beta_{j}^{\prime}=\begin{cases}\sum_{\lambda_{k}=\lambda}A_{j,k}\beta_{k}&\lambda_{j}=\lambda\\\
\beta_{j}&\lambda_{j}\neq\lambda\end{cases},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ A\in
O(s_{\lambda},\mathbb{R}),$ (S71)

where $s_{\lambda}$ denotes the multiplicity of the singular value $\lambda$.
Since the above transformation is a symmetry of the equations of motion (S61),
the equations of motion are covariant with respect to this transformation:

$\displaystyle iu^{-1}\partial_{t}\beta_{j}^{\prime}$
$\displaystyle=2\lambda_{j}(\beta_{j}^{\prime})^{*}+\beta_{j}^{\prime}\bigg{(}2\sum_{k}|\beta_{k}^{\prime}|^{2}+1-\Delta_{\text{eff}}/u\bigg{)}.$
(S72)

Now let
$\delta\beta^{\prime}_{j}\equiv\beta_{j}^{\prime}-\beta_{j,\text{ss}}^{\prime}$
denote the fluctuations about $\vec{\beta}_{\text{ss}}$. Assuming these
fluctuations are small, we can obtain linearized equations of motion for these
fluctuations:

$\displaystyle iu^{-1}\partial_{t}\delta\beta_{j}^{\prime}$
$\displaystyle=2\lambda_{j}\bigg{(}(\delta\beta_{j}^{\prime})^{*}-\frac{\lambda}{\lambda_{j}}e^{-2i\theta}\delta\beta_{j}^{\prime}\bigg{)}+2\beta_{j,\text{ss}}^{\prime}\sum_{k}((\beta_{k,\text{ss}}^{\prime})^{*}\delta\beta_{k}^{\prime}+c.c.)+O(\delta\beta^{\prime
2})$ (S73)

where we have implicitly used (S64). We now wish to argue that, if
$\lambda\neq\lambda_{*}$, then the Hurwitz criterion fails, that is, the
associated dynamical matrix contains an eigenvalue with positive real part. It
suffices to examine the stability of the fluctuations
$\delta\beta_{j}^{\prime}$ for $j\neq 1$, which evolve within this linear
approximation as follows:

$\displaystyle iu^{-1}\partial_{t}\delta\beta_{j}^{\prime}$
$\displaystyle=2\lambda_{j}\bigg{(}(\delta\beta_{j}^{\prime})^{*}-\frac{\lambda}{\lambda_{j}}e^{-2i\theta}\delta\beta_{j}^{\prime}\bigg{)}$
(S74)
$\displaystyle=2\lambda_{j}(\delta\beta_{j}^{\prime})^{*}-\frac{\kappa}{2}\delta\beta_{j}^{\prime}-i\bigg{(}\Delta-u(2R^{2}_{\text{ss}}+1)\bigg{)}\delta\beta_{j}^{\prime},$
(S75)

where we have shown explicitly that this is the equation of motion for a
detuned parametric amplifier, with detuning modified by the presence of the
Hubbard interaction $u$. Proceeding with the calculation, we then split the
fluctuations into real and imaginary parts via
$e^{-i\theta}\delta\beta_{j}^{\prime}=\delta x_{j}+i\delta y_{j}$, and obtain
the equations of motion

$\displaystyle u^{-1}\partial_{t}(\delta x_{j}+i\delta
y_{j})=-2\lambda_{j}e^{-2i\theta}\bigg{(}(1+\lambda/\lambda_{j})\delta
y_{j}+i(1-\lambda/\lambda_{j})\delta x_{j}\bigg{)}.$ (S76)

The corresponding eigenvalues of the dynamical matrix, using (S68), are

$\displaystyle\gamma_{j}^{\pm}=-\kappa/2\pm\sqrt{(\kappa/2)^{2}-4u^{2}(\lambda^{2}-\lambda_{j}^{2})},$
(S77)

Therefore, if $\lambda_{j}>\lambda$, then the eigenvalue $\gamma_{j}^{+}$ has
positive real part. Therefore, if $\lambda\neq\lambda_{*}$, then
$\vec{\beta}_{\text{ss}}$ is unstable. We also have a partial converse
statement: if $\lambda_{j}<\lambda$, then the fluctuations
$\delta\beta_{j}^{\prime}$ are stable.

#### Stability of solutions with $\lambda=\lambda_{*}$

Let $\vec{\beta}_{\text{ss}}\in X$ be a nonzero fixed point with corresponding
singular value $\lambda_{*}$. We will now investigate the conditions under
which $\vec{\beta}_{\text{ss}}$ is stable. By (S77), the fluctuations
$\delta\beta_{j}^{\prime}$ for $j=2,3,\dots N$ are all appropriately damped.
We thus must investigate the stability of the remaining fluctuations
$\delta\beta_{1}^{\prime}$, which have the following linearized equations of
motion:

$\displaystyle
u^{-1}\partial_{t}(e^{-i\theta}\delta\beta_{1}^{\prime})=-4\lambda_{*}e^{-2i\theta}\delta
y_{1}-4iR_{\text{ss}}^{2}\delta x_{1}.$ (S78)

The corresponding eigenvalues of the dynamical matrix, using (S68-S69), are

$\displaystyle\gamma^{\pm}_{j}=-(\kappa/2)\pm\sqrt{(\kappa/2)^{2}-8u^{2}R_{\text{ss}}^{2}(2R_{\text{ss}}^{2}+1)+8u\Delta
R_{\text{ss}}^{2}}.$ (S79)

Therefore, $\vec{\beta}_{\text{ss}}$ is stable if and only if

$\displaystyle 8u^{2}R_{\text{ss}}^{2}(2R_{\text{ss}}^{2}+1)+8u\Delta
R_{\text{ss}}^{2}>0,$ (S80)

which happens if and only if $u^{2}(2R_{\text{ss}}^{2}+1)-u\Delta>0$. To see
whether this criterion is satisfied, we must use the fact that
$\vec{\beta}_{\text{ss}}$ is a fixed point in order to obtain an additional
constraint on $R_{\text{ss}}$. In particular, taking the absolute value
squared of (S64) yields a quadratic equation for $R_{\text{ss}}^{2}+1$, with
two possible solutions:

$\displaystyle
2R_{\text{ss}}^{2}+1=\Delta/u\pm\sqrt{(2\lambda_{*})^{2}-(\kappa/2)^{2}}.$
(S81)

Therefore, criterion (S80) is satisfied if and only if

$\displaystyle
R_{\text{ss}}=\sqrt{\frac{\Delta/u-1+\sqrt{(2\lambda_{*})^{2}-(\kappa/2)^{2}}}{2}}>0.$
(S82)

Finally, let $\beta_{j}$ for $j=1,2,\dots,s$ denote the eigenmodes
corresponding to the maximum singular value $\lambda_{*}$, i.e. the so-called
"max-pairing modes". Since rotations $A\in O(s,\mathbb{R})$ of the max-pairing
modes are symmetries of the equations of motion, any such rotation $A$ must
send a stable fixed point to another stable fixed point. Therefore, when the
inequality (S82) is satisfied, the space $X^{\text{stab.}}\subset X$ of
nonzero stable fixed points is a sphere:

$\displaystyle
X^{\text{stab.}}=e^{i\theta}\bigg{\\{}x_{j}\in\mathbb{R}^{s}:\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\sum_{\lambda_{j}=\lambda_{*}}x_{j}^{2}=R_{\text{ss}}^{2}\bigg{\\}}$ (S83)

In particular, when $s>1$, the fluctuations tangent to $X^{\text{stab.}}$ are
Goldstone modes, i.e. zero-modes for the linearized dynamics. We can verify
this explicitly by fixing a point $\vec{\beta}_{\text{ss}}\in
X^{\text{stab.}}$, and expanding the linearized equations of motion for the
resulting fluctuations $e^{-i\theta}\delta\beta_{j}^{\prime}=\delta
x_{j}+i\delta y_{j}$, in the coordinate frame adapted to
$\vec{\beta}_{\text{ss}}$:

$\displaystyle iu^{-1}\partial_{t}(\delta x_{j}+i\delta
y_{j})=-4i\lambda_{*}e^{-2i\theta}\delta y_{j},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ j=2,\dots,s$
(S84)

In particular, the real-components $\delta x_{j}\in
T_{\vec{\beta}_{\text{ss}}}X^{\text{stab.}}$ of the fluctuations are zero
modes of the dynamical matrix, as was expected.

##  VI Mathematical background

###  VI.1 Proof of the $SU(1,1)$ decomposition theorem

We now establish the decomposition (S35). To make our task easier, we
establish the desired decomposition for a dense subspace of the Hilbert space.
The corresponding decomposition for the full Hilbert space can then be proven
using standard functional-analytic techniques.

In particular, let $V^{(j)}$ be the local nonunitary $SU(1,1)$ representations
defined in (S34). We define $W^{(j)}\subset V^{(j)}$ to be the algebraic part
of each representation, that is, the part of the representation consisting of
finite linear superpositions of Fock staes:

$\displaystyle W^{(j)}:=\bigg{\\{}\sum_{n}a_{n}\hat{\beta}_{j,+}^{\dagger
n}|h_{0}\rangle\in V^{(j)},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{finitely many\leavevmode\nobreak\
}a_{n}\leavevmode\nobreak\ \text{nonzero}\bigg{\\}}$ (S85)

The above subspaces allow us to conveniently establish the following theorem:

Theorem. The global nonunitary $SU(1,1)$ representation $\otimes_{j}W^{(j)}$
decomposes into irreducible subrepresentations as follows:

$\displaystyle\otimes_{j}W^{(j)}\simeq\bigoplus_{l=0}^{\infty}\Bigg{(}\bigoplus_{p=1}^{d_{l}}W_{l}^{(p)}\Bigg{)},$
(S86)

where the spaces $W^{(p)}_{l}\subset V^{(p)}_{l}$ are defined analogously:

$\displaystyle
W^{(p)}_{l}:=\bigg{\\{}\sum_{n}a_{n}\hat{K}_{+}^{n}|h_{l}^{(p)}\rangle,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{finitely
many\leavevmode\nobreak\ }a_{n}\leavevmode\nobreak\ \text{nonzero}\bigg{\\}}.$
(S87)

Proof. This can be proved by going to the Segal-Bargmann representation.
Within this representation, a finite-boson number state is represented as a
polynomial:

$\displaystyle\sum_{n}a_{n}\hat{\beta}_{j,+}^{\dagger n}|h_{0}\rangle\in
W^{(j)}\to\sum_{n}a_{n}x_{j}^{n}\in\mathbb{C}[x_{j}],$ (S88)

with $\mathbb{C}[x_{j}]$ the univariate polynomial ring generated by $x_{j}$.
It then follows that the tensor product $\otimes_{j}W^{(j)}$ is the
multivariate polynomial ring $\mathbb{C}[x_{1},\dots,x_{N}]$ generated by the
indeterminates $x_{1},\dots x_{N}$. Finally, creation and annihilation
operators are represented by partial differential operators:

$\displaystyle\hat{\beta}_{j,+}\to\frac{\partial}{\partial
x_{j}},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{\beta}_{j,+}^{\dagger}\to x_{j}.$ (S89)

We first establish our decomposition under the assumption that the singular
values $\lambda_{j}\equiv\lambda$ are all completely degenerate, and then
generalize the argument to the generic non-degenerate case
$\lambda_{i}\neq\lambda_{j}$. In the degenerate case, the global $SU(1,1)$
representation takes the simple form

$\displaystyle\hat{K}_{+}=\frac{\lambda}{2}\sum_{j}x_{j}^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{-}=\frac{1}{2\lambda}\sum_{j}\frac{\partial^{2}}{\partial x_{j}^{2}}$
(S90)

Therefore, within the Segal-Bargmann representation, our decomposition theorem
is equivalent to the following decomposition of the polynomial ring
$\otimes_{j}W^{(j)}$:

$\displaystyle\mathbb{C}[x_{1},...,x_{N}]\simeq\bigoplus_{l=0}^{\infty}\Bigg{(}\bigoplus_{p=1}^{d_{l}}\mathbb{C}[R^{2}]h_{l}^{(p)}\Bigg{)},$
(S91)

where $R^{2}=\sum_{j}x_{j}^{2}$, and $h_{l}^{(1)},...,h_{l}^{(d_{l})}$ is some
orthonormal basis of the space of harmonic homogeneous polynomials of degree
$l$. When interpreted pointwise, the above isomorphism reads as a harmonic
expansion of a fixed polynomial:

$\displaystyle
p(x_{1},...,x_{N})=\sum_{l=0}^{\infty}\sum_{p=1}^{d_{l}}q^{(p)}_{l}(R^{2})h_{l}^{(p)}(\vec{x})$
(S92)

where the $q^{(p)}_{l}$ are univariate polynomials. The decomposition (S91) is
then equivalent to the statement that the above mapping, interpreted as a
mapping from the RHS to the LHS, is an isomorphism of
$SU(1,1)$-representations. The representation (S90) acts the same on both
sides, so it suffices to establish the isomorphism at the vector space level,
i.e. prove that the above mapping is both injective and surjective. For
injectivity, it suffices to show that the subspaces
$\mathbb{C}[R^{2}]h_{l}^{(p)}$ are all mutually nonintersecting (except at
$\\{0\\}$), which can be verified by direct calculation:

$\displaystyle\langle
R^{2n}h_{l}^{(p)}|R^{2m}h_{l^{\prime}}^{(q)}\rangle=\bigg{(}\frac{2}{\lambda}\bigg{)}^{4n}\delta_{p,q}\delta_{l,l^{\prime}}\delta_{n,m}n!(N/2+l)_{n}$
(S93)

To demonstrate surjectivity, we must demonstrate that, for any polynomial $p$
on the LHS of (S92), a decomposition of the form given on the RHS exists. This
seems considerably more challenging to establish. However, here we are helped
by a basic fact from harmonic analysis:

Lemma. let $p$ denote a homogeneous polynomial of degree $l$. Then

$\displaystyle p(x_{1},...,x_{N})=h(\vec{x})+R^{2}q(x_{1},...,x_{N}),$ (S94)

where $q$ is a homogeneous polynomial of degree $l-2$, and $h$ is a
homogeneous harmonic polynomial of degree $l$.

A proof is usually given in standard textbooks on harmonic function theory
(see, e.g. Axler  _et al._ (2001)). In any case, iterating the above lemma, we
obtain, for any homogeneous polynomial $p$ a decomposition

$\displaystyle p(x_{1},...,x_{N})=\sum_{m=0,1,...}R^{2m}h_{\deg
p-2m}(\vec{x}),$ (S95)

where $h_{l}$ denotes a homogeneous harmonic polynomial of degree $l$, so that
the expansion (S92) can be established simply by writing out the polynomial on
the LHS of (S92) as a sum of homogeneous components, and then expanding the
harmonic polynomials appearing on the RHS of (S95) into a basis.

The preceeding arguments constitute a proof of the theorem for the degenerate
case $\lambda\equiv\lambda_{j}$. Therefore, all that is left is to reproduce
the above proof in the non-degenerate case $\lambda_{i}\neq\lambda_{j}$.
Luckily, the proof for the non-degenerate case follows immediately from the
degenerate case. In particular, we can write the $SU(1,1)$ representation as

$\displaystyle\hat{K}_{+}=\frac{1}{2}\sum_{j}y_{j}^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\hat{K}_{-}=\frac{1}{2}\sum_{j}\frac{\partial^{2}}{\partial y_{j}^{2}}$ (S96)

where we have made the change of variables $x_{j}\to
y_{j}:=\lambda_{j}^{1/2}x_{j}$. In particular, just by making the replacements
$x_{j}\to y_{j}$ in (S91), we obtain the desired result:

$\displaystyle\mathbb{C}[y_{1},...,y_{N}]\simeq\bigoplus_{l=0}^{\infty}\Bigg{(}\bigoplus_{p=1}^{d_{l}}\mathbb{C}[R^{2}]h_{l}^{(p)}\Bigg{)},$
(S97)

with $R^{2}$ and $h_{l}^{(p)}$ defined just as in (S91), but with the
replacements $x_{j}\to y_{j}$. This completes the proof of our theorem in the
non-degenerate case.

###  VI.2 Exact solution for the steady-state Wigner function

We now compute a closed form for the Wigner function $W_{\text{ss}}$ of the
steady-state density matrix $\hat{\rho}_{\text{ss}}$. The arguments in Roberts
and Clerk (2020) generalize in a straightforward manner to the $N$-mode case.
In particular, the steady-state Wigner function for our system satisfies the
identity

$\displaystyle W_{\text{ss}}(\vec{\alpha})$
$\displaystyle=2^{N}Q_{|\Psi_{+}\rangle}\bigg{(}\sqrt{2}\vec{\alpha}\bigg{)},$
(S98)

where $Q_{|\Psi_{+}\rangle}$ is the Husimi-Q representation of the $N$-mode
pure state $|\Psi_{+}\rangle$. Therefore, to express the Wigner function
$W_{\text{ss}}$ of the steady state in closed form it suffices to express the
Husimi-Q representation of $|\Psi_{+}\rangle$ in closed form. This task is
solved easily via the Segal-Bargmann representation Bargmann (1961, 1967) of
$|\Psi_{+}\rangle$, which can be calculated in a straightforward manner from
(S18):

$\displaystyle\Psi_{\text{SB}}(\vec{\alpha})=\frac{\,{}_{0}F_{1}\big{(}\delta;-\frac{1}{2}\sum_{ij}\frac{M_{ij}}{u}\alpha_{i}\alpha_{j}\big{)}}{\sqrt{\mathcal{N}}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\mathcal{N}=\sum_{l=0}^{\infty}\frac{\Phi_{l}(\vec{\lambda};N)}{(\delta)_{l}(\delta^{*})_{l}}.$
(S99)

Now, the $Q$-function of a pure state $|\Psi\rangle$ is given by
$Q_{|\Psi\rangle}(\vec{\alpha})=\pi^{-N}|\Psi_{\text{SB}}(\vec{\alpha}^{*})|^{2}e^{-|\vec{\alpha}|^{2}}$,
where $\Psi_{\text{SB}}$ is the Segal-Bargmann representation of
$|\Psi\rangle$. Therefore, (S98) yields

$\displaystyle W_{\text{ss}}(\vec{\alpha})$
$\displaystyle=\frac{1}{\mathcal{N}}\bigg{(}\frac{2}{\pi}\bigg{)}^{N}\bigg{|}\,_{0}F_{1}\bigg{(}\delta;-\sum_{ij}\frac{M_{ij}}{u}\alpha_{i}^{*}\alpha_{j}^{*}\bigg{)}\bigg{|}^{2}e^{-2|\vec{\alpha}|^{2}},$
(S100)

Figure S5: Nonthermal nature of the steady state. The exact solution for the
Wigner function can be used to symbolically check that the Hamiltonian and
steady state do not commute. Here, we evaluate the phase-space commutator
$H\star W_{\text{ss}}-W_{\text{ss}}\star H$ in the case $G=\kappa=U,$ and
$\Lambda=0$, for the cases (a) $N=1$, in which case we choose to evaluate the
result at the phase space point $\alpha=1$ and (b) $N=2$, in which case we
choose to evaluate the result at the phase space point
$\alpha_{1}=\alpha_{2}=1$. Note that, since $H,W_{\text{ss}}$ are both purely
real, the result is always purely imaginary. As expected, the result always
vanishes in the limit $\kappa\to 0^{+}$.

#### Scaling limit for the Wigner function

We now investigate the expression for the Wigner function in the high-density
limit $M_{ij}/u\to\infty$. Note that the expression (S100) for the Wigner
function is non-negative, and thus can be interpreted as a bona-fide
probability measure. In the limit $M_{ij}/u\to\infty$, this probability
measure will become supported on larger and larger regions of phase space, and
so it is useful to re-scale the phase space in such a way that the resulting
rescaled distribution converges to a limit.

The correct scaling turns out to be
$\beta_{j}:=\sqrt{-\lambda_{*}}\tilde{\beta}_{j}$, where
$\lambda_{*}\equiv\sup_{j}\lambda_{j}$ is the maximum singular value appearing
in the Autonne-Takagi factorization $M/u=V\Sigma V^{T}$, and

$\displaystyle\beta_{j}^{*}\equiv\sum_{k=1}^{N}V_{j,k}\alpha_{k}^{*}$ (S101)

is the change-of-coordinates on the classical phase space of our system,
induced by the unitary $V$ appearing in the factorization. The explicit form
of $V$ can be recovered from the spectral decomposition of the Laplacian of
our underlying connectivity graph. One can show that, at least when
$\Delta=0$, $\kappa=0^{+}$, in which case the Bessel function sitting inside
the absolute value in (S100) becomes a hyperbolic cosine, the steady-state
Wigner function $W_{\text{ss}}(\tilde{\beta})$ limits to a uniform
distribution on the sphere

$\displaystyle
S\equiv\bigg{\\{}(\tilde{\beta}_{1},\dots,\tilde{\beta}_{s},0,\dots,0)\in\mathbb{R}^{N}:\leavevmode\nobreak\
\leavevmode\nobreak\ \sum_{j=1}^{s}\tilde{\beta}_{j}^{2}=1\bigg{\\}},$ (S102)

where $\beta_{1},\dots,\beta_{s}$ are the max-pairing modes. This can be
established by expanding the Wigner function (S100) into a sum of four
exponentials, and then solving the associated saddle-point equations.

###  VI.3 Using the Wigner function to verify the nonthermal character of the
steady state

The exact solution for the Wigner function can be used to explicitly showcase
the nonequilibrium character of the steady state. One non-thermal feature of
the steady state is that it does not commute with the Hamiltonian, that is,
$[H,\rho_{\text{ss}}]\neq 0$. Therefore, the steady state cannot be written as
$\exp(-\beta\hat{H})$ for some $\beta$. This can be verified efficiently when
$\Lambda=0$, in which case we can write down the closed-form solution

$\displaystyle W_{\text{ss}}(\vec{\alpha})$
$\displaystyle=\bigg{(}\frac{2}{\pi}\bigg{)}^{N}\frac{\big{|}\,_{0}F_{1}\big{(}\delta;-2\lambda\sum_{j}\alpha_{j}^{*2}\big{)}\big{|}^{2}}{\,{}_{1}F_{2}(N/2;\delta,\delta^{*};\lambda^{2})}e^{-2|\vec{\alpha}|^{2}},$
(S103)

with $\lambda=NG/U$ in this case corresponding to the unique singular value of
the pairing matrix. To show that the steady state and the Hamiltonian do not
commute, we pass to the phase-space formulation of quantum mechanics. In the
phase space formulation of quantum mechanics, the noncommutativity of the
steady state and the Hamiltonian is equivalent to the statement that

$\displaystyle W_{\text{ss}}\star H-H\star W_{\text{ss}}\neq 0,$ (S104)

where $\star$ is the Moyal product Moyal (1949), and $H$ denotes the
symmetrically-ordered (i.e. Weyl) symbol of the Hamiltonian. Because $H$ is a
polynomial, the following derivative expansion terminates at a finite order
and hence can be calculated symbolically in closed form using a simple
computer algebra program:

$\displaystyle(f\star
g)(\vec{\alpha})=\exp\bigg{(}-\frac{1}{2}\sum_{j=1}^{N}\bigg{(}\frac{\partial}{\partial
x_{j}^{*}}\frac{\partial}{\partial y_{j}}-\frac{\partial}{\partial
x_{j}}\frac{\partial}{\partial
y_{j}^{*}}\bigg{)}\bigg{)}f(\vec{x})g(\vec{y})\bigg{|}_{\vec{x}=\vec{y}=\vec{\alpha}}$
(S105)

Since both the Wigner function and Hamiltonian are generically smooth
functions of $\vec{\alpha}$, the phase-space function S104 is generically a
smooth function of $\vec{\alpha}$. In Fig. S5 we exhibit a single point where
this phase-space function is non-vanishing.

##  VII Experimental realization using superconducting circuits ($D=0$)

Figure S6: First iteration of the circuit (no driving). A Josephson tunnel
junction is placed in parallel with a chain of LC oscillators to provide a
global Hubbard interaction.

Our $D=0$ model is relatively easily realizable using a simple superconducting
circuit with only three nonlinear elements. To see this, we first attempt to
realize the $D=0$ Hamiltonian without the driving. We begin by placing a chain
of $N$ uncoupled LC oscillators in series. Via examination of Kirchoff’s
current laws, the Hamiltonian that describes the equations of motion for the
chain is the following:

$\displaystyle
H=\sum_{j}\bigg{(}\frac{Q_{j}^{2}}{2C_{j}}+\frac{\Phi_{j}^{2}}{2L_{j}}\bigg{)},$
(S106)

where here, $Q_{j}$ denotes the charge stored on the $j$th capacitor, and
$\Phi_{j}$ denotes the time-integrated voltage across each inductor. One may
diagonalize this Hamiltonian by defining dimensionless creation- and
annihilation operators

$\displaystyle\hat{\Phi}_{j}:=\Phi_{j}^{\text{zpf}}(\hat{a}_{j}+\hat{a}_{j}^{\dagger}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{Q}_{j}:=-iQ_{j}^{\text{zpf}}(\hat{a}_{j}-\hat{a}_{j}^{\dagger}),$ (S107)

where here, $Q_{j}^{\text{zpf}},\Phi_{j}^{\text{zpf}}$ are the zero-point
vacuum fluctuations of the charge and phase across the inductive and
capacitive branches of each LC oscillator.

To add an infinite-range Bose-Hubbard interaction, we place a Josephson
junction in parallel with the chain of oscillators (c.f. Figure S6). The
Josephson junction, in the limit of extremely weak junction capacitance
$C_{J}\ll C_{j}$, can be modelled accurately to leading order via the
following interaction Hamiltonian:

$\displaystyle\hat{H}_{\text{int}}=E_{J}\cos\bigg{(}\sum_{j}\varphi_{j}(\hat{a}_{j}+\hat{a}_{j}^{\dagger})\bigg{)},$
(S108)

where $\varphi_{j}:=\Phi_{j}^{\text{zpf}}/2\pi\Phi_{0}$. We now tune the
dimensionless phase fluctuations $\varphi_{j}$ to be parametrically small and
uniform across all modes, that is $\varphi_{j}\equiv\varphi\ll 1$. Note that
this can be done without constraining the resonant frequencies of each
oscillator. Taylor’s theorem then says that

$\displaystyle\hat{H}_{\text{int}}=-\frac{E_{J}\varphi^{2}}{\hbar
2!}\bigg{(}\sum_{j}\hat{a}_{j}+h.c.\bigg{)}^{2}+\frac{E_{J}\varphi^{4}}{\hbar
4!}\bigg{(}\sum_{j}\hat{a}_{j}+h.c.\bigg{)}^{4}+O(\varphi^{6}).$ (S109)

We then go into a rotating frame with respect to the free Hamiltonian
$H_{\text{free}}:=\hbar\sum_{j}\omega_{j}\hat{n}_{j}$, where $\omega_{j}$ is
the bare resonance frequency of each LC resonator in the chain. Now we choose
the resonant frequencies of each mode so that the fundamental frequency
$\Omega$ of the rotating-frame Hamiltonian is much larger than the rate
$\hbar^{-1}E_{J}\varphi^{2}$. In this regime, the rotating-wave approximation
is valid and yields

$\displaystyle\hat{H}_{\text{RWA}}=-E_{J}\varphi^{2}\sum_{j}\hat{n}_{j}+\frac{E_{J}\varphi^{4}}{2}\bigg{(}\sum_{j}\hat{n}_{j}\bigg{)}^{2}.$
(S110)

###  VII.1 Adding a two-photon drive

Figure S7: Circuit incorporating coherent two-photon driving. The global
parametric drive is supplied by a flux-tunable transmon (blue shaded region).

To obtain the $D=0$ Hamiltonian for our model, we just have to add two-photon
driving to the above scheme. To do this, we play the same trick as above: this
time we add a symmetric SQUID in parallel with the oscillator chain. Assuming
that the junction capacitances in the SQUID are also much smaller than the
capacitances present in the oscillator chain, the new interaction Hamiltonian
is simply

$\displaystyle\hat{H}_{\text{int}}=E_{L}\cos\bigg{(}\varphi\sum_{j}\hat{a}_{j}+h.c.\bigg{)}+2E_{R}\cos\bigg{(}\frac{\Phi_{e}}{2\pi\Phi_{0}}\bigg{)}\cos\bigg{(}\varphi\sum_{j}\hat{a}_{j}+h.c.\bigg{)},$
(S111)

where $E_{R}=E_{R,1},E_{R,2}$ is the Josephson energy of each junction in the
symmetric SQUID. We then choose to drive the SQUID in such a way that

$\displaystyle\frac{\Phi_{e}}{2\pi\Phi_{0}}=\frac{\pi}{2}-\epsilon_{p}(t).$
(S112)

We also assume that $E_{R}\ll E_{L}$. As a result, we can truncate the
expansion of the SQUID potential at quadratic order, while continuing to
truncate the expansion of the left junction at quartic order:

$\displaystyle\hat{H}_{\text{int}}=-E_{R}\varphi^{2}\epsilon_{p}(t)\bigg{(}\sum_{j}\hat{a}_{j}+h.c.\bigg{)}^{2}-\frac{E_{L}\varphi^{2}}{2!}\bigg{(}\sum_{j}\hat{a}_{j}+h.c.\bigg{)}^{2}+\frac{E_{L}\varphi^{4}}{4!}\bigg{(}\sum_{j}\hat{a}_{j}+h.c.\bigg{)}^{4}+O(E_{L}\varphi^{6})+O(E_{R}\varphi^{4}).$
(S113)

By modulating the pump amplitude appropriately via

$\displaystyle\epsilon_{p}(t):=\epsilon_{0}\sum_{j}\cos(2\omega_{j}-2\omega_{p})t,$
(S114)

and going into a rotating frame with respect to the free Hamiltonian
$\hat{H}_{\text{free}}=\hbar\sum_{j}(\omega_{j}-\omega_{p})\hat{n}_{j}$, and
again assuming that the mode frequencies are all appropriately detuned from
each other, we obtain the following rotating-wave Hamiltonian:

$\displaystyle\hat{H}_{\text{RWA}}=\frac{E_{L}\varphi^{4}}{2}\bigg{(}\sum_{j}\hat{n}_{j}\bigg{)}^{2}+\sum_{j}(\hbar\omega_{p}-E_{L}\varphi^{2})\hat{n}_{j}-E_{R}\varphi^{2}\epsilon_{0}\sum_{j}(\hat{a}_{j}^{2}+h.c.)$
(S115)

We thus obtain the exact parameters of our solvable model in the regime $D=0$
(in SI units!):

$\displaystyle U$
$\displaystyle=\frac{NE_{L}\varphi^{4}}{2\hbar},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
G=-\frac{E_{R}\varphi^{2}\epsilon_{0}}{\hbar},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Delta=\omega_{p}-\frac{E_{L}\varphi^{2}}{\hbar},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\varphi=\sqrt{\frac{\hbar}{2}}\frac{(L_{j}/C_{j})^{1/4}}{2\pi\Phi_{0}}\equiv\text{const.}$
(S116)

###  VII.2 Effect of junction capacitances

We now compute the corrections to the Hamiltonian due to junction
capacitances, and demonstrate that these capacitances can be neglected when
they are much smaller than the capacitances in the oscillator chain. To
simplify the analysis we assume the junction capacitances in the symmetric
SQUID are the same, and define a new parameter $C_{\Sigma}:=C_{L}+2C_{R}$
corresponding to the sum of the three junction capacitances. The Maxwell
capacitance relation of the circuit, assuming the junction capacitances are
all zero, is

$\displaystyle\vec{q}=C_{0}\dot{\vec{\phi}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ C_{0}:=\begin{pmatrix}C_{1}&-C_{2}&0&\cdots&0\\\
-C_{2}&C_{2}+C_{3}&-C_{3}&&\vdots\\\ 0&-C_{3}&C_{3}+C_{4}&-C_{4}&\\\
\vdots&&-C_{4}&\ddots&\\\ &&&&-C_{N-1}\\\
0&\cdots&&-C_{N-1}&C_{N}\end{pmatrix},$ (S117)

where the nodal charges $q_{j}$ can be expressed in terms of the charges
$Q_{j}$ in the capacitance chain via

$\displaystyle Q_{N-j}=\sum_{k=0}^{j}q_{N-k},$ (S118)

With the junction capacitances included, the new capacitance matrix is
obtained in a very simple manner from $C_{0}$:

$\displaystyle C=C_{0}+\begin{pmatrix}0&\cdots&0\\\ \vdots&&\vdots\\\
0&\cdots&C_{\Sigma}\end{pmatrix}.$ (S119)

To obtain the corrections to the Hamiltonian that were neglected in the
previous analysis, we use the Sherman-Morrison formula, which is exact:

$\displaystyle
C_{k,k}^{-1}=(C_{0}^{-1})_{k,k^{\prime}}+\frac{C_{\Sigma}(C_{0}^{-1})_{k,N}(C_{0}^{-1})_{N,k^{\prime}}}{1+C_{\Sigma}(C_{0}^{-1})_{N,N}},$
(S120)

so that the correction to the Hamiltonian is rigorously

$\displaystyle\delta\hat{H}=\frac{C_{\Sigma}}{2}\sum_{k,k^{\prime}}\frac{(C_{0}^{-1})_{k,N}(C_{0}^{-1})_{N,k^{\prime}}}{1+C_{\Sigma}(C_{0}^{-1})_{N,N}}q_{k}q_{k^{\prime}},$
(S121)

which goes to zero as $C_{\Sigma}(C_{0}^{-1})_{i,j}\to 0$, as expected.

## References

  * Buča and Prosen (2012) B. Buča and T. Prosen, New Journal of Physics 14, 073007 (2012).
  * Fagnola and Umanità (2010) F. Fagnola and V. Umanità, Communications in Mathematical Physics 298, 523 (2010).
  * Autonne (1915) L. Autonne, _Sur les matrices hypohermitiennes et sur les matrices unitaires,_ (A. Rey, Lyon, 1915) oCLC: 3213505.
  * Bargmann (1961) V. Bargmann, Communications on Pure and Applied Mathematics 14, 187 (1961).
  * Bargmann (1967) V. Bargmann, Communications on Pure and Applied Mathematics 20, 1 (1967).
  * Axler (2020) S. Axler, “HFT.m,”  (2020).
  * Barut and Girardello (1971) A. O. Barut and L. Girardello, Communications in Mathematical Physics 21, 41 (1971).
  * Wang and Fazio (2021) P. Wang and R. Fazio, Phys. Rev. A 103, 013306 (2021).
  * Scarlatella  _et al._ (2021) O. Scarlatella, A. A. Clerk, R. Fazio, and M. Schiró, Physical Review X 11, 031018 (2021), publisher: American Physical Society.
  * Strand  _et al._ (2015) H. U. Strand, M. Eckstein, and P. Werner, Physical Review X 5, 011038 (2015), publisher: American Physical Society.
  * Axler  _et al._ (2001) S. Axler, P. Bourdon, and W. Ramey,  _Harmonic Function Theory_, 2nd ed., Graduate Texts in Mathematics, Vol. 137 (Springer, 2001).
  * Roberts and Clerk (2020) D. Roberts and A. A. Clerk, Phys. Rev. X 10, 021022 (2020).
  * Moyal (1949) J. E. Moyal, Mathematical Proceedings of the Cambridge Philosophical Society 45, 99–124 (1949).