#  Cavity-based reservoir engineering for Floquet-engineered superconducting
circuits

Francesco Petiziol<EMAIL_ADDRESS>Technische Universität Berlin,
Institut für Theoretische Physik, Hardenbergstraße 36, Berlin 10623, Germany
André Eckardt  Technische Universität Berlin, Institut für Theoretische
Physik, Hardenbergstraße 36, Berlin 10623, Germany

###### Abstract

Considering the example of superconducting circuits, we show how Floquet
engineering can be combined with reservoir engineering for the controlled
preparation of target states. Floquet engineering refers to the control of a
quantum system by means of time-periodic forcing, typically in the high-
frequency regime, so that the system is governed effectively by a time-
independent Floquet Hamiltonian with novel interesting properties. Reservoir
engineering, on the other hand, can be achieved in superconducting circuits by
coupling a system of artificial atoms (or qubits) dispersively to pumped leaky
cavities, so that the induced dissipation guides the system into a desired
target state. It is not obvious that the two approaches can be combined, since
reaching the dispersive regime, in which system and cavities exchange
excitations only virtually, can be spoiled by driving-induced resonant
transitions. However, working in the extended Floquet space and treating both
system-cavity coupling as well as driving-induced excitation processes on the
same footing perturbatively, we identify regimes, where reservoir engineering
of targeted Floquet states is possible and accurately described by an
effective time-independent master equation. We successfully benchmark our
approach for the preparation of the ground state in a system of interacting
bosons subjected to Floquet engineered magnetic fields in different lattice
geometries.

††preprint: APS/123-QED Figure 1: $M$ artificial atoms are coupled to $L$
pumped and leaky cavities. The atoms are periodically driven to Floquet
engineer effective Hamiltonians with desired properties, while the cavities
induce dissipation in a controlled way.

Floquet engineering is a powerful tool for quantum simulation, where time-
periodic driving is applied to manipulate the properties of a quantum system
[1, 2, 3, 4, 5, 6]. It has been applied successfully to engineered quantum
systems, such as ultracold atoms in optical lattices [7, 8, 9, 10, 11, 12, 13,
14, 15], photons in optical waveguides [16, 17, 18], and superconducting
circuits [19]. Recently, the question has been addressed, whether it is
possible to prepare ground (or low-temperature) states of effective Floquet-
engineered Hamiltonians by coupling them to a thermal environment [20, 21, 22,
23, 24, 25, 26]. Here we propose an alternative strategy for dissipative state
preparation in Floquet systems based on cavity-assisted reservoir engineering,
as it can be realized in superconducting circuits by coupling artificial atoms
to pumped leaky resonators. Such an approach is not straightforward, since
cavity-based reservoir engineering relies on the so-called dispersive regime,
in which the system exchanges excitations with the cavities only virtually.
Thus, one has to identify a regime, in which such a virtual change is not
spoiled by driving-induced resonant processes. Understanding and avoiding the
breakdown of dispersive regimes in circuit QED systems under periodic
modulation is a central challenge, as investigated, e.g., in the realization
of fast gates and controlled nonlinearities for high-Q cavity modes coupled
via driven transmons [27, 28]. To solve this problem, we describe both atoms
and cavities in the extended Floquet space. In this framework, we treat both
the system-bath coupling as well as driving-induced excitation processes on
equal footing using degenerate perturbation theory. This combined approach
contains both the standard perturbative treatment of the dispersive coupling
in non-driven systems as well as the high-frequency expansion of isolated
periodically-driven systems as limiting cases. But it also includes the
interplay of both processes, which has a crucial (potentially detrimental)
impact on the open driven dynamics and the preparation of target states. Based
on this theory, we formulate driven-dissipative schemes for the preparation of
non-trivial states in finite Floquet-engineered flux ladders (exhibiting
chiral ground state currents and frustration-induced localization effects).
The approach is confirmed via simulations of the full driven-dissipative
evolution of atoms and cavities.

We consider a two-dimensional array of $M$ artificial atoms (Fig. 1) in a
superconducting circuit [29, 30, 31, 19, 32] described by the Bose-Hubbard
Hamiltonian

$\displaystyle\hat{H}_{\mathrm{S}}(t)=\frac{U}{2}\sum_{j=1}^{M}\hat{n}_{j}(\hat{n}_{j}-1)-J\sum_{\langle
j,j^{\prime}\rangle}e^{i\theta_{j^{\prime}j}(t)}\hat{a}_{j^{\prime}}^{\dagger}\hat{a}_{j}.$
(1)

Here $\hat{a}_{j}$ and $\hat{n}_{j}=\hat{a}^{\dagger}_{j}\hat{a}_{j}$ denote
the bosonic annihilation and number operators for an excitation on site $j$.
The excitations experience an attractive on-site potential $U<0$,
corresponding to a level anharmonicity. Moreover, they can tunnel between
neighbouring sites with matrix elements of amplitude $J>0$. The time-periodic
Peierls phases $\theta_{ij}(t)=\theta_{i}(t)-\theta_{j}(t)$ describe a time-
dependent force, which, in a non-rotating reference frame, is described by on-
site potentials
$v_{j}=\hbar\dot{\theta}_{j}=\Delta+\nu_{j}\hbar\omega+\lambda\sin(\omega
t-\varphi_{j})$ with integer $\nu_{j}$. The transition to the rotating frame
adopted in Eq. (1) and considered in the following, is accomplished by
replacing $\hat{a}_{j}\to e^{-i\theta_{j}(t)}\hat{a}_{j}$ with
$\theta_{j}(t)=\Delta t/\hbar+\nu_{j}\omega t-\lambda\cos(\omega
t-\varphi_{j})/\hbar\omega$.

For isolated systems, the motivation for applying such periodic forcing is
that in the high-frequency regime, $\hbar\omega\gg|U|,J$ the dynamics is
approximately described by an effective time-independent Hamiltonian
$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}$, with new properties. In leading order,
it is obtained as time average
$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}=\frac{1}{T}\int_{0}^{T}\\!dt\hat{H}_{\mathrm{S}}(t)$
over one driving period $T=2\pi/\omega$ (rotating-wave approximation). This
gives rise to effective tunneling matrix elements
$-J_{jj^{\prime}}^{\mathrm{eff}}e^{i\theta^{\mathrm{eff}}_{jj^{\prime}}}$,
with amplitude
$J_{jj^{\prime}}^{\mathrm{eff}}=J\mathcal{J}_{\nu_{j}-\nu_{j^{\prime}}}(2\lambda\sin[(\varphi_{j^{\prime}}-\varphi_{j})/2]/\hbar\omega)$,
where $\mathcal{J}_{n}(\cdot)$ denotes a Bessel function, and with Peierls
phases
$\theta^{\mathrm{eff}}_{jj^{\prime}}=(\nu_{j^{\prime}}-\nu_{j})(\varphi_{j}+\varphi_{j^{\prime}})/2$,
which can describe artificial magnetic fields [3]. Such (and similar) Floquet
engineering has been employed successfully to experimentally engineer and
study interaction-driven phase transitions [8], kinetic frustration [9],
topological band structures [11, 12, 15], their chiral edge modes [16, 17,
19], Aharonov-Bohm cages [33], and two qubit-gates [34, 32] in systems of
ultracold atoms in optical lattices, optical waveguides, and superconducting
circuits.

We will investigate, whether it is possible to employ reservoir engineering
for cooling the system into the ground state of
$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}$. To this end, some of the atoms shall be
coupled individually to driven-damped cavities. The open dynamics of the whole
system is described by the master equation
$d\hat{\rho}/dt=-i[\hat{H}(t),\hat{\rho}]/\hbar+\sum_{\ell=1}^{L}\kappa_{\ell}\mathcal{D}[\hat{c}_{\ell}]\hat{\rho},$
where
$\mathcal{D}[\hat{c}_{\ell}]\hat{\rho}=\hat{c}_{\ell}\hat{\rho}\hat{c}_{\ell}^{\dagger}-\frac{1}{2}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}$
is a Lindblad dissipator and where the Hamiltonian
$\hat{H}(t)=\hat{H}_{\mathrm{S}}(t)+\hat{H}_{\mathrm{SC}}(t)+\hat{H}_{\mathrm{C}}(t)$
comprises the terms
$\hat{H}_{\mathrm{C}}(t)=\sum_{j=1}^{L}\big{[}\delta_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\mathcal{E}_{j}\hat{c}_{j}^{\dagger}e^{-i\omega_{j}t}+\mathcal{E}_{j}^{*}\hat{c}_{j}e^{i\omega_{j}t}]$
and
$\hat{H}_{\mathrm{SC}}(t)=\sum_{j=1}^{L}g_{j}\big{[}e^{-i\theta_{j}(t)}\hat{a}_{j}\hat{c}_{j}^{\dagger}+\mathrm{H.c.}]$
that describe the cavities in a frame rotating at frequency $\Delta/\hbar$ and
their coupling to the system [35, 32]. The $L\leq M$ cavities are described by
bosonic annihilation operators $\hat{c}_{j}$, are detuned by the atoms by
$\delta_{j}$, pumped with strength $\mathcal{E}_{j}$ at a frequency
$\omega_{j}$, and leak photons at a rate $\kappa_{j}$. The atoms are
enumerated so that the $j$th cavity couples to the $j$th atom with strength
$g_{j}$. The cavity leakage is assumed to be much larger than other decay and
dephasing rates in the array, which we thus neglect in the following. The
Floquet drive also dresses the array-cavity tunnelling, which hence acquires
the phase $\theta_{j}(t)$.

The atom-cavity Hamiltonian $\hat{H}(t)\equiv\sum_{m}\hat{H}_{m}e^{im\omega
t}$ gives rise to Floquet states $|\psi_{n}(t)\rangle=|u_{n}(t)\rangle
e^{-i\varepsilon_{n}t/\hbar}=|u_{n\mu}(t)\rangle
e^{-i\varepsilon_{n\mu}t/\hbar}$ with quasienergies
$\varepsilon_{n\mu}=\varepsilon_{n}+\mu\hbar\omega$ that are defined up to
integer multiples $\mu$ of $\hbar\omega$ and time-periodic Floquet modes
$|u_{n\mu}(t)\rangle=|u_{n}(t)\rangle e^{i\mu\omega
t}=|u_{n\mu}(t+T)\rangle\equiv\sum_{\bm{n},\bm{p},m}u_{n\mu}^{(\bm{n}\bm{p}m)}|\bm{n}\bm{p}\rangle
e^{im\omega t}$. Here $|\bm{n}\bm{p}\rangle$ denote Fock states with respect
to the occupation numbers $\bm{n}=(n_{1},\dots,n_{M})$ of system excitations
and $\bm{p}=(p_{1},\dots,p_{L})$ of cavity photons and $m$ is an integer
Fourier index. In Floquet space [36, 3], spanned by the Floquet-Fock states
$|\bm{n}\bm{p}m\rangle\\!\rangle$ (representing $|\bm{n}\bm{p}\rangle
e^{im\omega t}$ in the original space), the Floquet modes
$|u_{n\mu}\rangle\\!\rangle=\sum_{\bm{n},\bm{p},m}u_{n\mu}^{(\bm{n}\bm{p}m)}|\bm{n}\bm{p}m\rangle\\!\rangle$
are eigenstates with eigenvalue $\varepsilon_{n\mu}$ of the generalized
Hamiltonian $\hat{\mathcal{H}}$ [representing $\hat{H}(t)-i\hbar\partial_{t}$]
with matrix elements

$\langle\\!\langle\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}|\hat{\mathcal{H}}|\bm{n}\bm{p}m\rangle\\!\rangle=\langle\bm{n}^{\prime}\bm{p}^{\prime}|(\hat{H}_{m^{\prime}-m}+\delta_{m^{\prime}m}m\hbar\omega)|\bm{n}\bm{p}\rangle.$
(2)

The change of $m$ by $\Delta m=m^{\prime}-m\neq 0$, as it is described by the
Fourier components $\hat{H}_{\Delta m}$, corresponds to a resonant transition,
where the energy changes by $\Delta m$ driving quanta $\hbar\omega$. In turn,
processes without such driving-induced energy change are captured by the time-
averaged Hamiltonian $\hat{H}_{\Delta m=0}$, which contains
$\hat{H}^{\text{eff}}_{\text{S}}$.

For Floquet engineering, we aim at a regime, where $m$ remains a good quantum
number, such that the array coherent dynamics is indeed ruled by
$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}$. For reservoir engineering, we aim at a
regime, where the total excitation number $N$ remains a good quantum number,
since we strive for quantum simulation at conserved particle number. We, thus,
aim at a situation, where energy and excitation-number-changing processes can
be treated using degenerate perturbation theory, allowing us to block
diagonalize $\hat{\mathcal{H}}$ both with respect to $m$ and $N$. This is
challenging. While resonances, where $m$ or $N$ change individually are
routinely suppressed by making the associated energy costs, the frequency
$\hbar\omega$ or the atom-cavity detuning $\delta_{j}$, respectively, large,
it is now also required to avoid resonances where both $m$ and $N$ change.
Similar processes are sketched in Fig. 2(a)-(b), where an excitation can
escape [(a)] or be injected [(b)] in the array by exchanging energy
$m\hbar\omega$ with the Floquet drive. Although for a weak transverse single-
qubit drive as in Ref. [37] these are very high-order processes, they strongly
challenge usual dispersive-regime treatments here and cannot be neglected.
Nonetheless, as we will see, the desired regime can be found.

Figure 2: (a)-(b) Examples of unwanted excitation-photon conversion assisted
by the Floquet drive. In (a), an excitation tunnels from the array to a cavity
releasing an energy $m\hbar\omega$ to the drive, and leaks to the environment
being lost. In (b), a photon tunnels from the cavity to the array. (c)
Fundamental energy scales in the system.

A suitable description is achieved by performing a two-step block-
diagonalization of $\hat{\mathcal{H}}$ using van Vleck-type degenerate
perturbation theory in Floquet space [38, 39]. This approach captures
corrections arising from processes changing both $m$ and $N$ in second order
(except for second-order processes with respect to $m$ alone, which would
yield the first correction of $\hat{H}_{\mathrm{S}}^{\mathrm{eff}}$ in a high-
frequency expansion and are not relevant here). Such mixed processes would be
discarded when time averaging the total Hamiltonian before or after having
performed a dispersive-regime transformation, as it is done in Ref. [37]. This
makes our approach decisive for obtaining correct effective parameters both
for the coherent and incoherent dynamics, as confirmed by numerical
simulations. We arrive at the effective atom-cavity Hamiltonian [39]

$\hat{H}_{\mathrm{eff}}=\hat{H}_{\mathrm{S}}^{\mathrm{eff}}+\hat{H}_{\mathrm{C}}+\sum_{j=1}^{L}\big{[}\xi_{j}(\hat{n}_{j}-1)\hat{n}_{j}+\chi_{j}(\hat{n}_{j})\hat{c}_{j}^{\dagger}\hat{c}_{j}\big{]},$
(3)

with functions
$\chi_{j}(\hat{n}_{j})=\widetilde{\chi}_{j}(\hat{n}_{j})-\xi_{j}(\hat{n}_{j})$
and

$\displaystyle\widetilde{\chi}_{j}(\hat{n}_{j})=\sum_{m=-\infty}^{+\infty}\frac{g_{j}^{2}U\mathcal{J}_{m}^{2}(\lambda/\hbar\omega)\hat{n}_{j}}{[\delta_{j}+U(\hat{n}_{j}-1)-m\hbar\omega][\delta_{j}+U\hat{n}_{j}-m\hbar\omega]},$
$\displaystyle\xi_{j}(\hat{n}_{j})=\sum_{m=-\infty}^{+\infty}\frac{g_{j}^{2}\mathcal{J}_{m}^{2}(\lambda/\hbar\omega)}{\delta_{j}+U\hat{n}_{j}-m\hbar\omega}.$
(4)

The coupling to the cavities in Eq. (3) involves only operators preserving the
total number of atomic excitations. The prominent role played by the Floquet
drive is reflected in the effective tunnelling rate and atom-cavity coupling,
which explicitly depend on the driving frequency $\omega$ and amplitude
$\lambda$. Equation (3) can only be valid as long as resonances are avoided
that make the denominators of Eq. (Cavity-based reservoir engineering for
Floquet-engineered superconducting circuits) small. A sketch of the relation
between different energy scales in the system, that permit the approximations
adopted while being experimentally realistic, is shown in Fig. 2(c). The
starting point is the dispersive regime, where $|\delta_{j}|\gg g_{j}$, with a
strong nonlinearity $|U|\gg g_{j}$ that enhances the virtual coupling, see Eq.
(Cavity-based reservoir engineering for Floquet-engineered superconducting
circuits). Next $\hbar\omega,\lambda\gg J,g_{j}$ is chosen to enable Floquet
engineering of the tunnelling dynamics, while avoiding resonances.

Applying the perturbative treatment to the full master equation [40, 39], an
effective dissipator $\mathcal{D}_{\mathrm{eff}}(\hat{\rho})$ is obtained,
$\mathcal{D}_{\mathrm{eff}}(\hat{\rho})=\sum_{j=1}^{L}\kappa_{j}\sum_{m=-\infty}^{+\infty}\mathcal{D}[\hat{c}_{j,m}](\hat{\rho})$,
with

$\hat{c}_{j,m}=\hat{c}_{j}\delta_{m,0}+\frac{g_{j}\mathcal{J}_{m}(\lambda/\hbar\omega)}{\delta_{j}+U\hat{n}_{j}-m\hbar\omega}\hat{a}_{j}\
.$ (5)

Since the block diagonalization of $\hat{\mathcal{H}}$ mixes atom and cavity
degrees of freedom, $\mathcal{D}_{\mathrm{eff}}(\hat{\rho})$ does not involve
cavity decay only, but also a small perturbative term $\propto\hat{a}_{j}$,
describing excitation loss from the system. This is analogous to the undriven
case in the dispersive regime [32], except for additional driving-induced
decay channels with $m\neq 0$. The excitation loss is weak in the perturbative
regime assumed here, and can be counteracted by postselection [39]. In a
proof-of-principle implementation, where state tomography in the relevant
subspace is accessible, this can be done directly from the estimated density
matrix. In large systems, we consider observables that, while being key
signatures of the desired effects, carry also information about the total
excitation number allowing for postselection. These are site occupations and
excitation currents. The former are extracted directly by dispersive readout
of the atomic excitations; the latter can be detected as done in cold atom
experiments [41, 11], by biasing the on-site potentials $v_{j}(t)$ of pairs of
neighbouring sites (which is an excitation-conserving process) and measuring
the time evolution of site occupations.

The form of the array-cavity coupling in $\hat{H}_{\mathrm{eff}}$ puts us in a
position, where reservoir engineering for the atoms can be realized. The
underlying mechanism is that incoming pump photons that are detuned from the
cavity resonance, need to take (or give) energy from (to) the system, before
being emitted at the cavity frequency [42, 37, 35]. If the detuning matches an
energy gap in the atomic system, this induces a corresponding “dissipative”
transition. This type of quantum bath engineering has been implemented
experimentally with superconducting circuits, e.g., for a resonantly driven
two-level system and for a three-sites undriven Bose-Hubbard chain [37, 35].
In our setup, since the relevant atomic Hamiltonian is the Floquet-engineered
Hamiltonian $\hat{H}_{\mathrm{S}}^{\mathrm{eff}}$, we exploit this mechanism
to address transitions among effective eigenstates of the system with
artificial magnetic flux.

Figure 3: (a) Ladder geometry and coupling to the cavities. (b) Single-
excitation (left) and two-excitation (right) level structure and transition
energies addressed by the cavities (arrows). (c) Stroboscopic evolution of the
populations $p_{\eta}(t)=\bra{\eta}\hat{\rho}(t)\ket{\eta}$ of the eigenstates
$\ket{\eta}$ of $\hat{H}_{\text{S}}^{\text{eff}}$, for the full driven master
equation (solid) and the effective master equation (dashed). The grey shaded
area indicates the fraction of discarded states in postselection (black dashed
line for the effective model). The parameters are $\Phi=\pi/2$,
$\hbar\omega=20J$, $\bm{\delta}/\hbar\omega=(1.76,1.7)$, $U=8J$,
$\kappa_{1}=\kappa_{2}=0.1J$, $\bm{\mathcal{E}}/J=(1.2,0.5)$, $g_{1}=g_{2}=J$.
Inset: The size of the arrows and circles reproduces the current pattern and
excitation density in the final state. (d) As (c) but for two excitations in
the hard-core bosons subspace. The parameters are $\Phi=\pi/2$,
$\hbar\omega=20J$, $\bm{\delta}/\hbar\omega=(1.69,1.68,1.7)$, $U=8J$,
$\bm{\kappa}/J=(0.05,0.05,0.05)$, $\bm{\mathcal{E}}/J=(0.6,1.6,0.4)$,
$\bm{g}/J=(1,1,1)$.

To design a dissipative path, driving the atoms towards a target eigenstate of
$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}=\sum_{\eta}\varepsilon_{\eta}\ket{\eta}\\!\bra{\eta}$,
the cavity-pump detunings $d_{j}$ are set to match different effective energy
gaps $\varepsilon_{\eta}-\varepsilon_{\eta^{\prime}}$ in the array. Addressing
single gaps is possible provided such gaps are larger than the effective atom-
cavity coupling
$\chi_{j}^{(\eta\eta^{\prime})}\sqrt{\bar{n}_{\mathrm{ph},j}}$, where
$\bar{n}_{\mathrm{ph},j}$ is the mean photon number in the $j$th cavity and
$\chi_{j}^{(\eta\eta^{\prime})}=\bra{\eta}\chi_{j}(\hat{n}_{j})\ket{{\eta^{\prime}}}$
is the matrix element of the atomic coupling operator. The resonant transition
rate produced by the $j$th cavity is derived from the effective master
equation [39] as

$\Gamma_{\eta\to\eta^{\prime}}^{(j)}=4\bar{n}_{\mathrm{ph},j}|\chi_{j}^{(\eta\eta^{\prime})}|^{2}/\hbar^{2}\kappa_{j}\
,$ (6)

under the ‘bad-cavity’ condition that the photon leakage is strong compared to
the effective coupling,
$\hbar\kappa_{j}\gg\chi_{j}^{(\eta\eta^{\prime})}\sqrt{\bar{n}_{\mathrm{ph},j}}$.
This condition guarantees a clean monotonic exponential decay from
$\ket{\eta}$ to $\ket{\eta^{\prime}}$, but is not strictly essential: a
smaller $\kappa_{j}$ can still achieve the desired transition, though
producing more complex decay dynamics due the atom-cavity system reaching an
effective strong-coupling regime [37].

We test the proposed scheme for different systems with Floquet-engineered
artificial fluxes. We will compare our effective theory to a full simulation
of the time-dependent master equation for both system and cavities, using
typical circuit-QED parameters. Given the latter, the determination of the
parameters ensuring a successful protocol can be done systematically: while
$\hbar\omega$ and the shift $\Delta$ are chosen to be much larger than $J$ and
fine-tuned to avoid unwanted resonances, the pump frequencies on the cavities
realize the detuning conditions for reservoir engineering. The pump amplitudes
are adjusted accordingly to have an average number of one-to-two photons in
the cavities, which gave best results in simulations. Further details on
parameter choices are given in [39].

We start by considering one excitation in a ladder-type array [Fig. 3(a)]. The
ground state for non-zero magnetic flux $\Phi$ can exhibit chiral currents,
flowing unidirectionally along the edges of the ladder [11, 43, 44, 45, 46].
Although only one excitation is considered, including subspaces with several
excitations and their mutual interactions $U$ is essential, since it
influences virtual excitation-number fluctuations exploited for reservoir
engineering, see Eqs. (3) and (Cavity-based reservoir engineering for Floquet-
engineered superconducting circuits). In the two-plaquettes system of Fig.
3(a), the ground state is prepared using two cavities which realize a
“cascaded” cooling configuration along the effective spectrum as depicted in
Fig. 3(b) (left). The driving amplitude $\lambda$ is chosen to give the same
effective tunnelling rate $J_{\mathrm{eff}}$ along every lattice bond, while
the choice of driving phases and potential offsets yields Peierls phases
$r\Phi$ at the $r$th rung [39]. The simulated build up of population in the
ground state in time is shown in Fig. 3(c), leading to the final current and
density pattern depicted in the inset. The cooling is independent from the
initial state. While in Fig. 3(c) the initial state features the excitation
sitting in one atom, an equivalent population build-up is observed also for
mixed initial states, such as an infinite-temperature state in the single-
excitation manifold. This confirms that the cooling mechanism can reduce the
system entropy, in addition to lowering the energy. A cooling scheme in the
case of a three-plaquette ladder is reported in [39].

The protocol is effective also for multiple excitations. When a second
excitation is injected, due to the large on-site interaction $U$, singly- and
doubly-excited atoms define two essentially uncoupled subspaces [Fig.
3(b)-right], corresponding to hard-core bosons and a tightly bound pair,
respectively. The former subspace is most interesting (as, for instance, hard-
core bosons in finite two-dimensional lattices with homogeneous flux are
predicted to give rise to fractional-quantum-Hall(FQH)-type ground states [47,
48, 49]) and we can use three cavities to cool the system into the hard-core-
boson ground state. Starting with two excitations localized in different
sites, achievable by exciting two atoms to state
$a_{j^{\prime}}^{\dagger}a_{j}^{\dagger}\ket{0}$ with quick resonant pulses
before launching the Floquet drives, the population build-up in time is shown
in Fig. 3(d).

From Fig. 3(c) and (d), one can see that the effective master equation
reproduces the results of the full model very accurately, that the fraction of
discarded data in postselection remains small and does not render the
experiment inefficient, and that the target states can be prepared with rather
high fidelities. In the Supplemental Material [39], we discuss the geometry of
a diamond-chain of corner-sharing rhombic plaquettes with flux $\Phi=\pi$ and
the preparation of so called Aharonov Bohm cages— single-excitation states
that are localized in a subsystem via destructive interference [50, 33, 51].

Figure 4: (a) Bands in a 100 sites ladder and (b) interband cooling of one
excitation using two cavities. Parameters $\Phi=0.8\pi$, $\hbar\omega=20J$,
$\lambda\simeq 1.16\hbar\omega$, $d_{1}=d_{2}=J$,
$\delta_{1}=\delta_{2}=34.4J$, $g_{1}=g_{2}=J$,
$\mathcal{E}_{1}=\mathcal{E}_{2}=1.4J$. The shaded region indicates the
fraction of discarded records in postselection.

Until now we have addressed the preparation of (zero-temperature) ground
states in systems of moderate size. The essential ingredients demonstrated in
these examples also provide a perspective for potential applications in larger
systems. Namely, the preparation of a gapped ground state, like a Floquet-
engineered (topological) band insulator or a correlated fractional Chern
insulator, appear possible. Although for a large system, the number of
cavities cannot scale as quickly as the number of transitions, multiple
transitions can be controlled with a single cavity when they lie in an energy
window comparable with the cavity linewidth $\hbar\kappa$ [39]. This effect
can already be appreciated from the two-excitation example of Fig. 3(b)-(d),
where three cavities only are enough to work effectively in the two-excitation
subspace (21 states). Another issue with large systems are resonant
excitations. Namely the gapped state will most likely be embedded into a
continuum of excited states to which it couples resonantly. Although this
might limit the capability to resolve the exact ground state, playing a role
similar to an effective non-zero temperature, the preparation of low-energy
and low-entropy states can still be addressed.

We further exemplify the control of multiple transitions in a larger system,
by considering a 100-site square ladder with flux $\Phi=0.8\pi$ and tunnelling
along the rungs five times faster than along the legs, whose effective
spectrum features two well separated narrow bands [Fig. 4(a)]. One excitation
initialized at a site overlaps with all states in both bands, but can be
dissipatively pushed into the lower band using two cavities only, coupled to
the two ends of the upper leg. We demonstrate this successfully by performing
simulations with the effective master equation [Fig. 4(b)]. This example
further highlights that Floquet-dissipative schemes can also be used to
project a system into a (quasi)energetically well separated subspace, opening
potential applications for autonomous quantum error correction [52, 53] and
the preparation of gapped ground states (such as FQH states [47, 48, 54, 49]).
In view of the latter, particularly exciting is the combination of artificial
gauge fields [1, 4] or geometric frustration [55, 9], as it can be achieved
using Floquet engineering, with an interaction-induced hard-core constraint.
Such hard-core interactions can not only be used to mimic fermionic behavior
in 1D, but are also predicted to stabilize fractional-Chern-insulator states
in topologically non-trivial lattice systems as well as spin-liquid-like
states. Finally, the theory developed here can find straightforward
application also in different platforms where Floquet engineering has been
employed, such as quantum gas microscopes, whenever coupling to driven cavity-
like modes is possible.

This research was funded by the German Research Foundation (DFG) via the
Collaborative Research Center (SFB) 910, under project number 163436311, and
the Research Unit FOR 2414, under Project No. 277974659.

## References

  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Periodically Driven Quantum Systems: Effective Hamiltonians and Engineered Gauge Fields, Phys. Rev. X 4, 031027 (2014).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Adv. Phys. 64, 139 (2015).
  * Eckardt [2017] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017).
  * Aidelsburger  _et al._ [2018] M. Aidelsburger, S. Nascimbene, and N. Goldman, Artificial gauge fields in materials and engineered systems, C R Phys. 19, 394 (2018).
  * Oka and Kitamura [2019] T. Oka and S. Kitamura, Floquet Engineering of Quantum Materials, Annu. Rev. Condens. Matter Phys. 10, 387 (2019).
  * Rudner and Lindner [2020] M. S. Rudner and N. H. Lindner, Band structure engineering and non-equilibrium dynamics in Floquet topological insulators, Nat. Rev. Phys. 2, 229 (2020).
  * Lignier  _et al._ [2007] H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zenesini, O. Morsch, and E. Arimondo, Dynamical Control of Matter-Wave Tunneling in Periodic Potentials, Phys. Rev. Lett. 99, 220403 (2007).
  * Zenesini  _et al._ [2009] A. Zenesini, H. Lignier, D. Ciampini, O. Morsch, and E. Arimondo, Coherent Control of Dressed Matter Waves, Phys. Rev. Lett. 102, 100403 (2009).
  * Struck  _et al._ [2011] J. Struck, C. Ölschläger, R. L. Targat, P. Soltan-Panahi, A. Eckardt, M. Lewenstein, P. Windpassinger, and K. Sengstock, Quantum Simulation of Frustrated Classical Magnetism in Triangular Optical Lattices, Science 333, 996 (2011).
  * Struck  _et al._ [2012] J. Struck, C. Ölschläger, M. Weinberg, P. Hauke, J. Simonet, A. Eckardt, M. Lewenstein, K. Sengstock, and P. Windpassinger, Tunable gauge potential for neutral and spinless particles in driven optical lattices, Phys. Rev. Lett. 108, 225304 (2012).
  * Atala  _et al._ [2014] M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Observation of chiral currents with ultracold atoms in bosonic ladders, Nat. Phys. 10, 588 (2014).
  * Jotzu  _et al._ [2014] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Experimental realization of the topological Haldane model with ultracold fermions, Nature 515, 237 (2014).
  * Aidelsburger  _et al._ [2015] M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T. Barreiro, S. Nascimbène, N. R. Cooper, I. Bloch, and N. Goldman, Measuring the Chern number of Hofstadter bands with ultracold bosonic atoms, Nat. Phys. 11, 162 (2015).
  * Tai  _et al._ [2017] M. E. Tai, A. Lukin, M. Rispoli, R. S. A. Menke, D. Borgnia, P. M. Preiss, F. Grusdt, A. M. Kaufman, and M. Greiner, Microscopy of the interacting Harper–Hofstadter model in the two-body limit, Nature 546, 519 (2017).
  * Wintersperger  _et al._ [2020] K. Wintersperger, C. Braun, F. N. Ünal, A. Eckardt, M. Di Liberto, N. Goldman, I. Bloch, and M. Aidelsburger, Realization of an anomalous Floquet topological system with ultracold atoms, Nat. Phys. 16, 1058–1063 (2020).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Photonic Floquet topological insulators, Nature 496, 196–200 (2013).
  * Maczewsky  _et al._ [2017] L. J. Maczewsky, J. M. Zeuner, S. Nolte, and A. Szameit, Observation of photonic anomalous Floquet topological insulators, Nat. Commun. 8, 13756 (2017).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * Roushan  _et al._ [2017a] P. Roushan, C. Neill, A. Megrant, Y. Chen, R. Babbush, R. Barends, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, P. J. J. O’Malley, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. N. E. Kapit, and J. Martinis, Chiral ground-state currents of interacting photons in a synthetic magnetic field, Nat. Phys. 13, 146–151 (2017a).
  * Breuer and Holthaus [1991] H. Breuer and M. Holthaus, A semiclassical theory of quasienergies and Floquet wave functions, Ann. Phys. 211, 249 (1991).
  * Ketzmerick and Wustmann [2010] R. Ketzmerick and W. Wustmann, Statistical mechanics of Floquet systems with regular and chaotic states, Phys. Rev. E 82, 021114 (2010).
  * Shirai  _et al._ [2015] T. Shirai, T. Mori, and S. Miyashita, Condition for emergence of the Floquet-Gibbs state in periodically driven open systems, Phys. Rev. E 91, 030101 (2015).
  * Shirai  _et al._ [2016] T. Shirai, J. Thingna, T. Mori, S. Denisov, P. Hänggi, and S. Miyashita, Effective Floquet–Gibbs states for dissipative quantum systems, New J. Phys. 18, 053008 (2016).
  * Mori [2022] T. Mori, Floquet States in Open Quantum Systems, arXiv:2203.16358 (2022).
  * Seetharam  _et al._ [2015] K. I. Seetharam, C.-E. Bardyn, N. H. Lindner, M. S. Rudner, and G. Refael, Controlled Population of Floquet-Bloch States via Coupling to Bose and Fermi Baths, Phys. Rev. X 5, 041050 (2015).
  * Iadecola  _et al._ [2015] T. Iadecola, T. Neupert, and C. Chamon, Occupation of topological Floquet bands in open systems, Phys. Rev. B 91, 235133 (2015).
  * Zhang  _et al._ [2019] Y. Zhang, B. J. Lester, Y. Y. Gao, L. Jiang, R. J. Schoelkopf, and S. M. Girvin, Engineering bilinear mode coupling in circuit QED: Theory and experiment, Phys. Rev. A 99, 012314 (2019).
  * Zhang  _et al._ [2022] Y. Zhang, J. C. Curtis, C. S. Wang, R. J. Schoelkopf, and S. M. Girvin, Drive-induced nonlinearities of cavity modes coupled to a transmon ancilla, Phys. Rev. A 105, 022423 (2022).
  * Roushan  _et al._ [2017b] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis, Spectroscopic signatures of localization with interacting photons in superconducting qubits, Science 358, 1175 (2017b).
  * Ma  _et al._ [2019] R. Ma, B. Saxberg, C. Owens, N. Leung, Y. Lu, J. Simon, and D. I. Schuster, A dissipatively stabilized Mott insulator of photons, Nature 566, 51 (2019).
  * Yan  _et al._ [2019] Z. Yan, Y.-R. Zhang, M. Gong, Y. Wu, Y. Zheng, S. Li, C. Wang, F. Liang, J. Lin, Y. Xu, C. Guo, L. Sun, C.-Z. Peng, K. Xia, H. Deng, H. Rong, J. Q. You, F. Nori, H. Fan, X. Zhu, and J.-W. Pan, Strongly correlated quantum walks with a 12-qubit superconducting processor, Science 364, 753 (2019).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * Mukherjee  _et al._ [2018] S. Mukherjee, M. Di Liberto, P. Öhberg, R. R. Thomson, and N. Goldman, Experimental Observation of Aharonov-Bohm Cages in Photonic Lattices, Phys. Rev. Lett. 121, 075502 (2018).
  * Strand  _et al._ [2013] J. D. Strand, M. Ware, F. Beaudoin, T. A. Ohki, B. R. Johnson, A. Blais, and B. L. T. Plourde, First-order sideband transitions with flux-driven asymmetric transmon qubits, Phys. Rev. B 87, 220505 (2013).
  * Hacohen-Gourgy  _et al._ [2015] S. Hacohen-Gourgy, V. V. Ramasesh, C. De Grandi, I. Siddiqi, and S. M. Girvin, Cooling and Autonomous Feedback in a Bose-Hubbard Chain with Attractive Interactions, Phys. Rev. Lett. 115, 240501 (2015).
  * Sambe [1973] H. Sambe, Steady States and Quasienergies of a Quantum-Mechanical System in an Oscillating Field, Phys. Rev. A 7, 2203 (1973).
  * Murch  _et al._ [2012] K. W. Murch, U. Vool, D. Zhou, S. J. Weber, S. M. Girvin, and I. Siddiqi, Cavity-Assisted Quantum Bath Engineering, Phys. Rev. Lett. 109, 183602 (2012).
  * Eckardt and Anisimovas [2015] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a Floquet-space perspective, New. J. Phys. 17, 093039 (2015).
  * [39] See Supplemental Material which includes Refs. [56, 57, 58, 59, 60] and (i) the derivation of the effective Hamiltonian and master equation, with related driving schemes used; (ii) supplemental details on the derivation of the dissipative transition rates attained by reservoir engineering; (iii) details on the numerical simulations and postselection, and the preparation of localized states; (iv) supplemental details on the choice of parameters and an illustrative example of protocol design for Floquet-dissipative state preparation .
  * Schnell  _et al._ [2021] A. Schnell, S. Denisov, and A. Eckardt, High-frequency expansions for time-periodic Lindblad generators, Phys. Rev. B 104, 165414 (2021).
  * Keßler and Marquardt [2014] S. Keßler and F. Marquardt, Single-site-resolved measurement of the current statistics in optical lattices, Phys. Rev. A 89, 061601 (2014).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Introduction to quantum noise, measurement, and amplification, Rev. Mod. Phys. 82, 1155 (2010).
  * Hügel and Paredes [2014] D. Hügel and B. Paredes, Chiral ladders and the edges of quantum Hall insulators, Phys. Rev. A 89, 023619 (2014).
  * Piraud  _et al._ [2015] M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, S. Greschner, T. Vekua, and U. Schollwöck, Vortex and Meissner phases of strongly interacting bosons on a two-leg ladder, Phys. Rev. B 91, 140406 (2015).
  * Greschner  _et al._ [2015] S. Greschner, M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, U. Schollwöck, and T. Vekua, Spontaneous Increase of Magnetic Flux and Chiral-Current Reversal in Bosonic Ladders: Swimming against the Tide, Phys. Rev. Lett. 115, 190402 (2015).
  * Wang  _et al._ [2021] B. Wang, X.-Y. Dong, F. N. Ünal, and A. Eckardt, Robust and ultrafast state preparation by ramping artificial gauge potentials, New J. Phys. 23 (2021).
  * Sørensen  _et al._ [2005] A. S. Sørensen, E. Demler, and M. D. Lukin, Fractional Quantum Hall States of Atoms in Optical Lattices, Phys. Rev. Lett. 94, 086803 (2005).
  * Repellin  _et al._ [2020] C. Repellin, J. Léonard, and N. Goldman, Fractional Chern insulators of few bosons in a box: Hall plateaus from center-of-mass drifts and density profiles, Phys. Rev. A 102, 063316 (2020).
  * Wang  _et al._ [2022] B. Wang, X.-Y. Dong, and A. Eckardt, Measurable signatures of bosonic fractional Chern insulator states and their fractional excitations in a quantum-gas microscope, SciPost Phys. 12, 95 (2022).
  * Vidal  _et al._ [2000] J. Vidal, B. Douçot, R. Mosseri, and P. Butaud, Interaction Induced Delocalization for Two Particles in a Periodic Potential, Phys. Rev. Lett. 85, 3906 (2000).
  * Kremer  _et al._ [2020] M. Kremer, I. Petrides, E. Meyer, M. Heinrich, and O. Zilberberg, A square-root topological insulator with non-quantized indices realized with photonic Aharonov-Bohm cages, Nat. Commun. 11, 907 (2020).
  * Terhal [2015] B. M. Terhal, Quantum error correction for quantum memories, Rev. Mod. Phys. 87, 307 (2015).
  * Gertler  _et al._ [2021] J. M. Gertler, B. Baker, J. Li, S. Shiroi, J. Koch, and C. Wang, Protecting a bosonic qubit with autonomous quantum error correction, Nature 590, 243 (2021).
  * Liu  _et al._ [2021] Z. Liu, E. J. Bergholtz, and J. C. Budich, Dissipative preparation of fractional Chern insulators, Phys. Rev. Research 3, 043119 (2021).
  * Eckardt  _et al._ [2010] A. Eckardt, P. Hauke, P. Soltan-Panahi, C. Becker, K. Sengstock, and M. Lewenstein, Frustrated quantum antiferromagnetism with ultracold bosons in a triangular lattice, EPL 89, 10010 (2010).
  * Shavitt and Redmon [1980] I. Shavitt and L. T. Redmon, Quasidegenerate perturbation theories. A canonical van Vleck formalism and its relationship to other approaches, J. Chem. Phys. 73, 5711 (1980).
  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy,  _Quantum Optics, Chap. IV, Sec. 9.5_ (Cambridge University Press, Cambridge (UK), 1997).
  * Mølmer  _et al._ [1993] K. Mølmer, Y. Castin, and J. Dalibard, Monte carlo wave-function method in quantum optics, J. Opt. Soc. Am. B 10, 524 (1993).
  * Johansson  _et al._ [2012] J. Johansson, P. Nation, and F. Nori, QuTiP: An open-source Python framework for the dynamics of open quantum systems, Comp. Phys. Comm. 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Comp. Phys. Commun. 184, 1234 (2013).

Supplemental Material to

Cavity-based reservoir engineering for Floquet-engineered superconducting
circuits

Francesco Petiziol and André Eckardt

Technische Universität Berlin, Institut für Theoretische Physik,
Hardenbergstraße 36, Berlin 10623, Germany

##  I Effective master equation.

In this section, the effective Hamiltonian of Eq. (3) in the main text and the
effective master equation are derived.

###  I.1 Rotating frames.

We start by deriving the form of the Hamiltonian given in Eq. (1) in the main
text. The Hamiltonian of the atom array reads
$\hat{H}_{\mathrm{S}}(t)=\hat{H}_{\mathrm{op}}(t)+\hat{H}_{\mathrm{int}}+\hat{H}_{\mathrm{tun}}$
with

$\displaystyle\hat{H}_{\mathrm{op}}(t)=$
$\displaystyle\sum_{\ell=1}^{M}\hbar\dot{\theta}_{\ell}(t)\hat{n}_{\ell},\quad\hat{H}_{\mathrm{int}}=\frac{U}{2}\sum_{\ell=1}^{M}\hat{n}_{\ell}(\hat{n}_{\ell}-1),\quad\hat{H}_{\mathrm{tun}}=-J\sum_{\langle\ell,\ell^{\prime}\rangle}(\hat{a}_{\ell^{\prime}}^{\dagger}\hat{a}_{\ell}+\hat{a}_{\ell^{\prime}}\hat{a}_{\ell}^{\dagger}),$
(S7)

where $\langle\ell,\ell^{\prime}\rangle$ indicates ordered pairs of nearest
neighbours in the lattice. The periodically modulated site potential has
sinusoidal shape,
$\hbar\dot{\theta}_{\ell}(t)=\nu_{\ell}\hbar\omega+\lambda\sin(\omega
t-\varphi_{\ell}).$ The offsets $\nu_{\ell}\hbar\omega$, where
$\nu_{\ell}\in\\{0,1\\}$ in our applications, are needed for Floquet
engineering artificial magnetic fields. The Hamiltonian $H_{\mathrm{S}}(t)$ is
written in a frame rotating at frequency $\Delta/\hbar$, where $\Delta$ is a
static atomic potential equal for all atoms. The Hamiltonian for the $L$
cavities in the frame rotating at the same frequency $\Delta/\hbar$ reads

$\hat{H}_{\mathrm{C}}(t)=\sum_{\ell=1}^{L}\Big{(}\delta_{\ell}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}+\mathcal{E}_{\ell}\hat{c}_{\ell}^{\dagger}e^{-i\omega_{\ell}t}+\mathcal{E}^{*}_{\ell}\hat{c}_{\ell}e^{i\omega_{\ell}t}\Big{)},$
(S8)

such that $\delta_{\ell}$ represents the detuning between the cavity
frequencies and the static atomic potential $\Delta$. The quantity
$\omega_{\ell}$ is then the detuning between the frequency of the cavity pump
on the $\ell$th cavity from $\Delta$. The atom-cavity interaction is given by
the Hamiltonian
$\hat{H}_{\mathrm{SC}}=\sum_{\ell=1}^{L}g_{\ell}(\hat{a}_{\ell}\hat{c}_{\ell}^{\dagger}+\hat{a}_{\ell}^{\dagger}\hat{c}_{\ell})$,
characterized by the coupling strengths $g_{\ell}$. The dynamics of the
combined cavities-atoms system is then described by the Lindblad master
equation

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}[\hat{H}_{\mathrm{S}}(t)+\hat{H}_{\mathrm{C}}(t)+\hat{H}_{\mathrm{SC}},\hat{\rho}]+\sum_{\ell=1}^{L}\kappa_{\ell}\mathcal{D}[\hat{c}_{\ell}]\hat{\rho}\
,$ (S9)

where $\mathcal{D}[\hat{c}_{\ell}]\hat{\rho}$ is a Lindblad dissipator
describing cavity photon leakage,
$\mathcal{D}[\hat{c}_{\ell}]\hat{\rho}=\hat{c}_{\ell}\hat{\rho}\hat{c}_{\ell}^{\dagger}-\frac{1}{2}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{c}_{\ell}^{\dagger}\hat{c}_{\ell}.$
As discussed in Sec. III.1, this is the master equation used to simulate the
full driven dynamics of system and cavities. To arrive at Eq. (1), we
represent the master equation in the time-dependent frame comoving with the
Floquet drive and the potential offsets $\nu_{\ell}\hbar\omega$, i.e., such
that these terms are canceled from the comoving-frame Hamiltonian. The
comoving frame is defined by the unitary transformation

$\displaystyle R(t)$
$\displaystyle=\exp\Big{\\{}-i\sum_{\ell=1}^{M}\int_{0}^{t}\dot{\theta}_{\ell}(t)\hat{n}_{\ell}\Big{\\}},$
$\displaystyle=\exp\Big{\\{}-i\sum_{\ell=1}^{M}\Big{(}\nu_{\ell}\omega
t+\frac{\lambda}{\hbar\omega}\big{[}\cos(\varphi_{\ell})-\cos(\omega
t-\varphi_{\ell})\big{]}\Big{)}\hat{n}_{\ell}\Big{\\}}.$ (S10a)

Note that $R(0)=\mathbb{1}$ and that $R(t+T)=R(t)$ is time-periodic. Hence,
$R(t)$ maps to the identity at stroboscopic times $nT$, $n\in\mathbb{N}$, and
thus the rotating frame coincides with the laboratory frame at such times. The
comoving-frame Hamiltonian $R^{\dagger}(t)\hat{H}(t)R(t)-i\hbar
R^{\dagger}(t)\partial_{t}R(t)$ reads

$\hat{H}(t)=\hat{H}_{\mathrm{int}}+\hat{H}_{\mathrm{C}}(t)+\hat{H}_{\mathrm{tun}}(t)+\hat{H}_{\mathrm{SC}}(t),$
(S11)

with

$\qquad\hat{H}_{\mathrm{tun}}(t)=-J\sum_{\langle\ell,\ell^{\prime}\rangle}e^{i\theta_{\ell^{\prime}\ell}(t)}\hat{a}^{\dagger}_{\ell^{\prime}}\hat{a}_{\ell}+\mathrm{h.c.},\qquad\hat{H}_{\mathrm{SC}}(t)=\sum_{\ell=1}^{L}g_{\ell}e^{-i\theta_{\ell}(t)}\hat{a}_{\ell}\hat{c}^{\dagger}_{\ell}+\mathrm{h.c.},$
(S12)

The time-dependent Peierls phases $\theta_{\ell}(t)$ and
$\theta_{\ell\ell^{\prime}}(t)=\theta_{\ell}(t)-\theta_{\ell^{\prime}}(t)$ are
given by

$\displaystyle\theta_{\ell}(t)=\nu_{\ell}\omega
t-\frac{\lambda}{\hbar\omega}\cos(\omega t-\varphi_{\ell}),$ (S13a)
$\displaystyle\theta_{\ell\ell^{\prime}}(t)=\frac{2\lambda}{\hbar\omega}\sin\left(\frac{\varphi_{\ell}-\varphi_{\ell^{\prime}}}{2}\right)\sin\left(\omega
t-\frac{\varphi_{\ell}+\varphi_{\ell^{\prime}}}{2}\right)+\nu_{\ell\ell^{\prime}}\omega
t.$ (S13b)

with $\nu_{\ell\ell^{\prime}}=\nu_{\ell}-\nu_{\ell^{\prime}}$. Note that,
slightly differently from the main text, the frequency $\Delta/\hbar$ is not
included into $\theta_{\ell}(t)$ here, since Eq. (S7) was already written in a
frame rotating at $\Delta/\hbar$. The Hamiltonian $\hat{H}(t)$ of Eq. (S11)
thus corresponds to
$\hat{H}_{\mathrm{S}}(t)+\hat{H}_{\mathrm{C}}(t)+\hat{H}_{\mathrm{SC}}(t)$ in
the main text. In the transformation, we have absorbed time-independent phase
factors $e^{-i\lambda\cos(\varphi_{\ell})/\hbar\omega}$ into the operators,

$e^{i\lambda\cos(\varphi_{\ell})/\hbar\omega}\hat{a}_{\ell}^{\dagger}\to\hat{a}_{\ell}^{\dagger},$
(S14)

making a gauge transformation. In the following, it will be assumed that the
cavities are coupled to sites for which $\nu_{\ell}=0$, so we will set
$\nu_{\ell}=0$ in $\hat{H}_{\mathrm{SC}}(t)$. The dissipator remains
unmodified under the change of rotating frames.

Figure S1: Driving phases and resulting coupling patterns. The phases
$\varphi_{\ell}$ of the Floquet drives are indicated at the corresponding
lattice sites. Blue (red) links indicate a tunneling rate renormalized with
the Bessel function $\mathcal{J}_{0}$ ($\mathcal{J}_{1}$). Empty circles
indicate sites where there is a potential offset of $\hbar\omega$. Each arrow
indicates a phase factor $e^{i\Phi}$ representing the Peierls phases
engineered with the Floquet drives.

###  I.2  Artificial magnetic flux.

The Floquet engineering protocols considered in this work aim at realizing
effective magnetic fields, which enable the study of time-reversal symmetry
breaking and frustration with system excitations, which would not respond to a
physical magnetic field. We briefly review here, how the periodic drives yield
an artificial magnetic flux in our setup [3, 4]. Neglecting the cavities for
the moment, high-frequency Floquet engineering of the tunneling dynamics is
possible if the Floquet drives act on a timescale which is much faster than
the tunneling rate $J$, namely $\hbar\omega\gg J$ and
$\lambda\sim\hbar\omega$. The artificial flux is introduced by shaping with
the drives the time averaged tunneling Hamiltonian,
$\frac{1}{T}\int_{0}^{T}\\!dt\ \hat{H}_{\mathrm{tun}}(t).$ From Eq. (S12) and
using the generating function relation for the Bessel functions of the first
kind,
$e^{\frac{z}{2}(t-t^{-1})}=\sum_{m=-\infty}^{\infty}\mathcal{J}_{m}(z)t^{m}$,
for $z\in\mathbb{C},t\in\mathbb{C}\backslash\\{0\\}$, this time average gives

$\frac{1}{T}\int_{0}^{T}\\!dt\
\hat{H}_{\mathrm{tun}}(t)=-\sum_{\langle\ell,\ell^{\prime}\rangle}J^{\mathrm{eff}}_{\ell\ell^{\prime}}e^{i\theta^{\mathrm{eff}}_{\ell\ell^{\prime}}}\hat{a}_{\ell^{\prime}}^{\dagger}\hat{a}_{\ell}+\mathrm{h.c.}\
,$ (S15)

with effective tunneling rate as given in the main text,

$J^{\mathrm{eff}}_{\ell\ell^{\prime}}=J\mathcal{J}_{\nu_{\ell\ell^{\prime}}}\left(\frac{2\lambda}{\hbar\omega}\sin\big{[}(\varphi_{\ell^{\prime}}-\varphi_{\ell})/2\big{]}\right),$
(S16)

and effective phases
$\theta^{\mathrm{eff}}_{\ell\ell^{\prime}}=\nu_{\ell^{\prime}\ell}\frac{\varphi_{\ell^{\prime}}+\varphi_{\ell}}{2}.$
Beyond renormalizing the tunneling rate, the Floquet drives then permit one to
shape the complex tunneling phases
$\theta^{\mathrm{eff}}_{\ell\ell^{\prime}}$. In particular, note that when two
neighbouring sites $\ell,\ell^{\prime}$ have equal on-site potential offset,
$\nu_{\ell\ell^{\prime}}=0$, then the effective phase is zero, while it is
non-zero whenever $\nu_{\ell\ell^{\prime}}\neq 0$. In the examples used in the
numerical simulations and presented in the main text, the driving phases and
amplitude are chosen to give the same renormalized tunneling strength, at
chosen artificial flux $\Phi$. This is

$J^{\mathrm{eff}}_{\ell\ell^{\prime}}/J=\mathcal{J}_{0}\left(\frac{2\lambda}{\hbar\omega}\sin(\Phi/2)\right)=\mathcal{J}_{1}\left(\frac{2\lambda}{\hbar\omega}\sin(\Phi/2)\right),$
(S17)

with $\lambda$ given, by the condition that $\mathcal{J}_{0}$ and
$\mathcal{J}_{1}$ have the same value, $\lambda/\hbar\omega\simeq
0.72/\sin(\Phi/2)$. For flux $\Phi=\pi/2$, this yields the effective rate
$J_{\ell\ell^{\prime}}^{\mathrm{eff}}\simeq 0.55J$. The pattern of driving
phases and the resulting pattern of couplings and Peierls phases for the two-
plaquette ladder and rhombic lattices discussed in the article is shown in
Fig. S1.

###  I.3 “Floquet-dispersive” regime.

In the first perturbative step discussed in the text, the goal is to block-
diagonalize the system-cavity coupling in the Floquet Hamiltonian with respect
to the number $N$ of atomic excitations, thus extending the usual dispersive
regime of circuit QED [32] to the driven problem. This allows us to identify
regimes where system and cavities exchange energy without exchanging
excitations. To this end, the Hamiltonian of Eq. (S11) is split into an
unperturbed part $\hat{H}_{\mathrm{unp}}$, which commutes with
$\hat{N}=\sum_{j}\hat{n}_{j}$, and a perturbation $\hat{V}$ which we will
treat using van Vleck perturbation theory [56]. The unperturbed part is given
by the time independent terms,

$\hat{H}_{\mathrm{unp}}=\frac{U}{2}\sum_{j=1}^{M}\hat{n}_{j}(\hat{n}_{j}-1)+\sum_{j=1}^{L}\delta_{j}\hat{c}^{\dagger}_{j}\hat{c}_{j}\
,$ (S18)

where for the moment we do not include the cavity pump terms in
$\hat{H}_{\mathrm{C}}$ (see Section “Cavity pumps” below). The perturbation
reads

$\hat{V}(t)=\hat{H}_{\mathrm{tun}}(t)+\hat{H}_{\mathrm{SC}}(t)\ .$ (S19)

As discussed in the main text, the perturbation contains unwanted Floquet
$m$-photon processes that can lead to excitation loss or creation in the
system. The Floquet Hamiltonian $\hat{\mathcal{H}}$ in the extended space [36,
38] then reads

$\hat{\mathcal{H}}=\hat{\mathcal{H}}_{\mathrm{unp}}+\hat{\mathcal{V}}\ ,$
(S20)

with $\hat{\mathcal{H}}_{\mathrm{unp}}$ and $\hat{\mathcal{V}}$ defined as
follows. The unperturbed part $\hat{\mathcal{H}}_{\mathrm{unp}}$, representing
$\hat{H}_{\mathrm{unp}}$ of Eq. (S18) in the extended space, is defined as

$\hat{\mathcal{H}}_{\mathrm{unp}}=\sum_{\bm{np}}\sum_{m,m^{\prime}=-\infty}^{+\infty}|\bm{n}\bm{p}m\rangle\\!\rangle\varepsilon_{\bm{n}\bm{p}m}\langle\\!\langle\bm{n}\bm{p}m|\
,$ (S21)

where the unperturbed eigenstates $|\bm{n}\bm{p}m\rangle\\!\rangle$ are
labeled by: (i) the number $m$ of drive quanta $\hbar\omega$ in the Floquet
drive (ii) $\bm{n}=\\{n_{s}\\}_{s=1,\ldots,M}$, the set of excitation numbers
in the array (iii) $\bm{p}=\\{p_{s}\\}_{s=1,\ldots,L}$, the number of photons
in each cavity mode. In this section, we will use a calligraphic notation
$\hat{\mathcal{O}}$ to indicate the extended-space representation of an
operator $\hat{O}$. The unperturbed quasienergies explicitly read

$\displaystyle\varepsilon_{\bm{n}\bm{p}m}=$
$\displaystyle\bra{\bm{n}\bm{p}}\hat{H}_{\mathrm{unp}}\ket{\bm{n}\bm{p}}+m\hbar\omega\
,$ $\displaystyle=$
$\displaystyle\frac{U}{2}\sum_{\ell=1}^{M}n_{\ell}(n_{\ell}-1)+\sum_{\ell=1}^{L}\delta_{\ell}p_{\ell}+m\hbar\omega\
.$ (S22)

The perturbation $\hat{\mathcal{V}}$ reads

$\hat{\mathcal{V}}=\sum_{m,m^{\prime}}\ket{m}_{\\!F}\Big{(}\hat{H}_{\mathrm{tun},m-m^{\prime}}+\hat{H}_{\mathrm{SC},m-m^{\prime}}\Big{)}{}_{F}\\!\bra{m^{\prime}},$
(S23)

where we introduced the convenient short-hands

$\displaystyle\ket{m}_{\\!F}\equiv\sum_{\bm{p}\bm{n}}|\bm{p}\bm{n}m\rangle\\!\rangle\langle\bm{p}\bm{n}|,\qquad{}_{F}\\!\bra{m^{\prime}}\equiv\sum_{\bm{p}^{\prime}\bm{n}^{\prime}}|\bm{p}^{\prime}\bm{n}^{\prime}\rangle\langle\\!\langle\bm{p}^{\prime}\bm{n}^{\prime}m^{\prime}|.$
(S24)

From Eqs. (S12) and (S13), the Fourier components appearing in Eq. (S23) are

$\displaystyle\hat{H}_{\mathrm{tun},m}=$
$\displaystyle-\sum_{\langle\ell,\ell^{\prime}\rangle}\big{(}J_{\ell\ell^{\prime},m}\hat{a}_{\ell^{\prime}}^{\dagger}\hat{a}_{\ell}+J_{\ell^{\prime}\ell,m}\hat{a}_{\ell^{\prime}}\hat{a}_{\ell}^{\dagger}\big{)}\
,$ (S25a) $\displaystyle\hat{H}_{\mathrm{SC},m}=$
$\displaystyle\sum_{j=1}^{L}g_{j,m}[\hat{a}_{j}\hat{c}_{j}^{\dagger}+(-1)^{m}\hat{a}_{j}^{\dagger}\hat{c}_{j}]\
,$ (S25b)

where we have defined

$\displaystyle J_{\ell\ell^{\prime},m}=$ $\displaystyle
J\mathcal{J}_{m-\nu_{\ell^{\prime}\ell}}\left(\frac{2\lambda}{\hbar\omega}\sin\big{[}(\varphi_{\ell^{\prime}}-\varphi_{\ell})/2\big{]}\right)e^{-i(m-\nu_{\ell^{\prime}\ell})(\varphi_{\ell^{\prime}}+\varphi_{\ell})/2},$
(S26a) $\displaystyle g_{j,m}=$ $\displaystyle
g_{j}\mathcal{J}_{m}\Big{(}\frac{\lambda}{\hbar\omega}\Big{)}e^{im(\pi/2-\varphi_{j})}.$
(S26b)

We now decompose
$\hat{\mathcal{V}}=\hat{\mathcal{V}}_{\mathrm{D}}+\hat{\mathcal{V}}_{\mathrm{OD}}$
into its block-diagonal ($\hat{\mathcal{V}}_{\mathrm{D}}$) and block-off-
diagonal ($\hat{\mathcal{V}}_{\mathrm{OD}}$) part with respect to the total
excitation number $N$ in the array. The two components are then

$\displaystyle\hat{\mathcal{V}}_{\mathrm{D}}=\sum_{m,m^{\prime}=-\infty}^{+\infty}\ket{m}_{\\!F}\Big{(}\hat{H}_{\mathrm{tun},m-m^{\prime}}\Big{)}{}_{F}\\!\bra{m^{\prime}}\
,$ (S27a)
$\displaystyle\hat{\mathcal{V}}_{\mathrm{OD}}=\sum_{m,m^{\prime}=-\infty}^{+\infty}\ket{m}_{\\!F}\Big{(}\hat{H}_{\mathrm{SC},m-m^{\prime}}\Big{)}{}_{F}\\!\bra{m^{\prime}}\
.$ (S27b)

A transformation
$\hat{\mathcal{U}}=\exp(-\hat{\mathcal{G}})=\exp\left(-\sum_{k=1}^{+\infty}\hat{\mathcal{G}}^{(k)}\right)$
is searched which block-diagonalizes the array-cavity coupling in
$\hat{\mathcal{H}}$ with respect to $N$,

$\hat{\mathcal{U}}^{\dagger}\hat{\mathcal{H}}\hat{\mathcal{U}}=\hat{\mathcal{W}}\
,$ (S28)

with $\hat{\mathcal{W}}$ block-diagonal. Assuming that a perturbative
expansion is justified, namely if the unperturbed gaps are much larger than
the coupling to the perturbation, each generator $\hat{\mathcal{G}}^{(k)}$
will be of order $k$ in the perturbative parameter and one can impose the
block-diagonalization up to a desired order [38]. The relevant coupling in
$\hat{\mathcal{V}}_{\mathrm{D}}$ is $J$ while the one in
$\hat{\mathcal{V}}_{\mathrm{OD}}$ is $g_{j}$. Recall that these quantities,
and also the pump strengths $\mathcal{E}_{j}$, are assumed to be of similar
magnitude (see Fig. 2 in the main text). By expanding the expression (S28)
with the help of Baker-Campbell-Hausdorf formula, one finds that the first
order generator $\hat{\mathcal{G}}^{(1)}$ must satisfy

$\hat{\mathcal{V}}+[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{H}}_{\mathrm{unp}}]=0\
.$ (S29)

Choosing the generators to be fully off-diagonal with respect to the total
excitation number,

$\displaystyle\langle\\!\langle\bm{n}\bm{p}m|\hat{\mathcal{G}}^{(k)}|\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}\rangle\\!\rangle$
$\displaystyle=0,\quad\mbox{ for all }k\mbox{ if
}\sum_{j}n_{j}=\sum_{j}n_{j}^{\prime}\ ,$ (S30)

the first-order generator $\hat{\mathcal{G}}^{(1)}$ is then given by (for
$\sum_{j}n_{j}\neq\sum_{j}n_{j}^{\prime}$)

$\displaystyle\langle\\!\langle\bm{n}\bm{p}m|\hat{\mathcal{G}}^{(1)}|\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}\rangle\\!\rangle$
$\displaystyle=\frac{\langle\\!\langle\bm{n}\bm{p}m|\hat{\mathcal{V}}_{\mathrm{OD}}|\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}\rangle\\!\rangle}{\varepsilon_{\bm{n}\bm{p}m}-\varepsilon_{\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}}}.$
(S31)

The operator $\hat{\mathcal{V}}_{\mathrm{OD}}$, given by Eqs. (S27b) and
(S25b), is a sum of terms that exchange a single excitation between the system
and each cavity. Therefore, the relevant energy gaps in the denominator of Eq.
(S31), given occupation numbers $\\{\bm{n}\bm{p}\\}$, will only involve
$\\{\bm{n}^{\prime}\bm{p}^{\prime}\\}$ of the form
$\\{\bm{n}^{\prime}\bm{p}^{\prime}\\}=\\{n_{0},\ldots,n_{j}\pm
1\ldots,n_{M},p_{0},\dots,p_{j}\mp 1,\ldots p_{L}\\}$ and thus, using Eq.
(I.3),

$\varepsilon_{\bm{n}\bm{p}m}-\varepsilon_{\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}}=-\delta_{j}-Un_{j}+(m-m^{\prime})\hbar\omega,\quad\mbox{
or
}\quad\varepsilon_{\bm{n}\bm{p}m}-\varepsilon_{\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}}=\delta_{j}+U(n_{j}-1)+(m-m^{\prime})\hbar\omega,$
(S32)

The perturbative expansion is thus valid as long as these gaps remain much
larger than the couplings to the cavities $g_{j}$, for all $n_{j}$ and $m$,
$|\delta_{j}+Un_{j}+m\hbar\omega|\gg g_{j}.$ More precisely, this is true
under the above-mentioned assumption that the tunnelling rate $J$ in the
system and the strength $\mathcal{E}_{j}$ of the cavity pumps are of magnitude
comparable with $g_{j}$, such that products like
$Jg_{j}/(\varepsilon_{\bm{n}\bm{p}m}-\varepsilon_{\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}})$
or
$\mathcal{E}_{j}g_{j}/(\varepsilon_{\bm{n}\bm{p}m}-\varepsilon_{\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}})$
remain small and do not make the expansion diverge. For the denominators to be
sufficiently large, it is crucial that multiples $m\hbar\omega$ of Floquet
driving quanta must never be close to match the bare gaps, as discussed in the
main text and illustrated in Fig. 2. The second order generator
$\hat{\mathcal{G}}^{(2)}$ is found by imposing

$[\hat{\mathcal{G}}^{(2)},\hat{\mathcal{H}}_{\mathrm{unp}}]_{\mathrm{OD}}+[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{V}}]_{\mathrm{OD}}+\frac{1}{2}[\hat{\mathcal{G}}^{(1)},[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{H}}_{\mathrm{unp}}]]_{\mathrm{OD}}=0\
.$ (S33)

With the given choices of $\hat{\mathcal{G}}^{(1)}$ and
$\hat{\mathcal{G}}^{(2)}$, the block-diagonal operator $\hat{\mathcal{W}}$ up
to second order can be expressed as [38]

$\displaystyle\hat{\mathcal{W}}=$
$\displaystyle\hat{\mathcal{H}}_{\mathrm{unp}}+\hat{\mathcal{V}}_{\mathrm{D}}+\frac{1}{2}[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{V}}_{\mathrm{OD}}]_{\mathrm{D}}\
.$ (S34)

Therefore, what needs to be computed explicitly for determining
$\hat{\mathcal{W}}$ is the first-order generator $\hat{\mathcal{G}}^{(1)}$
only. Using the explicit expressions for $\hat{\mathcal{V}}$ from Eq. (S27)
and the unperturbed energies of Eq. (I.3) one finds

$\displaystyle\langle\\!\langle\bm{n}\bm{p}m|\hat{\mathcal{V}}_{\mathrm{OD}}|\bm{n}^{\prime}\bm{p}^{\prime}m^{\prime}\rangle\\!\rangle=$
$\displaystyle\sum_{j=1}^{L}g_{j,m-m^{\prime}}\Big{[}\sqrt{n_{j}+1}\sqrt{p_{j}}\delta_{n_{j},n_{j}^{\prime}-1}\delta_{p_{j},p_{j}^{\prime}+1}$
(S35)
$\displaystyle+(-1)^{m-m^{\prime}}\sqrt{n_{j}}\sqrt{p_{j}+1}\delta_{n_{j},n_{j}^{\prime}+1}\delta_{p_{j},p_{j}^{\prime}-1}\Big{]}\prod_{k\neq
j}\delta_{n_{j},n_{k}^{\prime}}\delta_{p_{j},p_{k}^{\prime}}.$ (S36)

Inserting into (S31), this gives then the “first-order” generator

$\hat{\mathcal{G}}^{(1)}=\sum_{j=1}^{L}\sum_{m,m^{\prime}=-\infty}^{+\infty}g_{j,m-m^{\prime}}\left[\hat{a}_{j}^{\dagger}\frac{(-1)^{m-m^{\prime}}}{\delta_{j}+U\hat{n}_{j}+\hbar\omega(m-m^{\prime})}\hat{c}_{j}-\frac{1}{\delta_{j}+U\hat{n}_{j}-\hbar\omega(m-m^{\prime})}\hat{a}_{j}\hat{c}_{j}^{\dagger}\right]\ket{m}_{\\!F}\\!\\!\bra{m^{\prime}}.$
(S37)

Returning to the non-extended space, the generator $\hat{G}^{(1)}(t)$ is time
periodic and reads

$\hat{G}^{(1)}(t)=\sum_{j=1}^{L}\sum_{m=-\infty}^{+\infty}g_{j,m}\left[\hat{a}_{j}^{\dagger}\frac{(-1)^{m}}{\delta_{j}+U\hat{n}_{j}+\hbar\omega
m}\hat{c}_{j}-\frac{1}{\delta_{j}+U\hat{n}_{j}-\hbar\omega
m}\hat{a}_{j}\hat{c}_{j}^{\dagger}\right]e^{im\omega t}.$ (S38)

The fact that $\hat{G}^{(1)}(t)$ is skew-hermitian,
$[\hat{G}^{(1)}(t)]^{\dagger}=-\hat{G}^{(1)}(t)$, can be seen by noting that
$g_{j,-m}^{*}=(-1)^{m}g_{j,m}$. For completing the computation of
$\hat{\mathcal{W}}$ of Eq. (S34), we further need the commutator
$[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{V}}_{\mathrm{OD}}]_{\mathrm{D}}$.
Introducing

$\displaystyle\hat{\zeta}_{j,m}(\hat{n}_{j})=\frac{1}{\delta_{j}+U\hat{n}_{j}+m\hbar\omega}\
,$ (S39)

the commutator gives

$\displaystyle\frac{1}{2}[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{V}}_{\mathrm{OD}}]_{\mathrm{D}}=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{L}\sum_{m,m^{\prime},\mu=-\infty}^{+\infty}C_{j,m-m^{\prime},\mu}\Big{\\{}(-1)^{m-m^{\prime}}\Big{(}\hat{\zeta}_{j,m-m^{\prime}+\mu}(\hat{n}_{j}-1)+\hat{\zeta}_{j,\mu}(\hat{n}_{j}-1)\Big{)}\hat{n}_{j}(1+\hat{c}_{j}^{\dagger}\hat{c}_{j})$
$\displaystyle-\Big{(}\hat{\zeta}_{j,m^{\prime}-m-\mu}(\hat{n}_{j})+\hat{\zeta}_{j,-\mu}(\hat{n}_{j})\Big{)}(\hat{n}_{j}+1)\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{\\}}\ket{m}_{\\!F}\\!\\!\bra{m^{\prime}}\
,$ (S40)

with
$C_{j,m,\mu}=g_{j}^{2}i^{m}e^{-im\varphi_{j}}\mathcal{J}_{m+\mu}\left(\frac{\lambda}{\hbar\omega}\right)\mathcal{J}_{\mu}\left(\frac{\lambda}{\hbar\omega}\right).$
Terms that would involve quadratic operators of the form
$\hat{c}_{j}^{\dagger}\hat{c}_{j}^{\dagger}$ and $\hat{c}_{j}\hat{c}_{j}$ are
part of the off-diagonal contribution
$[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{V}}_{\mathrm{OD}}]_{\mathrm{OD}}$
which is canceled by the second-order generator $\hat{\mathcal{G}}^{(2)}$ in
Eq. (S33). Terms involving two different cavities, e.g. involving operators of
the form
$\hat{a}_{j}^{\dagger}\hat{a}_{j^{\prime}}\hat{c}_{j}\hat{c}_{j^{\prime}}^{\dagger}$
with $j\neq j^{\prime}$, are included in
$[\hat{\mathcal{G}}^{(1)},\hat{\mathcal{V}}_{\mathrm{OD}}]_{\mathrm{D}}$ but
they ultimately vanish. The final effective Floquet Hamiltonian in the
extended space can then be reconstructed inserting Eqs. (S21), (S27a) and
(I.3) into Eq. (S34). Explicitly, it reads

$\displaystyle\hat{\mathcal{W}}=$
$\displaystyle\sum_{m,m^{\prime}=-\infty}^{+\infty}\Big{(}(\hat{H}_{\mathrm{unp}}+m\hbar\omega)\delta_{m,m^{\prime}}+\hat{H}_{\mathrm{tun},m-m^{\prime}}\Big{)}\ket{m}_{\\!F}\\!\\!\bra{m^{\prime}}$
$\displaystyle+\frac{1}{2}\sum_{m,m^{\prime}=-\infty}^{+\infty}\sum_{j=1}^{L}\sum_{\mu=-\infty}^{+\infty}C_{j,m-m^{\prime},\mu}\Big{\\{}(-1)^{m-m^{\prime}}\Big{(}\hat{\zeta}_{m-m^{\prime}+\mu}(\hat{n}_{j}-1)+\hat{\zeta}_{j,\mu}(\hat{n}_{j}-1)\Big{)}\hat{n}_{j}$
$\displaystyle+\Big{[}(-1)^{m-m^{\prime}}\Big{(}\hat{\zeta}_{m-m^{\prime}+\mu}(\hat{n}_{j}-1)+\hat{\zeta}_{j,\mu}(\hat{n}_{j}-1)\Big{)}\hat{n}_{j}-\Big{(}\hat{\zeta}_{j,m^{\prime}-m-\mu}(\hat{n}_{j})+\hat{\zeta}_{j,-\mu}(\hat{n}_{j})\Big{)}(\hat{n}_{j}+1)\Big{]}\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{\\}}\ket{m}_{\\!F}\\!\\!\bra{m^{\prime}}\
.$ (S41)

Rearranging terms and returning to the non-extended space, this leads to the
time-dependent Hamiltonian

$\displaystyle\hat{H}^{\prime}(t)=\hat{H}_{\mathrm{op}}(t)+\hat{H}_{\mathrm{int}}+\hat{H}_{C}+\hat{H}_{\mathrm{tun}}(t)+\hat{H}_{\mathrm{SC}}^{\prime}(t),$
(S42)

where

$\displaystyle\hat{H}_{\mathrm{op}}(t)=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{L}\sum_{m,\mu=-\infty}^{+\infty}C_{j,m,\mu}(-1)^{m}\Big{(}\hat{\zeta}_{j,m+\mu}(\hat{n}_{j}-1)+\hat{\zeta}_{j,\mu}(\hat{n}_{j}-1)\Big{)}\hat{n}_{j}e^{im\omega
t},$ (S43) $\displaystyle\hat{H}_{\mathrm{SC}}^{\prime}(t)=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{L}\sum_{m,\mu=-\infty}^{+\infty}C_{j,m,\mu}\Big{\\{}\Big{[}(-1)^{m}\Big{(}\hat{\zeta}_{j,m+\mu}(\hat{n}_{j}-1)+\hat{\zeta}_{j,\mu}(\hat{n}_{j}-1)\Big{)}$
$\displaystyle-\Big{(}\hat{\zeta}_{j,-m-\mu}(\hat{n}_{j})+\hat{\zeta}_{j,-\mu}(\hat{n}_{j})\Big{)}\Big{]}\hat{n}_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}-\Big{(}\hat{\zeta}_{j,-m-\mu}(\hat{n}_{j})+\zeta_{j,-\mu}(\hat{n}_{j})\Big{)}\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{\\}}e^{im\omega
t}.$ (S44)

Additional terms describing the transformed cavity pumps are discussed below.

###  I.4 High-frequency expansion

Under the assumptions that justify the Floquet-dispersive regime described in
the previous section, the dynamics can be studied separately within subspaces
with different total number of excitations $N$ in the array. We now conclude
the computation of a time-independent effective Hamiltonian by performing a
high-frequency expansion for $\hbar\omega\gg J$. In the extended-space
formalism, this corresponds to further block-diagonalizing the quasienergy
operator $\hat{\mathcal{W}}$, now with respect to the Floquet index. As
discussed, e.g., in Ref. [38], the result of this block-diagonalization is
that the Hamiltonian in the leading order is, in the non-extended space, the
time averaged Hamiltonian. From the extended space representation, this can be
extracted from the Floquet Hamiltonian simply by considering any diagonal
block with $m=m^{\prime}$ (modulo $\hbar\omega$),
$\hat{H}_{\mathrm{eff}}=\bra{m}\hat{\mathcal{W}}\ket{m}_{\\!F}.$ One then
obtains

$\displaystyle\hat{H}_{\mathrm{eff}}=$
$\displaystyle\hat{H}_{\mathrm{unp}}+\hat{H}_{\mathrm{tun},0}+\sum_{j=1}^{L}\sum_{m=-\infty}^{+\infty}\frac{g_{j}^{2}\mathcal{J}_{m}^{2}\left(\lambda/\hbar\omega\right)}{\delta_{j}+U(\hat{n}_{j}-1)+\hbar\omega
m}\hat{n}_{j}+$
$\displaystyle\sum_{j=1}^{L}\sum_{m=-\infty}^{+\infty}g_{j}^{2}\mathcal{J}_{m}^{2}\left(\frac{\lambda}{\hbar\omega}\right)\Big{[}\frac{1}{\delta_{j}+U(\hat{n}_{j}-1)+\hbar\omega
m}\hat{n}_{j}-\frac{1}{\delta_{j}+U\hat{n}_{j}-\hbar\omega
m}(\hat{n}_{j}+1)\Big{]}\hat{c}_{j}^{\dagger}\hat{c}_{j},$

which can be rewritten in the form of Eq. (3) of the main text,

$\hat{H}_{\mathrm{eff}}=\hat{H}_{\mathrm{unp}}-\sum_{\langle
j,j^{\prime}\rangle}\Big{(}J^{\mathrm{eff}}_{jj^{\prime}}e^{i\theta^{\mathrm{eff}}_{jj^{\prime}}}\hat{a}_{j^{\prime}}^{\dagger}\hat{a}_{j}+\mathrm{h.c.}\Big{)}+\sum_{j=1}^{L}\Big{(}\xi_{j}(\hat{n}_{j}-1)\hat{n}_{j}+\chi_{j}(\hat{n}_{j})\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{)},$
(S45)

with the functions $\xi_{j}(\hat{n}_{j})$ and $\chi_{j}(\hat{n}_{j})$ given in
Eq. (4), namely

$\displaystyle\xi_{j}(\hat{n}_{j})=\sum_{m=-\infty}^{+\infty}\frac{g_{j}^{2}\mathcal{J}_{m}^{2}\left(\lambda/\hbar\omega\right)}{\delta_{j}+U(\hat{n}_{j}-1)+\hbar\omega
m},$ (S46)
$\displaystyle\chi_{j}(\hat{n}_{j})=\sum_{m=-\infty}^{+\infty}g_{j}^{2}\mathcal{J}_{m}^{2}\left(\frac{\lambda}{\hbar\omega}\right)\Big{[}\frac{U}{(\delta_{j}+U(\hat{n}_{j}-1)-\hbar\omega
m)(\delta_{j}+U\hat{n}_{j}-\hbar\omega
m)}\hat{n}_{j}-\frac{1}{\delta_{j}+U\hat{n}_{j}-\hbar\omega m}\Big{]}.$ (S47)

In deriving the last expression, we have used the fact that
$\mathcal{J}_{m}^{2}\left(\lambda/\hbar\omega\right)=\mathcal{J}_{-m}^{2}\left(\lambda/\hbar\omega\right)$
to change $m\to-m$ in one of the terms. The Hamiltonian
$\hat{H}_{\mathrm{eff}}$ is then the result of our perturbative treatment for
the coherent dynamics.

For potential interest in the study of the single excitation dynamics, we
restrict the effective Hamiltonian of Eq. (S45) to the single-excitation
subspace, spanned by the states $\hat{a}_{j}^{\dagger}\ket{0}$. For the terms
$\xi_{j}(\hat{n}_{j}-1)\hat{n}_{j}$ and $\chi_{j}(\hat{n}_{j})$, this simply
amounts to substituting $\hat{n}_{j}\to 1$ in the denominators. The term
$\xi_{j}(\hat{n}_{j})$ gives instead

$\displaystyle\sum_{\ell,\ell^{\prime}}\sum_{j=1}^{L}\hat{a}^{\dagger}_{\ell}\ket{0}\bra{0}\hat{a}_{\ell}\frac{1}{\delta_{j}+U\hat{n}_{j}-\hbar\omega
m}\hat{a}^{\dagger}_{\ell^{\prime}}\ket{0}\bra{0}\hat{a}_{\ell^{\prime}}=$
$\displaystyle\sum_{j=1}^{L}\frac{1}{\delta_{j}-\hbar\omega
m}\mathbb{1}_{\mathrm{1p}}+\sum_{j=1}^{L}\frac{-U}{(\delta_{j}+U-\hbar\omega
m)(\delta_{j}-\hbar\omega m)}\hat{n}_{j}\ ,$ (S48)

where
$\mathbb{1}_{\mathrm{1p}}=\sum_{\ell}\hat{a}_{\ell}^{\dagger}\ket{0}\bra{0}\hat{a}_{\ell}$
is the identity operator in the single excitation subspace. The single-
excitation effective Hamiltonian $\hat{H}_{\mathrm{eff}}^{(1)}$ is thus

$\hat{H}_{\mathrm{eff}}^{(1)}=\hat{H}_{\mathrm{unp}}+\hat{H}_{\mathrm{tun},0}+\sum_{j=1}^{L}\Big{[}\xi_{j}^{(1)}(\hat{n}_{j}-\hat{c}_{j}^{\dagger}\hat{c}_{j})+\chi_{j}^{(1)}\hat{n}_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{]}\
,$ (S49)

with

$\xi_{j}^{(1)}=\sum_{m=-\infty}^{+\infty}g_{j}^{2}\frac{\mathcal{J}_{m}^{2}\left(\lambda/\hbar\omega\right)}{\delta_{j}-\hbar\omega
m},\qquad\chi_{j}^{(1)}=\sum_{m=-\infty}^{+\infty}g_{j}^{2}\frac{2U\mathcal{J}_{m}^{2}(\lambda/\hbar\omega)}{(\delta_{j}-\hbar\omega
m)(\delta_{j}+U-\hbar\omega m)}\ .$ (S50)

###  I.5 Cavity pumps

We did not yet consider the impact of the block-diagonalization with respect
to $N$ on the cavity pump term, such that block-diagonalizing the array-cavity
coupling with respect to $N$ was equivalent to block-diagonalizing the whole
quasienergy operator. In order to introduce the cavity pumps at the level of
the effective Hamiltonian, one needs to compute how these terms are modified
by the perturbative transformation. In summary, the transformation introduces
off-diagonal driving terms for the artificial atoms that couple to a cavity.
These terms are largely off-resonant and can thus be neglected in our
applications.

Focusing on each cavity pump at a time, it is convenient to work in a frame
rotating at the corresponding pump frequency $\omega_{j}$ for the whole system
(similarly to how one would approach a multilevel system subject to multiple
pulses, where the effect of each pulse can be better understood in a frame
rotating at the corresponding drive frequency). In this frame, the pump term
becomes time independent,

$\hat{H}_{j,\mathrm{pump}}=\mathcal{E}_{j}\hat{c}^{\dagger}_{j}+\mathcal{E}_{j}^{*}\hat{c}_{j}\
.$ (S51)

and the atoms and cavities energies acquire a potential shift
$-\hbar\omega_{j}$. We can then perform the Floquet-dispersive regime
transformation and high-frequency expansion for the $j$th cavity in this
frame. The resulting effective Hamiltonian terms are the same as in Eq. (S45),
with the additional potential shift $-\hbar\omega_{j}$ for atoms and cavities.

The Floquet-dispersive regime transformation $U(t)$ then acts on the pump
according to

$\displaystyle\hat{U}^{\dagger}(t)\hat{H}_{j,\mathrm{pump}}\hat{U}(t)=\sum_{j=1}^{L}[\mathcal{E}_{j}\hat{U}^{\dagger}(t)\hat{c}^{\dagger}_{j}\hat{U}(t)+\mathcal{E}_{j}^{*}U^{\dagger}(t)\hat{c}_{j}\hat{U}(t)]\
.$ (S52)

The transformed cavity operators can be computed perturbatively through
Hausdorf expansion, finding

$\displaystyle\hat{c}_{j}(t)\equiv
U^{\dagger}(t)\hat{c}_{j}U(t)\simeq\hat{c}_{j}+\sum_{m=-\infty}^{+\infty}\frac{g_{j,m}}{\delta_{j}+U\hat{n}_{j}-m\hbar\omega}\hat{a}_{j}e^{im\omega
t},$ (S53a) $\displaystyle\hat{c}_{j}^{\dagger}(t)\equiv
U^{\dagger}(t)\hat{c}^{\dagger}_{j}U(t)\simeq\hat{c}_{j}^{\dagger}+\sum_{m=-\infty}^{+\infty}\hat{a}^{\dagger}_{j}\frac{g^{*}_{j,-m}}{\delta_{j}+U\hat{n}_{j}+m\hbar\omega}e^{im\omega
t},$ (S53b)

with $g_{j,m}$ given in Eq. (S26). Note that we have kept only the first
commutator term in the Hausdorf expansion since $\mathcal{E}_{j}\sim g_{j}$ in
our model, and thus the single-commutator term is already of second order in
the perturbative expansion when inserted into Eq. (S52). The leading order in
the high-frequency expansion is then the time average of these operators,
namely

$\displaystyle\hat{c}_{j,0}=\hat{c}_{j}+\frac{g_{j}\mathcal{J}_{0}(\lambda/\hbar\omega)}{\delta_{j}+U\hat{n}_{j}}\hat{a}_{j},\qquad(\hat{c}_{j}^{\dagger})_{0}=\hat{c}_{j}^{\dagger}+\hat{a}^{\dagger}_{j}\frac{g_{j}\mathcal{J}_{0}(\lambda/\hbar\omega)}{\delta_{j}+U\hat{n}_{j}}.$
(S54)

The corrections to $\hat{c}_{j}$ and $\hat{c}_{j}^{\dagger}$ thus describe an
extra drive of the artificial atoms connecting manifolds with different
excitation number. The whole transformed pump Hamiltonian then becomes

$\hat{H}_{j,\mathrm{pump}}^{\prime}=\hat{H}_{j,\mathrm{pump}}+\sum_{j=1}^{L}\Big{(}\frac{\widetilde{\mathcal{E}}_{j}}{\delta_{j}+U\hat{n}_{j}}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\frac{\widetilde{\mathcal{E}}_{j}^{*}}{\delta_{j}+U\hat{n}_{j}}\Big{)},$
(S55)

with
$\widetilde{\mathcal{E}}_{j}=g_{j}\mathcal{E}_{j}\mathcal{J}_{0}(\lambda/\hbar\omega)$.
In the undriven system, this effective driving Hamiltonian for the artificial
atoms is used for exciting the array (implementing single-qubit rotations in
quantum information processing applications) if the pump frequency is resonant
with the atomic transitions. Here, since the cavity pumps are far detuned and
under the assumptions of validity of the Floquet-dispersive regime and high-
frequency expansion, these extra driving terms on the array are far off-
resonant and can be neglected. Returning to the frame rotating at
$\Delta/\hbar$, the cavity pump term in conclusion remains unaltered,

$\hat{H}_{j,\mathrm{pump}}^{\mathrm{eff}}(t)=\hat{H}_{j,\mathrm{pump}}(t)=\mathcal{E}_{j}\hat{c}^{\dagger}_{j}e^{-i\omega_{j}t}+\mathcal{E}_{j}^{*}\hat{c}_{j}e^{i\omega_{j}t}\
.$ (S56)

###  I.6 Effective dissipator

The perturbation theory discussed for the Hamiltonian can be reformulated at
the superoperator level to compute an effective dissipator via the same
procedure. The master equation (S9) in the frame where the Floquet drive has
been integrated out can be rewritten as
$d\hat{\rho}/dt=\mathcal{L}(t)\hat{\rho}$, with Lindbladian superoperator

$\mathcal{L}(t)=-\frac{i}{\hbar}[\hat{H}(t),\cdot]+\sum_{j}\kappa_{j}\mathcal{D}[\hat{c}_{j}](\cdot),$
(S57)

and Hamiltonian $\hat{H}(t)$ given in Eq. (S11). The unitary Floquet-
dispersive regime transformation $\hat{U}(t)$ that block-diagonalizes the
Hamiltonian with respect to the total number of atom excitations transforms
the Lindbladian $\mathcal{L}(t)$ into

$\mathcal{L}^{\prime}(t)=\hat{U}^{\dagger}(t)\mathcal{L}(t)\hat{U}(t)=\frac{i}{\hbar}[\hat{H}^{\prime}(t),\cdot]+\sum_{j}\kappa_{j}\mathcal{D}[\hat{U}(t)^{\dagger}\hat{c}_{j}\hat{U}(t)](\cdot).$
(S58)

with
$\hat{H}^{\prime}(t)=\hat{U}^{\dagger}(t)\hat{H}(t)\hat{U}(t)-i\hbar\hat{U}^{\dagger}(t)\partial_{t}\hat{U}(t)$.
To second order, the transformed Hamiltonian $\hat{H}^{\prime}(t)$ is the
Floquet-dispersive-regime Hamiltonian given explicitly in Eq. (S42). Focusing
on the dissipative part, one can expand the operators
$\hat{c}_{j}(t)=\hat{U}^{\dagger}(t)\hat{c}_{j}\hat{U}(t)$ as a Fourier
series,

$\displaystyle\hat{c}_{j}(t)=\sum_{m=-\infty}^{+\infty}\hat{c}_{j,m}e^{im\omega
t},$ (S59)
$\displaystyle\hat{c}^{\dagger}_{j}(t)=\sum_{m=-\infty}^{+\infty}(\hat{c}_{j,m})^{\dagger}e^{-im\omega
t}.$ (S60)

The leading-order expressions for $\hat{c}_{j,m}$ are found straightforwardly
from Eq. (S53), leading to the expression given in the main text, Eq. (5). The
dissipators then read

$\mathcal{D}[\hat{c}_{j}(t)](\cdot)=\sum_{m,m^{\prime}=-\infty}^{+\infty}\Big{\\{}\hat{c}_{j,m}(\cdot)[\hat{c}_{j,m^{\prime}}]^{\dagger}-\frac{1}{2}[\hat{c}_{j,m^{\prime}}]^{\dagger}\hat{c}_{j,m}(\cdot)-\frac{1}{2}(\cdot)[\hat{c}_{j,m^{\prime}}]^{\dagger}\hat{c}_{j,m}\Big{\\}}e^{i(m-m^{\prime})\omega
t}.$ (S61)

Redefining $m^{\prime}=m-\mu$, this becomes

$\mathcal{D}[\hat{c}_{j}(t)](\cdot)=\sum_{m,\mu=-\infty}^{+\infty}\Big{\\{}\hat{c}_{j,m}(\cdot)[\hat{c}_{j,m-\mu}]^{\dagger}-\frac{1}{2}[\hat{c}_{j,m-\mu}]^{\dagger}\hat{c}_{j,m}(\cdot)-\frac{1}{2}(\cdot)[\hat{c}_{j,m-\mu}]^{\dagger}\hat{c}_{j,m}\Big{\\}}e^{i\mu\omega
t},$ (S62)

such that the $\mu$-th Fourier component of the dissipator can be immediately
identified. Finally, the time average of the transformed Lindbladian following
from the high-frequency expansion selects the zero-frequency component only,
which is

$\mathcal{D}_{\mathrm{eff}}(\cdot)=\sum_{j=1}^{L}\kappa_{j}\sum_{m=-\infty}^{+\infty}\Big{\\{}\hat{c}_{j,m}(\cdot)[\hat{c}_{j,m}]^{\dagger}-\frac{1}{2}[\hat{c}_{j,m}]^{\dagger}\hat{c}_{j,m}(\cdot)-\frac{1}{2}(\cdot)[\hat{c}_{j,m}]^{\dagger}\hat{c}_{j,m}\Big{\\}}.$
(S63)

The time averaged effective Lindbladian then reads

$\mathcal{L}_{\mathrm{eff}}=-\frac{i}{\hbar}[\hat{H}_{\mathrm{eff}},\cdot]+\mathcal{D}_{\mathrm{eff}}(\cdot),$
(S64)

with the effective Hamiltonian of Eq. (S45). Simulations of the effective
dynamics are performed using this master equation, including the cavity pumps.

##  II Reservoir engineering for the stroboscopic dynamics

We now discuss how reservoir engineering is realized at the level of the
effective stroboscopic Hamiltonian of the array shaped through the Floquet
drives. Our starting point is the effective master equation (S64) derived in
the previous Section, neglecting the effective excitation decay from the array
(see section “Postselection”). It reads

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}[\hat{H}_{\mathrm{eff}},\hat{\rho}]+\sum_{j}\kappa_{j}\mathcal{D}[\hat{c}_{j}](\hat{\rho}).$
(S65)

The theory developed here follows similar ones done in other contexts [37,
35], with the difference that here multiple cavities are considered and that
the relevant system Hamiltonian is the Floquet-engineered Hamiltonian with
artificial magnetic flux. The effective Hamiltonian can be rewritten in the
form

$\hat{H}_{\mathrm{eff}}=\hat{H}_{\mathrm{S}}^{\mathrm{eff}}+\sum_{j=1}^{L}\Big{[}\chi_{j}(\hat{n}_{j})\hat{c}_{j}^{\dagger}\hat{c}_{j}+\delta_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\mathcal{E}_{j}\hat{c}_{j}^{\dagger}e^{-i\omega_{j}t}+\mathcal{E}_{j}^{*}\hat{c}_{j}e^{i\omega_{j}t}\Big{]},$
(S66)

where the effective atomic Hamiltonian is

$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}=\sum_{j=1}^{M}\Big{\\{}\xi_{j}(\hat{n}_{j}-1)\hat{n}_{j}+\frac{U}{2}\hat{n}_{j}(\hat{n}_{j}-1)\Big{\\}}-\sum_{\langle
j,j^{\prime}\rangle}J^{\mathrm{eff}}_{jj^{\prime}}(e^{i\theta^{\mathrm{eff}}_{jj^{\prime}}}\hat{a}_{j}\hat{a}_{j^{\prime}}^{\dagger}+\mathrm{h.c.}).$
(S67)

Note that we have included the effective term
$\sum_{j}\xi_{j}(\hat{n}_{j}-1)\hat{n}_{j}$ into the effective atomic
Hamiltonian, slightly differently from Eq. (3) of the main text. Since the
effective atom-cavity coupling is diagonal with respect to the cavity Fock
basis, we can make a gauge transformation for the cavities which eliminates
the time-dependence given by the cavity pump terms (corresponding to a frame
rotating at $\omega_{j}$ for the $j$th cavity) without altering the atomic
part. The resulting $\hat{H}_{\mathrm{eff}}$ reads

$\hat{H}_{\mathrm{eff}}=\hat{H}_{\mathrm{S}}^{\mathrm{eff}}+\sum_{j=1}^{L}\Big{[}\chi_{j}(\hat{n}_{j})\hat{c}_{j}^{\dagger}\hat{c}_{j}+d_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\mathcal{E}_{j}\hat{c}_{j}^{\dagger}+\mathcal{E}_{j}^{*}\hat{c}_{j}\Big{]},$
(S68)

where we have introduced the cavity-pump detunings
$d_{j}=\delta_{j}-\hbar\omega_{j}$. We displace the cavity fields in order to
remove the terms describing the cavity pumps. The displacement is described by
the operators
$\hat{D}_{j}(\alpha_{j})=e^{\alpha_{j}\hat{c}_{j}^{\dagger}-\alpha_{j}^{*}\hat{c}_{j}}$
such that
$\hat{D}_{j}^{\dagger}(\alpha_{j})\hat{c}_{j}\hat{D}_{j}(\alpha_{j})=\hat{c}_{j}+\alpha_{j}.$
In the absence of coupling to the array, the displaced cavity modes would look
like damped (but not pumped) harmonic oscillators relaxing to the vacuum. This
vacuum corresponds, for the non-displaced fields, to coherent states
$\ket{\alpha_{j}}=\hat{D}(\alpha_{j})\ket{0}$ having mean photon number
$\bar{n}_{\mathrm{ph},j}=|\alpha_{j}|^{2}$. The field displacement changes the
Hamiltonian (S66) according to

$\displaystyle\hat{H}_{\mathrm{eff}}=\hat{H}_{\mathrm{S}}^{\mathrm{eff}}+\sum_{j=1}^{L}\Big{\\{}\chi_{j}(\hat{n}_{j})\big{[}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\alpha_{j}^{*}\hat{c}_{j}+\alpha_{j}\hat{c}_{j}^{\dagger}+|\alpha_{j}|^{2}\big{]}+d_{j}\big{[}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\alpha_{j}^{*}\hat{c}_{j}+\alpha_{j}\hat{c}_{j}^{\dagger}+|\alpha_{j}|^{2}\big{]}+\mathcal{E}_{j}\hat{c}_{j}^{\dagger}+\mathcal{E}_{j}^{*}\hat{c}_{j}\Big{\\}},$
(S69)

where we have dropped constant energy shifts. Moreover, it changes the
dissipative part of the master equation as follows,

$\displaystyle\kappa_{j}\hat{D}_{j}^{\dagger}(\alpha_{j})\mathcal{D}[\hat{c}_{j}]\hat{\rho}\hat{D}_{j}(\alpha_{j})=-\frac{i}{\hbar}\Big{[}\hat{H}_{\alpha_{j}},\hat{\rho}\Big{]}+\kappa_{j}\mathcal{D}[\hat{c}_{j}]\hat{\rho},\qquad\hat{H}_{\alpha_{j}}=\frac{i\hbar\kappa_{j}}{2}(\alpha_{j}^{*}\hat{c}_{j}-\alpha_{j}\hat{c}_{j}^{\dagger}).$
(S70)

The amplitudes $\alpha_{j}$ are chosen such as to cancel all terms
proportional to $\hat{c}_{j}$ and $\hat{c}^{\dagger}_{j}$ that do not couple
to the system in Eq. (S69), which yields

$\alpha_{j}=-\frac{\mathcal{E}_{j}}{d_{j}-i\hbar\kappa_{j}/2}.$ (S71)

The master equation for the displaced density matrix is then of the form of
Eq. (S65) with displaced effective Hamiltonian

$\hat{H}_{\mathrm{eff}}=\hat{H}_{\mathrm{S}}^{\mathrm{eff}}+\sum_{j=1}^{L}\Big{\\{}\chi_{j}(\hat{n}_{j})\big{[}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\alpha_{j}^{*}\hat{c}_{j}+\alpha_{j}\hat{c}_{j}^{\dagger}+|\alpha_{j}|^{2}\big{]}+d_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{\\}}.$
(S72)

Let us recall that the Hamiltonian $\hat{H}_{\mathrm{S}}^{\mathrm{eff}}$ is
the atomic Hamiltonian that is Floquet engineered to feature artificial
magnetic flux and whose eigenstates we wish to prepare and stabilize. We
represent the Hamiltonian of Eq. (S72) in the basis of the eigenstates of
$\hat{H}_{\mathrm{S}}^{\mathrm{eff}}=\sum_{\mu}\varepsilon_{\mu}\ket{\mu}\\!\bra{\mu},$
obtaining

$\hat{H}_{\mathrm{eff}}=\sum_{\mu}\varepsilon_{\mu}\ket{\mu}\\!\bra{\mu}+\sum_{j=1}^{L}\sum_{\mu,\mu^{\prime}}\Big{\\{}\bra{\mu}\chi_{j}(\hat{n}_{j})\ket{{\mu^{\prime}}}\ket{\mu}\\!\bra{{\mu^{\prime}}}\Big{(}\hat{c}_{j}^{\dagger}\hat{c}_{j}+\alpha^{*}\hat{c}_{j}+\alpha_{j}\hat{c}_{j}^{\dagger}+|\alpha_{j}|^{2}\Big{)}+d_{j}\hat{c}_{j}^{\dagger}\hat{c}_{j}\Big{\\}}.$
(S73)

In interaction picture with respect to the first and last term, this further
becomes

$\hat{H}_{\mathrm{eff}}(t)=\sum_{j=1}^{L}\sum_{\mu,\mu^{\prime}}\bra{\mu}\chi_{j}(\hat{n}_{j})\ket{{\mu^{\prime}}}\ket{\mu}\\!\bra{{\mu^{\prime}}}\Big{[}(\hat{c}_{j}^{\dagger}\hat{c}_{j}+|\alpha_{j}|^{2})e^{i(\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}})t}+\alpha_{j}^{*}\hat{c}_{j}e^{i(\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}}-d_{j})t}+\alpha_{j}\hat{c}_{j}^{\dagger}e^{i(\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}}+d_{j})t}\Big{]}.$
(S74)

From this expression one can see that by choosing
$d_{j}=\pm(\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}})$, so by setting the
detuning of the pump from the $j$-th cavity resonance to match the effective
transition energy $\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}}$ of the
system, terms of the form
$\ket{\mu}\\!\bra{{\mu^{\prime}}}\hat{c}_{j},\,\ket{{\mu}}\\!\bra{{\mu^{\prime}}}\hat{c}_{j}^{\dagger},$
become resonant. These terms describe system transitions between the effective
eigenstates $\ket{{\mu^{\prime}}}$ and $\ket{{\mu}}$ which are accompanied by
the absorption or emission of a photon (in the displaced modes) in the $j$-th
cavity. For the other terms in the Hamiltonian not to compete with such
processes and thus being negligible in rotating-wave-approximation, it must
hold that the system effective energy gaps are much larger than the effective
atom-cavity coupling,

$|\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}}|\gg|\\!\bra{\mu}\chi_{j}(\hat{n}_{j})\ket{{\mu^{\prime}}}\\!|\sqrt{n_{\mathrm{ph},j}}\
.$ (S75)

Also, it must hold that the difference between the gaps of the resonant and
non-resonant processes are much larger than the coupling,

$|(\varepsilon_{\mu}-\varepsilon_{\mu^{\prime}})-(\varepsilon_{\eta}-\varepsilon_{\eta^{\prime}})|\gg|\\!\bra{\mu}\chi_{j}(\hat{n}_{j})\ket{{\mu^{\prime}}}\\!|\sqrt{n_{\mathrm{ph},j}},$
(S76)

where we have fixed $d_{j}=\varepsilon_{\eta}-\varepsilon_{\eta^{\prime}}$ for
some $\eta$. The latter condition is needed if one really aims at addressing a
single transition. However, in the applications studied the fact that multiple
effective transitions have a similar gap is advantageous, since they can all
be controlled with a single cavity. Under the above conditions, we can focus
on a single resonant transition, for clarity, and the cavity used to control
it. The Hamiltonian then reduces to

$\hat{H}_{\eta\eta^{\prime}j}=\overline{\chi}_{\eta\eta^{\prime}j}\ket{{\eta}}\\!\bra{{\eta^{\prime}}}\hat{c}_{j}^{\dagger}+\overline{\chi}^{*}_{\eta^{\prime}\eta
j}\ket{{\eta^{\prime}}}\\!\bra{{\eta}}\hat{c}_{j}\ ,$ (S77)

where we have defined the rate
$\overline{\chi}_{\eta\eta^{\prime}j}=\alpha_{j}\bra{\eta}\chi_{j}(\hat{n}_{j})\ket{{\eta^{\prime}}}$.
In this reduced model, the master equation now describes a Jaynes-Cummings
interaction between the two-level system comprising $\ket{\varepsilon_{\eta}}$
and $\ket{\varepsilon_{\eta^{\prime}}}$ and the cavity with photon leakage
from the cavity,

$\frac{d\hat{\rho}}{dt}=-\frac{i}{\hbar}[\hat{H}_{\eta\eta^{\prime}j},\hat{\rho}]+\kappa_{j}\mathcal{D}[\hat{c}_{j}](\hat{\rho})\
.$ (S78)

This equation can be solved in the “bad-cavity” limit,
$\hbar\kappa_{j}\gg\overline{\chi}_{\eta\eta^{\prime}j}$ [57]. The solution
features a dissipative transition from $\ket{{\eta}^{\prime}}$ to
$\ket{{\eta}}$ in which population is transferred exponentially fast,
$p_{\eta^{\prime}}(t)=e^{-\Gamma_{\eta^{\prime}\to\eta}^{(j)}t}p_{\eta^{\prime}}(0)$,
with rate

$\displaystyle\Gamma_{\eta^{\prime}\to\eta}^{(j)}=$ $\displaystyle
4\frac{\bar{n}_{\mathrm{ph},j}}{\hbar^{2}\kappa_{j}}\
\lvert\bra{{\eta}}\chi_{j}(\hat{n}_{j})\ket{{\eta}^{\prime}}\rvert^{2}.$ (S79)

The transition takes place also if the system departs from the bad-cavity
limit ($\hbar\kappa_{j}\lesssim\overline{\chi}_{\eta\eta^{\prime}j}$), but
with more complex transient dynamics [37]. This is due to the fact that, in
this case, the atom-cavity system enters an effective strong-coupling regime
and can exhibit damped coherent evolution, rather than purely monotonic decay.
In conclusion, each cavity can be exploited to induce a dissipative transition
at rate $\Gamma_{\eta^{\prime}\to\eta}^{(j)}$ between two eigenstates (more
generally, two energy resolvable subspaces) of the effective atomic
Hamiltonian obtained through Floquet engineering. In some of the examples
shown in the main text, a single cavity is used to control multiple
transitions. This is possible, and can be exploited, whenever the
corresponding transition energies differ by less than the cavity decay energy
$\hbar\kappa$. Indeed, the rate for a generic, non-resonant, transition can be
estimated via Fermi’s golden rule [35], and is proportional to the cavity
spectral density. For a pumped-damped cavity, the latter has Lorentzian shape
of width $\kappa$ centered at the cavity-pump detuning [42, 37, 35]. Hence, if
a transition is activated resonantly by a cavity, other transitions with an
energy gap in a nearby window of width $\hbar\kappa$ will be activated with a
sizeable rate, since they still fall within the width of the Lorentzian peak.

##  III Numerical simulations, postselection, and further examples

###  III.1 Numerical simulations

Numerical simulations of the full driven model and of the effective model are
done using the master equations (S9) and (S64), respectively, in both cases
after having performed an additional transformation. Indeed, in order to
reduce the number of cavity states to be included in the truncated Hilbert
space, a displacement of all the cavity modes is first performed to eliminate
linear terms in the cavity operators (not coupling to the array) from the
master equation, following the same steps done in Sec. II. When pumped, each
cavity initially relaxes to the coherent state described by the displacement
operator, corresponding to the vacuum of the displaced fields, over a
timescale much faster than that characterizing the dynamics we want to study.
Therefore, in the simulations each cavity is initialized directly in the
vacuum of the displaced field and states with up to three photons are included
per cavity. For the array, the driven dynamics with $N$ total excitations is
simulated including states with up to $N+1$ total excitations, since there are
no non-zero matrix elements associated with the simultaneous creation or loss
of more than a single extra excitation. Only states up to $N$ excitations
instead are kept for the effective master equation, since higher excitation
states are completely decoupled. All the Floquet-dressed coefficients entering
the effective Hamiltonian, Eq. (4) in the main text, are computed with the
truncation $m\in[-50,50]$. The same truncation applies to the effective jump
operators, Eq. (5). In all simulations expect for those reported in Fig. 3(d)
(concerning the two-excitation example), the master equation is directly
integrated numerically in time. The results of Fig. 3d are computationally
more demanding, and are obtained instead by using the Monte Carlo wavefunction
method [58] and averaging over 350 trajectories. The code for solving the
master equations in time makes extensive use of the QuTip Python library [59,
60].

###  III.2 Postselection

We here discuss further how it is possible, in practice, to probe and verify
interesting properties characterizing the states prepared, despite the weak
excitation loss due to Purcell decay that could in principle spoil the quantum
simulation, by simply postselecting on the measurement results. In general,
one can postselect by measuring observables that carry also information about
the total number of excitations, such that a record can be discarded if the
latter was not conserved. The site densities and currents that we present in
Figs. 3 and 4 belong indeed to this class, and are at the same time key
signatures characterizing the states stabilized. Site densities are measured
directly through standard dispersive readout of all atoms, from which the
total number of excitations is also determined allowing postselection. These
quantities immediately allow one to observe the typical density patterns of
flux ladders, which are more generally an important signature of quantum-Hall-
like states. Also, they allow observation of the signature interference
pattern in the Aharonov-Bohm cage example, Sec. III.4. Currents can be
measured using the method proposed in Ref. [41] and experimentally realized in
cold atoms, e.g., in [11]. This involves biasing the atom frequencies (which
is the control knob we use to perform Floquet engineering) before measuring
the excitation number of each atom via dispersive readout. More in detail, one
performs a sequence of measurements in which pairs of neighbouring atoms are
first detuned from other atoms, suppressing tunneling between different pairs,
and then site occupations are measured. The time evolution of the latter,
obtained by measuring after different idle times, will exhibit oscillations
which are directly determined by the current between the two sites, allowing
its estimation [41]. Since the frequency biasing does not change the atomic
excitations, and the measurements also give the total excitation number,
postselection can be performed while currents are measured.

In numerical simulations, where we can access the full density matrix
$\hat{\rho}$, postselection in reproduced by projecting the density matrix in
the subspace of conserved excitations (defined by the projector $\hat{P}$)
when computing observables,
$\rho_{\mathrm{ps}}=\hat{P}\rho\hat{P}/\mathrm{tr}[\hat{P}\hat{\rho}]$, which
identifies the discarded fraction as $p=1-\mathrm{tr}[\hat{P}\hat{\rho}]$. In
Fig. 3 and 4 we show the occupation of the effective eigenstates computed in
this way, despite being hardly measurable directly, in order to unambiguously
certify the success of the protocols, in addition to showing measurable
observables such as site occupations and currents.

Figure S2: Ground state preparation for an excitation in a three-plaquette
ladder. (a) System geometry and cavity couplings. (b) Effective energy
spectrum and transitions induced via reservoir engineering. (c) Build up of
population in the effective ground state (solid, for the full driven model,
dashed, for the effective master equation) induced by the Floquet-dissipative
scheme. The grey shaded area indicates the fraction of discarded data in
postselection in the full model, with the corresponding prediction given by
the effective master equation (black dashed) (d) Final current patterns (the
size of the arrows reproduces the current strength in arbitrary units) and
site densities (color scale) in the final state of panel (c) [upper panel] as
compared to the ideal ground state [lower panel]. The parameters used are
$\Phi=0.4\pi$, $\hbar\omega=20J$, $\lambda=0.72\omega$,
$\bm{\delta}/\hbar\omega=(1.74,1.71,1.764)$, $U=8J$,
$\bm{\kappa}/J=(0.09,0.09,0.09)$, $\bm{\mathcal{E}}/J=(1,1.7,0.4)$,
$\bm{g}/J=(1,1,1)$.

###  III.3 More on the stabilization in the bosonic ladder

Considering the ladder system of Fig. 3 of the main text, we label the sites
with indices $(\ell,r)$ indicating the leg ($\ell=1$ for the upper leg,
$\ell=2$ for the lower), and the rung, $r=1,2,3$. The current on the edges
along the legs ($j_{\ell;r,r+1}^{\leavevmode\hbox to5.32pt{\vbox
to4.47pt{\pgfpicture\makeatletter\hbox{\hskip 0.09999pt\lower-0.95354pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{3.41418pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.85355pt}{-0.85355pt}\pgfsys@lineto{5.12128pt}{2.56064pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}$) and along the
rungs ($j_{r}^{\perp}$) are computed as the expectation value of the operators
[43, 44, 11]

$\displaystyle\hat{j}_{\ell;r,r+1}^{\leavevmode\hbox to5.32pt{\vbox
to4.47pt{\pgfpicture\makeatletter\hbox{\hskip 0.09999pt\lower-0.95354pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{4.26773pt}{3.41418pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.2pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.85355pt}{-0.85355pt}\pgfsys@lineto{5.12128pt}{2.56064pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}=iJ(\hat{a}_{\ell,r+1}^{\dagger}\hat{a}_{\ell,r}-\hat{a}_{\ell,r+1}\hat{a}_{\ell,r}^{\dagger}),$
(S80a)
$\displaystyle\hat{j}_{r}^{\perp}=iJ(e^{-ir\Phi}\hat{a}_{1,r}^{\dagger}\hat{a}_{2,r}-e^{ir\Phi}\hat{a}_{1,r}\hat{a}_{2,r}^{\dagger}).$
(S80b)

A further example of ground state preparation in a small bosonic ladder is
reported in Fig. S2. In this case, a single excitation in a three-plaquettes
is considered, coupled to three cavities. A final current pattern in agreement
with the ideal one is achieved. Better agreement is obtained at longer
stabilization times.

###  III.4 Preparing Aharonov-Bohm cages

Figure S3: (a) Rhombic geometry and (b) single-excitation energy spectrum of
the effective array Hamiltonian. The cavity couples at the central site and
cools the state component in the most excited state to the ground state. (c)
Final localized site density exhibiting Aharonov-Bohm caging. The parameters
used are $\hbar\omega=25J$, $\delta_{1}=34.8J$, $U=12J$,
$\kappa_{1}=0.2J/\hbar$, $\mathcal{E}_{1}=2.7J$, $g_{1}=J$, final time
$Jt/\hbar=600$.

We here consider the diamond-chain of corner-sharing rhombic plaquettes with
flux $\Phi=\pi$ [Fig. S3(a)]. It exhibits so called Aharonov Bohm cages,
single-excitation states that are strictly localized on each x-shaped group of
five sites via destructive interference [50, 33, 51]. While these states form
a flat Bloch band in an extended chain, one of them exists already as the
ground state of a minimal model with two plaquettes. It can be prepared in a
controlled fashion by starting from a single excitation on the central site
(having equal overlap with both the most excited and the ground state) using a
single resonator, as shown in Fig. S3(b). The site occupations after the
preparation are depicted in Fig. S3(c), confirming the achievement of the
desired effect. In a large chain this preparation scheme can easily be scaled
up to prepare localized excitations at targeted positions.

##  IV Choosing the parameters and an illustrative example

In this section we discuss further and exemplify how the parameters for the
driven-dissipative stabilization protocols are chosen.

###  IV.1 General conditions on parameters

The values of the fundamental system parameters are chosen to be typical of
state-of-the-art circuit QED platforms. In particular, we probed values such
as $U\sim 8-15J$, $g_{j}=J$, $\kappa_{j}\sim 0.05-0.2J$, $\delta_{j}\sim
30-50J$ having in mind tunneling rates in the range $J/2\pi\hbar\sim 20-50$
MHz. We then impose that the constraints required by Floquet and reservoir
engineering discussed in the main text are satisfied. Summarizing, this means

  1. 1.

$\hbar\omega\gg J:$ The results shown in the article are obtained choosing
$\hbar\omega=20J$. This leads to single-excitation spectra of the Floquet-
engineered effective Hamiltonian for the atoms with bandwidth of a few units
of $J$.

  2. 2.

$|\delta_{j}+Un_{j}+m\hbar\omega|\gg g_{j}$. For the given values of
$\hbar\omega$ and $U$, the values of the atom-cavity detunings $\delta_{j}$
are fine tuned by changing the atomic potential $\Delta$ to ensure that
Floquet resonances are avoided.

  3. 3.

$|\varepsilon_{\eta}-\varepsilon_{\eta^{\prime}}|\gg\hbar\kappa_{j}\gg|\\!\bra{\eta}\chi_{j}(\hat{n}_{j})\ket{{\eta^{\prime}}}\\!|\sqrt{\bar{n}_{\mathrm{ph},j}}$:
for the engineered dissipative processes to resolve single effective
transitions, the effective gaps must be much larger than the decay energy
$\hbar\kappa_{j}$ and than the effective atom-cavity coupling. However, for
the cavities to act as a Markovian bath, the decay rate must in turn be much
larger than the effective coupling. In few-plaquette systems, the effective
gaps are of the order of the effective tunneling rate
$J^{\mathrm{eff}}_{\ell\ell^{\prime}}\simeq 0.5J$. The effective atom-cavity
coupling $\chi_{j}(\hat{n}_{j})$ is a second order process and is in general
smaller than these gaps, away from Floquet resonances. Therefore, at fixed
$\hbar\kappa_{j}\ll J$ and given that $d_{j}$ is fixed by the gap to be
controlled, the conditions stated here mainly constrain the pump amplitude
$\mathcal{E}_{j}$ for the cavities. These amplitudes should be as large as
possible to enhance the transition rate, but small enough to satisfy the above
conditions. In general, we find optimal transition rates for amplitudes
leading to $\bar{n}_{\mathrm{ph},j}\sim 1-2$, consistently with what was
observed in Ref. [35] for an undriven system.

  4. 4.

Another degree of freedom is given by the “position” of the cavities, namely
to which sites they are coupled. After inspecting the numerical values of the
matrix elements of the coupling operators
$\bra{\eta}\chi_{j}(\hat{n}_{j})\ket{{\eta^{\prime}}}$, the choice is done in
order to maximize the effective coupling for a desired transition.

Using these conditions as guiding lines, the search for good parameters is
done numerically. An example of this search and of how the dissipative control
protocols are designed is illustrated with a small system in the following.

Figure S4: (a) Sketch of the single-plaquette system and dissipative
transition induced by the cavity. (b) a process in which one photon tunnels to
the cavity by emitting two drive quanta $\hbar\omega$ in the Floquet drive and
leaks to the environment, leading to excitation loss. Such processes spoil the
efficiency of the Floquet and reservoir engineering as visible from panel (c).
(c) effective ground state occupation (color scale) as a function of time and
of the atom-cavity detuning $\delta$ for $U=15.5J$. Left and right panels
compare the results with the full master equation and the results predicted by
the effective master equation, respectively. The other parameters are
$\Phi=\pi/2$, $g=J$, $\kappa=0.1J$, $\mathcal{E}=1.6J$, $\hbar\omega=20J$,
$\lambda\simeq 0.72\omega$ (the value of $\lambda$ is such that the effective
tunneling strength along every edge is the same, see Section “Artificial
magnetic flux”). The upper panel corresponds to non-postselected results,
while the lower panel shows postselected results. (d) same as in panel (c) for
a different interaction value $U=8J$. Thin red lines in the panels “effective”
in (c) and (d) approximately indicate points where the perturbative expansion
diverges, namely they represent regions where one of the denominators in the
effective coupling becomes smaller than a numerical threshold $10^{-7}J$.

###  IV.2 Single dissipative transition

Consider one excitation circulating in an array of four artificial atoms
making up a single square plaquette pierced by an artificial magnetic flux
$\Phi$ and coupled to a single cavity, as depicted in Fig. S4(a). In this
case, due to the symmetry of the system, the position of the cavity is
irrelevant and leads to the same (absolute value of the) matrix elements for
the system-cavity coupling. Figure S4 characterizes an ideal transition from
the effective second-excited eigenstate $\ket{2}$ of the effective atomic
Hamiltonian to the ground state $\ket{0}$. This is realized by setting the
cavity-pump detuning to match the energy difference
$d=\varepsilon_{2}-\varepsilon_{0}$, then initializing the system in $\ket{2}$
in the simulation and studying the stroboscopic occupation of $\ket{0}$ as a
function of time for varying atom-cavity detuning $\delta$, with other
parameters fixed. Panels (c) and (d) compare the prediction of the effective
master equation with the full dynamics for two values of the interaction
parameter $U$, in (c) also comparing postselected and non-postselected
results. First of all, one can see that the effective model, which is much
less computationally demanding, can be reliably used for parameter
exploration. Low-probability (white) lines indicate divergences in the
perturbative expansion. These divergences correspond to drive-assisted
tunneling between the system and the cavity. For instance, as sketched in
panel (b) of Fig. S4, the white region at $\delta\sim 2\hbar\omega$ in panel
(c) corresponds to a process in which the excitation tunnels towards the
cavity while emitting two Floquet quanta in the drive. The excitation is
eventually lost due to cavity decay.

The effective model is derived using unperturbed gaps that do not include the
presence of the tunneling $J$ in the system, which is included instead in the
“diagonal” perturbation $\hat{\mathcal{V}}_{\mathrm{D}}$, see Eq. (S27a). The
resonances that are visible in the panels “effective” of Fig. S4(c)-(d) then
split into a larger number of resonances in the full model, which will still
be separated by a gap much smaller than $\delta$, since $J\ll\delta$.
Moreover, higher-order processes not captured in the effective model will
contribute to this effect and further introduce small energy shifts of the
resonances. All these effects yield the more complex pattern of “white
resonances” visible in the panels “full driven” of S4(c)-(d). Although these
resonances can make the search difficult, good parameter regimes can still be
found by carefully avoiding regions where the approximations underlying the
effective model break down.

###  IV.3 Autonomous stabilization

We now illustrate in the single-plaquette system how protocols for autonomous
stabilization are designed. In particular, we stabilize the first-excited
state $\ket{1}$ by coupling the system to a second cavity and exploiting both
“down-in-energy” and “up-in-energy” transitions, as depicted in Fig.
S5(a)-(b). The first cavity remains set to cool transitions with gap
$\varepsilon_{2}-\varepsilon_{0}$, as in the previous subsection. Due to
symmetries in the spectrum, the first cavity thus “cools” from
$\varepsilon_{2}$ to $\varepsilon_{0}$ and from $\varepsilon_{3}$ to
$\varepsilon_{1}$. The second cavity-pump detuning is set instead at the
(negative) gap $\varepsilon_{0}-\varepsilon_{1}$, thus inducing “heating”
transitions from $\varepsilon_{0}$ to $\varepsilon_{1}$ and from
$\varepsilon_{2}$ to $\varepsilon_{3}$. The fact that $\ket{1}$ is the unique
steady state can intuitively be seen by noting that no transition can occur
that leaves such a state, while a sequence of transitions can connect any
other state to $\ket{1}$. The stabilization is shown in Fig. S5(c), which
depicts the time evolution of the eigenstate occupations (using the full
master equation and postselecting – see Section “Numerical simulations and
postselection”) starting from an infinite-temperature state in the single-
excitation subspace, i.e.,
$\hat{\rho}_{0}=\frac{1}{M}\sum_{\ell=1}^{M}\hat{a}_{\ell}^{\dagger}\ket{0}\bra{0}\hat{a}_{\ell}=\mathbb{1}_{\mathrm{1p}}/M$.
The system successfully evolves towards $\ket{{1}}$. Note that such a
transition from a fully mixed state to a (nearly) pure state, which entails an
entropy decrease, could not be achieved in the absence of dissipative
processes.

Figure S5: Autonomous stabilization. (a) Two cavities are attached to the
upper two sites. (b) Dissipative transitions engineered using the two
cavities. The dissipative path pushes the system towards the first excited
state of the effective Hamiltonian. (c) Effective-eigenstate populations as a
function of time, starting from an infinite-temperature state in the single-
excitation manifold. The parameters used are $g_{j}=J$, $U=15.5J$,
$\Delta=2.25\omega$, $\hbar\omega=20J$, $\kappa_{j}=0.1J$,
$\bm{\mathcal{E}}=(1.8J,J)$.