# Bulk properties of honeycomb lattices of superconducting microwave
resonators

Alexis Morvan  Laboratoire de Physique des Solides, CNRS, Université Paris
Saclay, Orsay, France  Mathieu Féchant  Laboratoire de Physique des Solides,
CNRS, Université Paris Saclay, Orsay, France  Gianluca Aiello  Laboratoire de
Physique des Solides, CNRS, Université Paris Saclay, Orsay, France  Julien
Gabelli  Laboratoire de Physique des Solides, CNRS, Université Paris Saclay,
Orsay, France  Jérôme Estève  Laboratoire de Physique des Solides, CNRS,
Université Paris Saclay, Orsay, France

###### Abstract

We have realized different honeycomb lattices for microwave photons in the 4
to $8\text{\,}\mathrm{GHz}$ band using superconducting spiral resonators. Each
lattice comprises a few hundred sites. Two designs have been studied, one
leading to two bands touching at the Dirac points and one where a gap opens at
the Dirac points. Using a scanning laser technique to image the eigenmodes of
this new type of photonic lattices, we are able to reconstruct their band
structure. The measured bands are in excellent agreement with _ab initio_
models that combine numerical simulations of the electromagnetic properties of
the spiral resonator and analytical calculations.

Superconducting photonic lattices for microwave photons hold the promise to
simulate the behaviour of strongly interacting bosons in tailored 1D and 2D
lattices Koch and Le Hur (2009); Koch  _et al._ (2010); Schmidt and Koch
(2013). In such systems, the excitation are microwave photons with a typical
frequency around $6\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ stored in
superconducting resonators with high quality factor. The main interest of this
platform is the possibility to reach the strongly interacting regime using
Josephson junctions as a non-linear element. But, Josephson junctions are
prone to disorder, and the realization of large disorder free non-linear
lattices remains a very difficult task. Nevertheless, many-body effects have
already been demonstrated in small uni-dimensional lattices Roushan  _et al._
(2017); Fitzpatrick  _et al._ (2017); Ma  _et al._ (2018); Guo  _et al._
(2020); Fedorov  _et al._ (2020).

Independently of controlling the interactions in the lattice, it seems also
interesting to develop techniques allowing to probe lattices in the linear
regime and to understand their properties Underwood  _et al._ (2016); Kollár
_et al._ (2019). This research direction, towards the realization of
engineered band structures for photons, catches up with the rapidly growing
field of photonic lattices and topological photonics Lu  _et al._ (2014);
Khanikaev and Shvets (2017); Ozawa  _et al._ (2018). In comparison to other
photonic systems, the cryogenic environment and the low energy of the photons
bring some technical difficulties to characterize the properties of such
lattices.

In this manuscript, we use a laser scanning imaging technique in order to map
the spatial distribution of the resonant modes across the lattice. This
technique is similar to the laser scanning microscopy developed to observe the
current density in superconductor films Culbertson  _et al._ (1998) and has
already been used to characterize single a superconducting resonator Zhuravel
_et al._ (2012); Averkin  _et al._ (2013) or coupled resonators Wang  _et al._
(2019). It is an alternative to the scanning technique developed in Underwood
_et al._ (2016). From a Fourier analysis of the spatial distribution of the
modes, we reconstruct the dispersion relation of the lattice bands. We apply
this method to three honeycomb lattices with different designs, where the $A$
and $B$ sites in the unit cell correspond to identical or different
resonators. In one case, we observe two bands touching at the Dirac point,
while in the latter, a gap opens at the Dirac points. These properties of the
bands are similar to the ones for electrons in graphene, even though the two
systems are not described by the same wave equation. In order to understand
our lattices and pinpoint the difference with electronic systems, we have
developed two approaches to predict the properties of the lattice. The first
approach consists in projecting the Maxwell equations on the basis of the
resonator modes, which corresponds to a coupled mode theory (CMT), a technique
widely used in photonics Huang (1994). The band structure of the lattice is
then obtained in terms of overlap integrals between the electric and magnetic
fields of neighbouring resonators. The second approach relies on more
intensive numerical simulations to obtain an equivalent circuit to the
lattice, from which the band structure can be calculated. These two _ab
initio_ approaches allow us to reproduce the measured band structures with a
very good accuracy.

The manuscript is organized as follows. In the first part, we present the
properties of the spiral resonator that is used as the site of the different
lattices. In the second part, we present the measurement of the band structure
of the three lattice designs. And the last part details the two models that we
have developed to predict the band and mode structures of the lattices.

##  I Spiral resonator properties

The three different lattices studied in this manuscript are honeycomb lattices
where each site consists of a superconducting spiral resonator as shown in
figure 1. The spiral is made of a $4.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$
wide Nb wire that is winded in an hexagonal pattern with an overall size of
$300\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The lattices are patterned
through photo-lithography and reactive ion etching starting from a
$\sim$300\text{\,}\mathrm{nm}$$ thick sputtered Nb layer on top of a silicon
wafer. The design of the spiral was adjusted through numerical simulations
such that the resonance frequency $\omega_{0}$ of the fundamental mode is
close to $2\pi\times 6\text{\,}\mathrm{GHz}$. This frequency can be finely
tuned by adjusting the length of the wire at the center of the spiral. For the
experiments shown in this manuscript, we use the two designs shown in figure
1a,b that were chosen to obtain two sites with slightly different $\omega_{0}$
spaced by $2\pi\times 120\text{\,}\mathrm{MHz}$. The figure 1c shows the
simulated self-admittance of the two designs, from which we deduce
$\omega_{0}$ as well as the resonance frequencies of the higher order modes.
The second mode is expected to resonate around $14\text{\,}\mathrm{GHz}$ and
plays no role in the experiments presented here, where all measurements are
performed between 5 and $7\text{\,}\mathrm{GHz}$.

Figure 1: (a,b) Design of the two spirals used as the sites of our honeycomb
lattices. The spiral shown in (a) consists of a $8.6\text{\,}\mathrm{mm}$ long
wire with a width of $4.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and a gap
between adjacent turns of $8.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$. The
length of the (b) spiral is shorter resulting in a slightly higher resonance
frequency. (c) Self-admittance of the two spirals calculated for a port
located at the end of the wire at the center of the spiral. The resonances
correspond to the zero-crossings. The inset shows a zoom close to the
fundamental resonance.

When two spirals are approached at a distance $d$ as shown in figure 2, the
fundamental modes couple, giving rise to two resonances at frequencies
$\omega_{-}$ and $\omega_{+}$. If the coupling is not too strong, the coupled
modes can be expressed as a linear combination of the uncoupled modes. This
approach that consists in using a reduced set of well chosen modes as a basis
is known as the Coupled Mode Theory (CMT)Powell  _et al._ (2010); Lomanets
_et al._ (2012); Elnaggar  _et al._ (2015). Here, we use a basic CMT with only
one mode per resonator. The Maxwell equations projected in this basis lead to
a linear system, whose eigenvalues correspond to the resonance frequencies
$\omega_{-}$ and $\omega_{+}$. If the two resonators are identical with equal
$\omega_{0}$, one obtains:

$\omega_{\pm}=\sqrt{\frac{1\pm\kappa_{e}}{1\pm\kappa_{m}}}\ \omega_{0}$ (1)

where $\kappa_{e}$ ($\kappa_{m}$) is the electric (magnetic) coupling
constant. These coupling constants are proportional to the overlap of the
electric (magnetic) fields of the two resonators. We define

$\displaystyle D_{ij}$
$\displaystyle=\int\epsilon(\mathbf{r})\,\mathbf{E}_{i}(\mathbf{r})\cdot\mathbf{E}_{j}(\mathbf{r})d^{3}\mathbf{r}$
(2) $\displaystyle G_{ij}$
$\displaystyle=\int\mu_{0}\,\mathbf{H}_{i}(\mathbf{r})\cdot\mathbf{H}_{j}(\mathbf{r})d^{3}\mathbf{r}$
(3)

where $E_{i}$ ($H_{i}$) is the electric (magnetic) field spatial dependence of
the mode associated with the $i=1,2$ resonator. We assume these functions to
be real. We explain in Appendix A how we compute these integrals from the
charge and current distribution in the spiral that is obtained from the
numerical simulation of a single resonator. In particular, one has to take
into account charge and current images due to the dielectric interface at the
surface of the sample and to the metallic ground plane below the sample. The
coupling constants are then given by $\kappa_{e}=D_{12}/D_{11}$ and
$\kappa_{m}=G_{12}/G_{11}$. Figure 2a shows the evolution of $\kappa_{e}$ and
$\kappa_{m}$ as a function of $d$ for two spirals shown in fig. 1a. Because
$\kappa_{e}>0$ and $\kappa_{m}<0$, the magnetic and the electric couplings add
up to increase the mode splitting, while the mean frequency
$(\omega_{+}-\omega_{-})/2$ remains close to $\omega_{0}$. At large distances
($d>$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$$), the magnetic coupling
dominates, while at short distances both couplings are important.

Figure 2: Coupling between two spiral resonators as shown in 1a a) Electric
and magnetic coupling constants (see main text and Appendix A) as a function
of $d$. b) Resonance frequencies $\omega_{+}$ and $\omega_{-}$ of the two
coupled resonators as a function of $d$. Red dots correspond to the result of
a numerical simulation that computes the admittance matrix of the two coupled
resonators. The coupled mode theory prediction only uses the simulation of the
electric and magnetic field created by a single resonator. It fails at short
distances, where the resonator field is perturbed by the presence of the other
resonator.

We compare the CMT predictions to a full numerical simulation of the two
coupled resonators. From the simulation, we obtain the $2\times 2$ admittance
matrix $Y(\omega)$, which corresponds to the admittance matrix for two ports
located at the ends of the two spirals (see figure). The two resonances
$\omega_{+}$ and $\omega_{-}$ are then obtained as the zeros of $\det
Y(\omega)$. The CMT predictions coincide with the ones of the full numerical
simulation at large distances ($d>$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$$)
but fails at short distances as shown in figure 2b. This is because, at short
distances, the coupling is too strong and the coupled modes are not linear
superposition of the uncoupled modes. In the following, we show results for
lattices where $d=$5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ or
$d=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. We therefore expect lattices
with $d=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ to be well described by
the CMT method, while lattices with $d=$5\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$ will require a more advanced model including the simulation of coupled
sites. The half splitting $(\omega_{+}-\omega_{-})/2$ obtained from the
admittance matrix method is $2\pi\times 200\text{\,}\mathrm{MHz}$ for
$d=$5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and $2\pi\times
120\text{\,}\mathrm{MHz}$ for $d=$30\text{\,}\mathrm{\SIUnitSymbolMicro m}$$.
This number gives an estimate of the nearest neighbour coupling amplitude in a
tight-binding description of the lattice.

##  II Lattice measurements

###  II.1 Density of states

We have built honeycomb lattices consisting of a few hundred spiral resonators
that fit on a $20\times 10\text{\,}{\mathrm{mm}}^{2}$ sample, which is then
mounted on the $1\text{\,}\mathrm{K}$ stage of a dry dilution fridge. We have
characterized three designs of honeycomb lattice labeled G,SI and SII, as
shown in figure 3. The G design corresponds to the situation where all sites
are identical, as in graphene, and are occupied by the spiral shown in figure
1a. The other two designs, SI and SII, correspond to the more general
situation where the two nonequivalent A and B sites of the honeycomb lattice
are occupied by the two different spirals shown in figure 1a and 1b. Because,
to first approximation, the sites only differ by their onsite frequency, the
SI and SII designs realize a so-called Semenoff insulator Semenoff (1984). In
the graphene case, the expected band structure consists of two bands touching
at the two Dirac points, while in the Semenoff case a gap opens at the Dirac
points. In the two Semenoff designs, a horizontal boundary divides the sample
in two halves: The A (B) sites in the lower half are occupied by the
resonators that occupy the B (A) sites in the upper half. The two halves
correspond to the same infinite lattice and therefore have the same bulk
properties. The role of this boundary is to create valley Hall boundary states
as discussed in Morvan  _et al._ (2020).

Figure 3: Design of the three different honeycomb lattices studied in this
article. In the upper G design, all the sites are identical, while in the SI
and SII designs, the A and B sites of the lattice correspond to two different
resonators. Four coplanar waveguides are connected to single sites located on
the edges of the lattice to probe the sample.  Figure 4: Honeycomb lattice
with the two different sites A and B. The differents circles indicates the
1st, 2nd and 3rd neighbors sites.

We probe the transmission through the lattice using four microwave ports that
are connected to four coplanar waveguides. Each waveguide is coupled to a
single site located on the edge of the lattice (see figure 3). We first
characterize the bulk properties of the lattice by estimating the density of
states (DOS) that we deduce from transmission measurements. Figure 5 shows a
typical transmission spectrum obtained through a lattice of type SII. The
lattice modes are visible as sharp resonant peaks. Depending on the sample, we
identify between 43% (G design) and 80% (SI design) of the expected
resonances. This number is limited by our signal to noise ratio and by the
finite width of the peaks. Peaks that are too weakly coupled to the
measurement ports are not identified as resonances and peaks too close in
frequency are counted as a single peak. By counting the number of peaks in a
frequency window of 15 MHz for the design G and 15 MHz for the design SI and
SII, we can however obtain an estimate of the DOS as shown in figure 5. As
expected, we observe two bands for all three samples, with a clear gap for the
SI and SII samples. We can also identify the two van Hove singularities
corresponding to the two maxima in the DOS of each band. The two _ab initio_
models to which we compare our data in the figure 5 only have a global offset
frequency and are detailed in section LABEL:sec:models. The frequencies of the
remarkable points observed in the DOS (band minima, maxima and Dirac point)
are well reproduced, but the measured DOS is not in quantitative agreement
with the model predictions because a significant fraction of the modes is
missing from the measured DOS.

Figure 5: a) Transmission $|S_{xx}|^{2}$ for the SII sample. The lattice modes
appear as sharp resonant peaks. The transmission on resonance is way below one
indicating that the modes are under coupled: intrinsic loss dominate the
coupling loss to the measurement waveguides. b,c,d) The blue data show the
estimated DOS from transmission spectra similar as the one shown in a. The
data are compared to the predictions of the two _ab initio_ models detailed in
section III. The CMT model only uses the simulation of a single resonator,
while the $Y$ matrix model relies on the simulation of the admittance matrix
of a small lattice. The discrepancy between the data and the model is due to
the fact that we miss some peaks and underestimate the DOS.

###  II.2 Mode imaging and dispersion relation

In order to further characterize the lattice, we use a laser scanning
technique to map the spatial variation of the modes identified in the
transmission spectra. This imaging technique allows us to obtain a partial
information on the dispersion relation of the lattice modes. The measurement
consists in monitoring the transmission of one (or many) mode while scanning a
laser spot across the lattice. The experimental setup is shown in figure 6.
This method has been previously used to map the spatial profile of the
resonant modes of a single resonator Zhuravel  _et al._ (2012) or a chain of
resonator. The optical setup was designed to obtain a laser waist on the
sample of $60\text{\,}\mathrm{\SIUnitSymbolMicro m}$, which is much larger
than the size of the spiral wire width and spacing but smaller than the
overall resonator size. This allows us to average the mode distribution over
each site, while keeping sufficient resolution to resolve adjacent sites. We
observe that the main effect of the laser is to induce dissipation on the
illuminated site as observed in Zhuravel  _et al._ (2006, 2012). Because the
modes are under-coupled to the probe ports, the laser induced increase of the
loss results in a decrease of the mode transmission. To first order, the
transmission drop for a given mode is proportional to the mode squared
amplitude averaged over the illuminated area. In order to get rid of slow
drifts and improve the signal to noise ratio, the laser intensity is modulated
at a few k Hz . We then digitally demodulate the transmitted signal measured
with a VNA and record the amplitude of the in phase signal as a function of
the laser position.

Figure 6: Mapping the spatial dependence of the lattice modes a,b,c) Cryogenic
laser scanning setup. A collimated laser beam is mounted on a motorized mirror
mount outside the dry refrigerator (a). The beam is relayed by two lenses to
pass through openings in the $50\text{\,}\mathrm{K}$ and the
$4\text{\,}\mathrm{K}$ shields (b). A final lens focuses the beam onto the
sample that is clamped on the $1\text{\,}\mathrm{K}$ stage (c). The tilt
motion of the outer mirror mount results in a translation of the focused laser
spot on the sample. d) Fine scan of a single resonator in the lattice. The
image shows the transmission of one resonant peak as a function of the laser
position. We compare this image to a convolution of the current intensity of
the first mode of our spiral with a gaussian spot of waist $w$. Our waist size
allows to discriminate each spiral and homogenise the response when the laser
spot is on the resonator.

Figure 6d shows a fine scan of one lattice site, which appears as a blurred
hexagon with a central dip. This is consistent with the fact that the laser
induced loss is maximal where the current density is large Zhuravel  _et al._
(2006) as shown by the simulation in figure 5d. We attribute a single value
for the mode intensity per site by averaging over a few measurement points
well inside the hexagon surrounding one site. We have checked that the final
result is rather insensitive to the details of the averaging procedure. In
order to optimize the acquisition time, different scanning techniques have
been tested, continuous (as shown in figure 5e) or raster. We obtain best
results with a raster scan consisting of six measurement points per site,
while we monitor the transmission change of all the peaks in a frequency
window of about $1\text{\,}\mathrm{GHz}$. This scan method allows us to image
tens of modes in a single scan through the lattice.

Figure 7: Measured mode spatial dependence for the G sample. The black and
white images show the raw data coming from the measured transmission change at
the frequency indicated above each image, which corresponds to a resonance
peak identified in the transmission spectrum. This image is proportional to
the mode intensity ac cross the lattice. The red and blue image shows the
reconstructed signed mode amplitude (see main text and Appendix C). Dark blue
(red) indicates a large negative (positive) amplitude.  Figure 8: Measured
mode spatial dependence for the SI and SII samples. Here, we only show the
reconstructed signed mode amplitudes (see figure 7). The frequency of the mode
is indicated above each image.

Figures 7 and 8 show the results of the mode imaging technique applied to some
modes of $G$, $SI$ and $SII$ samples. In order to attribute a wavevector to
each measured mode, we have developed a reconstruction technique to deduce the
signed mode amplitude from the measured intensity. We do so by supposing that
the mode is a linear combination of three known basis modes. As a basis, we
choose the modes expected for a lattice having the same finite size geometry
as the measured and described by a tight-binding model with only nearest
neighbour coupling. We then look for the combination that matches best the
measured mode intensity and use its sign to attribute a sign to the data.
Details of the method are given in Appendix B. We then take the Fourier
transform of the reconstructed signed mode amplitude. In the reciprocal space,
several peaks lying on a circle are observed. This allows us to attribute a
single value $k$ corresponding to the norm of a wavevector to each mode.
Figure III shows the dispersion relation $\omega(k)$ that we obtain through
this analysis for the $G$ and $SI$ samples. The measured dispersion for the
$SII$ sample is almost identical to the one of the $SI$ sample and is not
shown here. We recover the two asymmetric bands observed in the DOS
estimation. The mode dispersion is quadratic at small $|k|$ leading to an
effective mass for the photons in the lattice which is on the order of () for
the lower (upper) band. At larger $k$, the dispersion relation clearly
deviates from a quadratic behaviour. For the G sample, we expect a linear
dispersion around the Dirac points but we are not able to image a sufficiently
large number of modes in this region to clearly reproduce this behaviour. This
is due to the finite size of the sample and also to the radial averaging over
the direction of the wavevector. For the SI sample, we observe that bands
curve again in the opposite direction close to the Dirac points.

Figure 9: Band structure of the G and SI samples. For a given resonance
frequency, the corresponding norm of the wavevector associated to the mode is
obtained from the Fourier analysis of the measured spatial dependence as shown
in figures 7 and 8. The measured dispersion relation is compared to the
predictions of the two _ab initio_ models detailed in section III (see also
figure 5). The theoretical predictions are plotted as shaded areas, because,
for a given $|k|$, the models predict different resonance frequencies
depending on the orientation of the wavevector.

##  III Lattice models

In comparison to the DOS, the measured dispersion shown in is less affected by
the fact that we are not able to probe all the modes. It allows us to
precisely compare our data with two _ab initio_ models. Following the analysis
of the coupling of two spiral resonators, we extend the CMT and the admittance
matrix model to the case of an infinite lattice. The two models have no free
parameters and differ in the following way: the CMT model solely relies on the
simulation of the charge and current distribution of a single spiral at
resonance, while the admittance matrix model simulates the admittance matrix
of a cluster of a few coupled sites. As explained in I, the CMT model gives a
clear physical picture of the coupling in the lattice in terms of overlap
integrals between the different sites but fails when the coupling is too
strong. As a consequence, the band structure of the S lattice is well
predicted by the CMT but observe a discrepancy for the G lattice, which is
more pronounced for the lower band. The $Y$ matrix calculation gives a more
exact description of the lattice at any coupling but requires a more intensive
numerical computation.

###  III.1 Coupled mode theory

We follow the derivation of Elnaggar  _et al._ (2015) and adapt it to the
specific case of coupled superconducting resonators. We only consider one mode
per resonator and look for a solution to the Maxwell equations at a frequency
$\omega$ as

$\displaystyle\mathbf{E}(\mathbf{r},t)$
$\displaystyle=\sum_{i}a_{i}(t)\mathbf{E}_{i}(\mathbf{r})$ (4)
$\displaystyle\mathbf{H}(\mathbf{r},t)$
$\displaystyle=\sum_{i}b_{i}(t)\mathbf{H}_{i}(\mathbf{r})$ (5)

where $\mathbf{E}_{i}(\mathbf{r})$ ($\mathbf{H}_{i}(\mathbf{r})$) is the
electric (magnetic) field of the mode associated to the resonator at site $i$.
The two fields verify

$\nabla\times\mathbf{E}_{i}=\mu_{0}\omega_{i}H_{i}$ (6)

where $\omega_{i}$ is the resonance frequency of the mode of the resonator at
site $i$. In addition to the overlap integrals defined in (3), we consider the
integral

$M_{ij}=\int_{S}[\mathbf{E}_{i}(\mathbf{r})\times\mathbf{H}_{j}(\mathbf{r})]\cdot
d\mathbf{S}$ (7)

where the surface $S$ corresponds to all the metallic boundaries in the
circuit, which consists of the resonators and the walls of the box enclosing
the sample. This integral can be rewritten in terms of the overlap integrals
$D_{ij}$ and $G_{ij}$ assuming that

$M_{ij}\approx\int_{S_{i}\cup
S_{j}}[\mathbf{E}_{i}(\mathbf{r})\times\mathbf{H}_{j}(\mathbf{r})]\cdot
d\mathbf{S}$ (8)

where $S_{i}$ is the surface of the resonator at site $i$. The integral over
$S_{i}$ is null because $\mathbf{E}_{i}(\mathbf{r})$ is normal to the surface
of the resonator. Using that the current density in the resonator at site $j$
is given by
$\mathbf{J}_{j}(\mathbf{r})=\mathbf{n}\times\mathbf{H}_{j}(\mathbf{r})$, where
$\mathbf{n}$ is the normal to the sample surface, we obtain

$M_{ij}\approx\int\mathbf{E}_{i}(\mathbf{r})\cdot\mathbf{J}_{j}(\mathbf{r})dV$
(9)

Using this approximation and considering the volume integral of
$\nabla\cdot(\mathbf{E}_{i}\times\mathbf{H}_{j})$, one obtains

$2M_{ij}=\omega_{i}G_{ij}-\omega_{j}D_{ij}$ (10)

In matrix form, the last equation writes

$2M=\Omega G-D\Omega$ (11)

where $\Omega$ is the diagonal matrix with elements $\Omega_{ii}=\omega_{i}$.
In the same way, the volume integrals of
$\nabla\cdot(\mathbf{H}_{i}\times\mathbf{E})$ and
$\nabla\cdot(\mathbf{E}_{i}\times\mathbf{H})$ lead to the following two
equations

$\displaystyle G\dot{b}$ $\displaystyle=-\Omega Da$ (12) $\displaystyle 2Mb$
$\displaystyle=\Omega Gb-\omega D\dot{a}$ (13)

Using the approximation (11) for $M$ and eliminating $b$, we look for periodic
solutions at frequency $\omega$ and obtain the following eigenvalue problem

$\omega^{2}a=\Omega G^{-1}\Omega Da$ (14)

In the case of two resonators only, we recover the result given in (1). The
matrices $G$, $\Omega$ and $D$ are real symmetric matrices, but, in general,
$\Omega G^{-1}\Omega D$ is not symmetric. But its eigenvalues are real and
positive with non-orthogonal eigenvectors.

Equation (14) can be used to find the coupled mode frequencies of any ensemble
of coupled resonators. In order to take advantage of the lattice periodicity,
we relabel the $a_{i}$ and $b_{i}$ amplitudes in a given lattice cell as
$a_{\mu}(\mathbf{R})$ and $b_{\mu}(\mathbf{R})$, where $\mathbf{R}$ is Bravais
lattice vector identifying the position of the cell in the lattice and the
$\mu$ index distinguishes the A and B sites. Using the same notation, we
define the following overlap integrals between neighbouring sites in the
lattice

$\displaystyle\langle\mathbf{E}_{\mu}\mathbf{E}_{\nu}\rangle_{\mathbf{R}}$
$\displaystyle=\int\epsilon(\mathbf{r}^{\prime})\,\mathbf{E}_{\mu}(r^{\prime})\cdot\mathbf{E}_{\nu}(\mathbf{r}^{\prime}-\mathbf{R})\
d^{3}\mathbf{r}^{\prime}$ (15)
$\displaystyle\langle\mathbf{H}_{\mu}\mathbf{H}_{\nu}\rangle_{\mathbf{R}}$
$\displaystyle=\int\mu_{0}\,\mathbf{H}_{\mu}(\mathbf{r}^{\prime})\cdot\mathbf{H}_{\nu}(\mathbf{r}^{\prime}-\mathbf{R})\
d^{3}\mathbf{r}^{\prime}$ (16)

We now look for periodic solutions over the lattice as

$a_{\mu}(R)=e^{i\mathbf{k}\cdot\mathbf{R}}a_{\mu}(\mathbf{k})\ \
b_{\mu}(R)=e^{i\mathbf{k}\cdot\mathbf{R}}b_{\mu}(\mathbf{k})$ (17)

and define the Fourier transform of the $D$ and $G$ overlap matrices

$\displaystyle D_{\mu\nu}(\mathbf{k})$
$\displaystyle=\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}\langle\mathbf{E}_{\mu}\mathbf{E}_{\nu}\rangle_{\mathbf{R}}$
(18) $\displaystyle G_{\mu\nu}(\mathbf{k})$
$\displaystyle=\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}\langle\mathbf{H}_{\mu}\mathbf{H}_{\nu}\rangle_{\mathbf{R}}$
(19)

We suppose that the lattice contains $N$ cells and we look for solutions with
periodic boundary conditions such that
$\sum_{\mathbf{R}}e^{j\mathbf{k}\cdot\mathbf{R}}=N\delta_{\mathbf{k},0}$,
where $\mathbf{k}$ can take $N$ values in the first Brillouin zone. Equations
(12) and (13) then lead to

$\Omega G^{-1}(\mathbf{k})\Omega
D(\mathbf{k})a(\mathbf{k})=\omega^{2}(\mathbf{k})a(\mathbf{k})$ (20)

where $\Omega$ is now the diagonal matrix with the resonance frequencies of
the A and B sites. The resulting eigenvalue problem (20) gives the dispersion
of the two bands as a function of $\mathbf{k}$.

We obtain the band structure of our lattices by including all couplings up to
the third neighbour coupling (see figure 4) in equations (18) and (19) and
solving (20). The results are shown in figure . In order to compare the
obtained band structure with the one of a tight-binding model with the same
range of coupling, we consider the case of the G lattice where the A and B
sites have the same resonance frequency. In order to obtain a simple
analytical formula, we neglect the dependence of the overlap integrals with
the direction of separation (e.g. we suppose that
$\langle\mathbf{E}_{A}\mathbf{E}_{B}\rangle_{\mathbf{a}_{1}}=\langle\mathbf{E}_{A}\mathbf{E}_{B}\rangle_{\mathbf{a}}$).
With this approximation, a straightforward calculation shows that
$D(\mathbf{k})$ and $G(\mathbf{k})$ are diagonal in the same basis, leading to
the following dispersion relation $\omega_{\pm}(\mathbf{k})$

$\omega_{\pm}(\mathbf{k})=\sqrt{\frac{1+\kappa_{e}^{(2)}\,f_{2}(\mathbf{k})\pm|\kappa_{e}^{(1)}\,f_{1}(\mathbf{k})+\kappa_{e}^{(3)}\,f_{3}(\mathbf{k})|}{1+\kappa_{m}^{(2)}\,f_{2}(\mathbf{k})\pm|\kappa_{m}^{(1)}\,f_{1}(\mathbf{k})+\kappa_{m}^{(3)}\,f_{3}(\mathbf{k})|}}\
\omega_{0}$ (21)

where the electric coupling constants are given by

$\displaystyle\kappa^{(1)}_{e}$
$\displaystyle=\langle\mathbf{E}_{A}\mathbf{E}_{B}\rangle_{\mathbf{0}}/\langle\mathbf{E}_{A}\mathbf{E}_{A}\rangle_{\mathbf{0}}$
$\displaystyle\kappa^{(2)}_{e}$
$\displaystyle=\langle\mathbf{E}_{A}\mathbf{E}_{A}\rangle_{\mathbf{a}_{1}}/\langle\mathbf{E}_{A}\mathbf{E}_{A}\rangle_{\mathbf{0}}$
$\displaystyle\kappa^{(3)}_{e}$
$\displaystyle=\langle\mathbf{E}_{A}\mathbf{E}_{A}\rangle_{\mathbf{a}_{1}+\mathbf{a}_{2}}/\langle\mathbf{E}_{A}\mathbf{E}_{A}\rangle_{\mathbf{0}}$

The magnetic couplings $\kappa_{m}^{(i)}$ are given by the same expressions in
terms of magnetic field overlaps and the link functions $f_{i}(\mathbf{k})$
correspond to

$\displaystyle f_{1}(\mathbf{k})$
$\displaystyle=1+e^{i\mathbf{k}\cdot\mathbf{a}_{1}}+e^{i\mathbf{k}\cdot\mathbf{a}_{2}}$
$\displaystyle f_{2}(\mathbf{k})$
$\displaystyle=2[\,\cos\mathbf{k}\cdot\mathbf{a}_{1}+\cos\mathbf{k}\cdot\mathbf{a}_{2}+\cos\mathbf{k}\cdot(\mathbf{a}_{1}-\mathbf{a}_{2})\,]$
$\displaystyle f_{3}(\mathbf{k})$
$\displaystyle=e^{i\mathbf{k}\cdot(\mathbf{a}_{1}+\mathbf{a}_{2})}+2\cos\mathbf{k}\cdot(\mathbf{a}_{1}-\mathbf{a}_{2})$

Expression (21) can be compared to the tight-binding dispersion relation

$\omega^{\rm
TB}_{\pm}(\mathbf{k})=\omega_{0}+t_{2}\,f_{2}(\mathbf{k})\pm|t_{1}\,f_{1}(\mathbf{k})+t_{3}\,f_{3}(\mathbf{k})|$
(22)

which comes from the diagonalisation of the tight-binding Hamiltonian

$H^{\rm
TB}(\mathbf{k})=\begin{pmatrix}\omega_{0}+t_{2}f_{2}(\mathbf{k})&t_{1}f_{1}(\mathbf{k})+t_{3}f_{3}(\mathbf{k})\\\
t_{1}f_{1}^{*}(\mathbf{k})+t_{3}f_{3}^{*}(\mathbf{k})&\omega_{0}+t_{2}f_{2}(\mathbf{k})\end{pmatrix}$
(23)

The coupling $t_{i}$ corresponds to the $i$-th nearest neighbour coupling. The
dispersion relation (21) and (22) match if one keeps only nearest neighbour
terms, performs a first order expansion in $\kappa_{e}^{(1)}$ and
$\kappa_{m}^{(1)}$ and identifies $t_{1}$ to
$\omega_{0}(\kappa_{e,1}-\kappa_{m,1})/2$. Including 2nd and 3rd nearest
neighbour overlap integrals lead to a dispersion relation that cannot be
identified to its tight-binding counterpart with the same coupling range.

As mentioned in the introduction, a strong motivation to build lattices of
superconducting resonators is to include non-linear elements in the lattice.
For example, in the case considered here, one could design spirals with an
empty area in the center where a transmon might be inserted and coupled to the
spiral. In this context, it can be interesting to obtain an effective tight-
binding Hamiltonian that describes the linear behaviour of the lattice. Using
the CMT approach, we obtained the two equations of motion

$\displaystyle\dot{a}$ $\displaystyle=\Omega b$ (24) $\displaystyle G\dot{b}$
$\displaystyle=-\Omega Da$ (25)

And the total energy in the system is given by

$H=\frac{1}{2}a^{T}D\,a+\frac{1}{2}b^{T}G\,b$ (26)

In order to identify $a$ and $b$ with the two quadratures of the resonator
modes, we normalize the basis fields to

$D_{ii}=G_{ii}=\hbar\omega_{i}$ (27)

With this definition, and in the absence of coupling between the resonators
($G=D=\hbar\Omega$), the equations of motion coincide with the Hamilton
equations derived from $H$ assuming standard commutation rules between the
mode quadratures, $[a_{m},a_{n}]=0$, $[b_{m},b_{n}]=0$ and
$[a_{m},b_{n}]=i\hbar\delta_{mn}$. In the following, we drop the $\hbar$
factors and set $\hbar=1$. When the overlap between adjacent sites are non-
zero, the equations of motion do not coincide anymore with the naive Hamilton
equations considering that the $a$’s and $b$’s commute. This comes from the
fact that the basis used to project the Maxwell equations is not orthonormal,
which modifies the commutation relations. This problem is well known in the
calculation of electronic band structure using the LCAO method. In order to
find an orthonormal basis, we adapt the Löwdin procedure used for electrons
Aiken  _et al._ (1980). We first apply the canonical transformation
$a\rightarrow\Omega^{1/2}a$ and $b\rightarrow\Omega^{-1/2}b$ to remove the
$\Omega$ dependence in the equation of motion, which transforms the
Hamiltonian to

$H=\frac{1}{2}a^{T}\tilde{D}\,a+\frac{1}{2}b^{T}\tilde{G}\,b$ (28)

with $\tilde{D}=\Omega^{1/2}D\Omega^{1/2}$ and
$\tilde{G}=\Omega^{-1/2}G\Omega^{-1/2}$. This step is not necessary if
$\Omega$ is proportional to the identity. The equations of motion become
$\dot{a}=b$ and $G\dot{b}=-Da$. We now apply the Löwdin transformation
$a\rightarrow\tilde{G}^{-1/2}a$ and $b\rightarrow\tilde{G}^{-1/2}b$ in order
to restore canonical commutation relations. The Hamiltonian becomes

$H=\frac{1}{2}a^{T}\tilde{G}^{-1/2}\tilde{D}\tilde{G}^{-1/2}\,a+\frac{1}{2}b^{T}\,b$
(29)

The equations of motion are modified to $\dot{a}=b$ and
$\dot{b}=-\tilde{G}^{-1/2}\tilde{D}\tilde{G}^{-1/2}a$, which now correspond to
the Hamilton equations obtained with $H$ assuming canonical commutation
relations. We finally apply the reverse canonical transformation
$a\rightarrow\Omega^{-1/2}a$ and $b\rightarrow\Omega^{1/2}b$ to explicitly
restore the $\Omega$ dependence. The final Hamiltonian is

$H=\frac{1}{2}a^{T}\Omega^{-1/2}\tilde{G}^{-1/2}\tilde{D}\tilde{G}^{-1/2}\Omega^{-1/2}\,a+\frac{1}{2}b^{T}\Omega\,b$
(30)

We can now introduce the bosonic operators $c=(a+ib)/\sqrt{2}$ and rewrite $H$
as

$H=\sum_{i}\frac{\omega_{i}}{2}(c_{i}^{\dagger}c_{i}+c_{i}c^{\dagger}_{i})+\frac{1}{2}(c+c^{\dagger})^{T}J(c+c^{\dagger})$
(31)

with
$2J=\Omega^{-1/2}\tilde{G}^{-1/2}\tilde{D}\tilde{G}^{-1/2}\Omega^{-1/2}-\Omega$.
A rotating wave approximation can then be performed to transform the last term
to $c^{T}(J/2)c^{\dagger}+{\rm hc}$.

###  III.2 Admittance matrix

In order to model more accurately the G lattice, we extend the admittance
matrix method introduced in section I to an infinite lattice of resonators and
we also take into account the possibility to have $N\geq 1$ microwave ports
per resonator. We introduce the voltage amplitude $V(\mathbf{r}^{\prime})$
that is a complex vector with $2N$ components corresponding to the voltage at
the $2N$ ports describing the voltage in the lattice cell at position
$\mathbf{r}^{\prime}$. At a given frequency $\omega$, the Kirchhoff’s circuit
laws describing the lattice can be written

$\sum_{\mathbf{r}}Y_{\mathbf{r}}(\omega)V(\mathbf{r}^{\prime}+\mathbf{r})=0$
(32)

where the sum over $\mathbf{r}$ is over the points of the Bravais lattice. The
matrix $Y_{0}(\omega)$ is the admittance matrix between the ports belonging to
the same cell, while $Y_{\mathbf{r}}(\omega)$ is the admittance matrix between
the ports corresponding to two cells separated by $\mathbf{r}$, it has the
property $Y_{-\mathbf{r}}(\omega)=Y_{\mathbf{r}}^{T}(\omega)$. We look for a
periodic solution $V(\mathbf{r})=V(\mathbf{k})e^{i\mathbf{k}\mathbf{r}}$ and
obtain

$\sum_{\mathbf{r}}Y_{\mathbf{r}}(\omega)e^{i\mathbf{k}\mathbf{r}}V(\mathbf{k})=0$
(33)

where $V(\mathbf{k})$ is a $2N$ element vector. The dispersion relation is
then obtained by solving

$\det\left(\sum_{r}Y_{\mathbf{r}}(\omega)e^{i\mathbf{k}\mathbf{r}}\right)=0$
(34)

The $Y_{\mathbf{r}}(\omega)$ matrix is obtained from the numerical simulation
of a finite size lattice. The CMT calculations indicate that one should
simulate a lattice with enough sites such that the central site is surrounded
by all neighbouring sites up to the 3rd nearest neighbour. We therefore choose
the geometry shown in figure with 16 sites. In order for the sum over
$\mathbf{r}$ to converge rapidly when the distance between sites increases,
the number of ports $N$ and their location on the spiral must be well chosen.
We observe from the simulations that $N=1$ is not sufficient, while results
obtained for $N>2$ and various port locations give the same results. The
results of the final simulations shown in are obtained with $N=3$ ports
located as shown in .

Compared to the CMT model, the admittance matrix is numerically more demanding
because it requires the simulation of a much larger (here 16) number of sites.
The dispersion obtained for the S lattice confirms the prediction of the CMT
model, which were already in good agreement with our experimental data. For
the G lattice, we observe a deviation compared to the CMT model leading to a
much better agreement with the measured dispersion. We can thus conclude that
our lattice are well under control and that _ab initio_ electromagnetic
simulation accurately describe their properties.

###### Acknowledgements.

The authors gratefully acknowledge the conversations and insights of Marco
Aprili.

## References

  * Koch and Le Hur (2009) J. Koch and K. Le Hur, Physical Review A 80, 023811 (2009).
  * Koch  _et al._ (2010) J. Koch, A. A. Houck, K. L. Hur, and S. M. Girvin, Physical Review A - Atomic, Molecular, and Optical Physics 82, 1 (2010), arXiv:1006.0762 .
  * Schmidt and Koch (2013) S. Schmidt and J. Koch, Annalen der Physik 525, 395 (2013).
  * Roushan  _et al._ (2017) 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, E. Kapit, H. Neven, and J. Martinis, Nature Physics 13, 146 (2017).
  * Fitzpatrick  _et al._ (2017) M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Physical Review X 7, 011016 (2017).
  * Ma  _et al._ (2018) R. Ma, B. Saxberg, C. Owens, N. Leung, Y. Lu, J. Simon, and D. I. Schuster, arXiv:1807.11342 [cond-mat, physics:quant-ph] (2018), arXiv:1807.11342 [cond-mat, physics:quant-ph] .
  * Guo  _et al._ (2020) Q. Guo, C. Cheng, H. Li, S. Xu, P. Zhang, Z. Wang, C. Song, W. Liu, W. Ren, H. Dong, R. Mondaini, and H. Wang, “Stark many-body localization on a superconducting quantum processor,” (2020), arXiv:2011.13895 [quant-ph] .
  * Fedorov  _et al._ (2020) G. P. Fedorov, S. V. Remizov, D. S. Shapiro, W. V. Pogosov, E. Egorova, I. Tsitsilin, M. Andronik, A. A. Dobronosova, I. A. Rodionov, O. V. Astafiev, and A. V. Ustinov, “Photon transport in a bose-hubbard chain of superconducting artificial atoms,” (2020), arXiv:2011.11454 [quant-ph] .
  * Underwood  _et al._ (2016) D. L. Underwood, W. E. Shanks, A. C. Y. Li, L. Ateshian, J. Koch, and A. A. Houck, Physical Review X 6, 021044 (2016).
  * Kollár  _et al._ (2019) A. J. Kollár, M. Fitzpatrick, and A. A. Houck, Nature 571, 45 (2019).
  * Lu  _et al._ (2014) L. Lu, J. D. Joannopoulos, and M. Soljačić, Nature Photonics 8, 821 (2014).
  * Khanikaev and Shvets (2017) A. B. Khanikaev and G. Shvets, Nature Photonics 11, 763 (2017).
  * Ozawa  _et al._ (2018) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, arXiv:1802.04173 [cond-mat, physics:physics, physics:quant-ph] (2018), arXiv:1802.04173 [cond-mat, physics:physics, physics:quant-ph] .
  * Culbertson  _et al._ (1998) J. C. Culbertson, H. S. Newman, and C. Wilker, Journal of Applied Physics 84, 2768 (1998).
  * Zhuravel  _et al._ (2012) A. P. Zhuravel, C. Kurter, A. V. Ustinov, and S. M. Anlage, Physical Review B 85, 134535 (2012).
  * Averkin  _et al._ (2013) A. S. Averkin, A. P. Zhuravel, A. Karpov, S. M. Anlage, and A. V. Ustinov, in  _2013 7th International Congress on Advanced Electromagnetic Materials in Microwaves and Optics_ (IEEE, Talence, France, 2013) pp. 142–144.
  * Wang  _et al._ (2019) H. Wang, A. Zhuravel, S. Indrajeet, B. Taketani, M. Hutchings, Y. Hao, F. Rouxinol, F. Wilhelm, M. LaHaye, A. Ustinov, and B. Plourde, Phys. Rev. Applied 11, 054062 (2019).
  * Huang (1994) W.-P. Huang, JOSA A 11, 963 (1994), publisher: Optical Society of America.
  * Powell  _et al._ (2010) D. A. Powell, M. Lapine, M. V. Gorkunov, I. V. Shadrivov, and Y. S. Kivshar, Physical Review B 82, 155128 (2010).
  * Lomanets  _et al._ (2012) V. Lomanets, O. Zhuromskyy, G. Onishchukov, and U. Peschel, Physical Review B 85, 125110 (2012).
  * Elnaggar  _et al._ (2015) S. Y. Elnaggar, R. Tervo, and S. M. Mattar, IEEE Transactions on Microwave Theory and Techniques 63, 2115 (2015), arXiv:1305.6085 .
  * Semenoff (1984) G. W. Semenoff, Physical Review Letters 53, 2449 (1984).
  * Morvan  _et al._ (2020) A. Morvan, M. Féchant, G. Aiello, J. Gabelli, and J. Estève, to be published (2020).
  * Zhuravel  _et al._ (2006) A. P. Zhuravel, A. G. Sivakov, O. G. Turutanov, A. N. Omelyanchouk, S. M. Anlage, A. Lukashenko, A. V. Ustinov, and D. Abraimov, Low Temperature Physics 32, 592 (2006).
  * Aiken  _et al._ (1980) J. Aiken, J. A. Erdos, and J. Goldstein, International Journal of Quantum Chemistry 18, 1101 (1980).


${\mbox{\boldmath$A$}}_{M+1,\,M}^{(\mathcal{I}_{j})}$ contains nonzero
elements in sub-matrix
${\mbox{\boldmath$A$}}_{\textbf{2},\,\textbf{1}}^{(\mathcal{I}_{j})}$, and all
other blocks in
${\mbox{\boldmath$A$}}_{\textbf{2},\,\textbf{1}}^{(\mathcal{I}_{j})}$ contain
only zero entries. From Cauchy–Schwarz inequality, we have:

$\left\|{\mbox{\boldmath$R$}}_{\textbf{2},\textbf{1}}^{(\mathcal{I}_{j})}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}\right\|_{1}=\left\|\bigl{[}\textrm{diag}(\mathbbm{1}^{T}{\mbox{\boldmath$A$}}_{M+1,M}^{(\mathcal{I}_{j})})\bigr{]}\bigl{[}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}(\mathcal{G}_{M})\bigr{]}\right\|_{1}\leq\left\|\textrm{diag}(\mathbbm{1}^{T}{\mbox{\boldmath$A$}}_{M+1,M}^{(\mathcal{I}_{j})})\right\|_{1}\cdot\left\|{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}(\mathcal{G}_{M})\right\|_{1},$

where
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}(\mathcal{G}_{M})$ is
the sub-vector corresponding to the state partition $\mathcal{G}_{M}$.
Furthermore, according to Lemma 15 and Eqn. (13) after replacing the subscript
$i$ in $\tilde{\pi}_{i}^{(\infty)}$ with $M$ and taking into consideration of
the equivalence of the infinite space $\Omega^{(\mathcal{I}_{j})}$ and
$\Omega^{(\infty)}$ in regard to truncation at $\mathcal{I}_{i}$, we have the
probability of the boundary block $\mathcal{G}_{M}$:
$\left\|{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}(\mathcal{G}_{M})\right\|_{1}\rightarrow
0$ when $M\rightarrow\infty$. When synthesis reactions are concentration
independent (zero order reactions) as usually the case [57], the norm
$\left\|\textrm{diag}(\mathbbm{1}^{T}{\mbox{\boldmath$A$}}_{M+1,M}^{(\mathcal{I}_{j})})\right\|_{1}$
is a constant representing the total synthesis rates over states in
$\mathcal{G}_{M}$. We have:
$\left\|{\mbox{\boldmath$R$}}_{\textbf{2},\textbf{1}}^{(\mathcal{I}_{j})}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}\right\|_{1}\rightarrow
0$ when $M\rightarrow\infty$. Therefore with Eqn. (59), we have:

$\lim_{M\rightarrow\infty}\left\|\tilde{{\mbox{\boldmath$A$}}}^{(\mathcal{I}_{i,j})}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}\right\|_{1}=0.$

Hence,

$\lim_{M\rightarrow\infty}\tilde{{\mbox{\boldmath$A$}}}^{(\mathcal{I}_{i,j})}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}=\textbf{0}.$

That is, when the maximum copy number limit of the $i$-th MEG is sufficiently
large, both ${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}$ and
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}$ are the steady
state solutions of
$\tilde{{\mbox{\boldmath$A$}}}^{(\mathcal{I}_{i,j})}{\mbox{\boldmath$y$}}=0$.
According to Perron–Frobenius theorem for the transition rate matrix of
continuous-time Markov chains [58], the dCME governed by
$\tilde{{\mbox{\boldmath$A$}}}^{(\mathcal{I}_{i,j})}$ has a globally unique
stationary distribution. In addition, by construction of the matrix, via
enumeration of the state space, matrix
$\tilde{{\mbox{\boldmath$A$}}}^{(\mathcal{I}_{i,j})}$ is irreducible, as all
microstates in the state space can be reached from the initial state.
Therefore, the matrix $\tilde{{\mbox{\boldmath$A$}}}^{(\mathcal{I}_{i,j})}$
has only one zero eigenvalue [58], both
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}$ and
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}$ are eigenvectors
corresponding to the eigenvalue $0$. Therefore, we have the relationship
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}=c{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}$,
where $c$ is an arbitrary real number. As both vectors are non-negative, and
$\mathbbm{1}^{T}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}\leq
1=\mathbbm{1}^{T}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}$, there must
exist an
$\epsilon=1-\mathbbm{1}^{T}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}\geq
0$, such that
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}=(1+\epsilon){\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}$.
According to Lemma 15,
$\mathbbm{1}^{T}{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}\rightarrow
1$, when the maximum copy number limit of the $i$-th MEG goes to infinity.
Therefore we have $\epsilon\rightarrow 0$ when $M\rightarrow\infty$.
Therefore, we have shown both
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}\geq{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}$
and
${\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{i,j})}\rightarrow{\mbox{\boldmath$\pi$}}^{(\mathcal{I}_{j})}_{\textbf{1}}$
component-wise, when the maximum copy number limit of the $i$-th MEG goes to
infinity.

## References

  * 1. Jacob Stewart-Ornstein and Hana El-Samad.  Stochastic modeling of cellular networks.  Computational Methods in Cell Biology, 110:111, 2012.
  * 2. Hong Qian.  Cooperativity in cellular biochemical processes: noise-enhanced sensitivity, fluctuating enzyme, bistability with nonlinear feedback, and other mechanisms for sigmoidal responses.  Annual Review of Biophysics, 41:179–204, 2012.
  * 3. H.H. McAdams and A. Arkin.  It’s a noisy business! Genetic regulation at the nanomolar scale.  Trends in Genetics, 15(2):65–69, 1999.
  * 4. Darren J Wilkinson.  Stochastic modelling for quantitative description of heterogeneous biological systems.  Nature Reviews Genetics, 10(2):122–133, 2009.
  * 5. Youfang Cao, Hsiao-Mei Lu, and Jie Liang.  Probability landscape of heritable and robust epigenetic state of lysogeny in phage lambda.  Proceedings of the National Academy of Sciences of the United States of America, 107(43):18445–18450, 2010.
  * 6. D. T. Gillespie.  Exact stochastic simulation of coupled chemical reactions.  Journal of Physical Chemistry, 81:2340–2361, 1977.
  * 7. Daniel T. Gillespie.  A rigorous derivation of the chemical master equation.  Physica A, 188:404–425, 1992.
  * 8. N.G. Van Kampen.  Stochastic processes in physics and chemistry, 3rd Edition.  Elsevier Science and Technology books, 2007.
  * 9. D.A. Beard and H. Qian.  Chemical biophysics: quantitative analysis of cellular systems.  Cambridge Univ Pr, 2008.
  * 10. Daniel T. Gillespie.  A diffusional bimolecular propensity function.  The Journal of Chemical Physics, 131(16):164109, 2009.
  * 11. I.G. Darvey, B.W. Ninham, and P.J. Staff.  Stochastic models for second order chemical reaction kinetics. the equilibrium state.  The Journal of Chemical Physics, 45:2145–2155, 1966.
  * 12. D.A. McQuarrie.  Stochastic approach to chemical kinetics.  Journal of Applied Probability, 4:413–478, 1967.
  * 13. H.M. Taylor and S. Karlin.  An Introduction to Stochastic Modeling, 3rd Ed. Academic Press, 1998.
  * 14. I.J. Laurenzi.  An analytical solution of the stochastic master equation for reversible bimolecular reaction kinetics.  The Journal of Chemical Physics, 113:3315–3322, 2000.
  * 15. Melissa Vellela and Hong Qian.  A quasistationary analysis of a stochastic chemical reaction: Keizer’s paradox.  Bulletin of Mathematical Biology, 69(5):1727–1746, 2007.
  * 16. Daniel. T. Gillespie.  The chemical langevin equation.  The Journal of Chemical Physics, 113:297–306, 2000.
  * 17. N. G. Van Kampen.  A power series expansion of the master equation.  Canadian Journal of Physics, 39(4):551–567, 1961.
  * 18. Daniel T. Gillespie.  The chemical Langevin and Fokker−Planck equations for the reversible isomerization reaction†.  The Journal of Physical Chemistry A, 106(20):5063–5071, 2002.
  * 19. Eric L. Haseltine and James B. Rawlings.  Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics.  The Journal of Chemical Physics, 117(15):6959–6969, 2002.
  * 20. C. W. Gardiner.  Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences.  Springer, New York, 2004.
  * 21. R. Grima, P. Thomas, and A. V. Straube.  How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations?  The Journal of Chemical Physics, 135(8):084103, Aug 2011.
  * 22. P. Thomas, H. Matuschek, and R. Grima.  How reliable is the linear noise approximation of gene regulatory networks?  BMC Genomics, 14 Suppl 4:S5, 2013.
  * 23. Brian Munsky and Mustafa Khammash.  The finite state projection algorithm for the solution of the chemical master equation.  The Journal of Chemical Physics, 124(4):044104, 2006.
  * 24. Youfang Cao and Jie Liang.  Optimal enumeration of state space of finitely buffered stochastic molecular networks and exact computation of steady state landscape probability.  BMC Systems Biology, 2(1):30, 2008.
  * 25. Shev MacNamara, Alberto M Bersani, Kevin Burrage, and Roger B Sidje.  Stochastic chemical kinetics and the total quasi-steady-state assumption: application to the stochastic simulation algorithm and chemical master equation.  The Journal of chemical physics, 129(9):095105, 2008.
  * 26. Shev MacNamara, Kevin Burrage, and Roger B Sidje.  Multiscale modeling of chemical kinetics via the master equation.  Multiscale Modeling & Simulation, 6(4):1146–1168, 2008.
  * 27. Verena Wolf, Rushil Goel, Maria Mateescu, and Thomas Henzinger.  Solving the chemical master equation using sliding windows.  BMC Systems Biology, 4(1):42, 2010.
  * 28. Tobias Jahnke.  On reduced models for the chemical master equation.  Multiscale Modeling & Simulation, 9(4):1646–1676, 2011.
  * 29. Roger B Sidje.  Expokit: a software package for computing matrix exponentials.  ACM Transactions on Mathematical Software (TOMS), 24(1):130–156, 1998.
  * 30. Brian Munsky and Mustafa Khammash.  A multiple time interval finite state projection algorithm for the solution to the chemical master equation.  Journal of Computational Physics, 226(1):818 – 835, 2007.
  * 31. Jianjun Paul Tian and D Kannan.  Lumpability and commutativity of Markov processes.  Stochastic analysis and Applications, 24(3):685–702, 2006.
  * 32. Laurent Truffet.  Near complete decomposability: bounding the error by a stochastic comparison method.  Advances in Applied Probability, pages 830–855, 1997.
  * 33. Peter Buchholz.  Exact and ordinary lumpability in finite Markov chains.  Journal of Applied Probability, pages 59–75, 1994.
  * 34. W.J. Stewart.  Introduction to the numerical solution of Markov chains.  Princeton University Press NJ, 1994.
  * 35. Hendrik Vantilborgh.  Aggregation with an error of o$(\epsilon^{2})$.  Journal of the ACM (JACM), 32(1):162–190, 1985.
  * 36. John G Kemeny and James Laurie Snell.  Finite Markov chains, volume 210.  Springer-Verlag New York, 1976.
  * 37. A Irle.  Stochastic ordering for continuous-time processes.  Journal of Applied Probability, pages 361–375, 2003.
  * 38. B.J. Daigle, M.K. Roh, D.T. Gillespie, and L.R. Petzold.  Automated estimation of rare event probabilities in biochemical systems.  The Journal of Chemical Physics, 134:044110, 2011.
  * 39. M.K. Roh, B.J. Daigle, D.T. Gillespie, and L.R. Petzold.  State-dependent doubly weighted stochastic simulation algorithm for automatic characterization of stochastic biochemical rare events.  Journal of Chemical Physics, 135(23):234108, 2011.
  * 40. Youfang Cao and Jie Liang.  Adaptively biased sequential importance sampling for rare events in reaction networks with comparison to exact solutions from finite buffer dCME method.  The Journal of Chemical Physics, 139(2):025101, 2013.
  * 41. Yuichi Taniguchi, Paul J Choi, Gene-Wei Li, Huiyi Chen, Mohan Babu, Jeremy Hearn, Andrew Emili, and X Sunney Xie.  Quantifying e. coli proteome and transcriptome with single-molecule sensitivity in single cells.  Science, 329(5991):533–538, 2010.
  * 42. Timothy S Gardner, Charles R Cantor, and James J Collins.  Construction of a genetic toggle switch in Escherichia coli.  Nature, 403(6767):339–342, 2000.
  * 43. Thomas B Kepler and Timothy C Elston.  Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations.  Biophysical Journal, 81(6):3116–3136, 2001.
  * 44. Keun-Young Kim and Jin Wang.  Potential energy landscape and robustness of a gene regulatory network: toggle switch.  PLoS Computational Biology, 3(3):e60, 2007.
  * 45. D. Schultz, J. N. Onuchic, and P. G. Wolynes.  Understanding stochastic simulations of the smallest genetic networks.  The Journal of Chemical Physics, 126(24):245102, 2007.
  * 46. Brian Munsky and Mustafa Khammash.  The finite state projection approach for the analysis of stochastic noise in gene networks.  Automatic Control, IEEE Transactions on, 53(Special Issue):201–214, 2008.
  * 47. Peter Deuflhard, Wilhelm Huisinga, T Jahnke, and Michael Wulkow.  Adaptive discrete Galerkin methods applied to the chemical master equation.  SIAM Journal on Scientific Computing, 30(6):2990–3011, 2008.
  * 48. Paul Sjöberg, Per Lötstedt, and Johan Elf.  Fokker–Planck approximation of the master equation in molecular biology.  Computing and Visualization in Science, 12(1):37–50, 2009.
  * 49. Vladimir Kazeev, Mustafa Khammash, Michael Nip, and Christoph Schwab.  Direct solution of the chemical master equation using quantized tensor trains.  PLoS Computational Biology, 10(3):e1003359, 03 2014.
  * 50. Adam Arkin, John Ross, and Harley H. McAdams.  Stochastic kinetic analysis of developmental pathway bifurcation in phage $\lambda$-infected Escherichia coli cells.  Genetics, 149(4):1633–1648, 1998.
  * 51. Erik Aurell, Stanley Brown, Johan Johanson, and Kim Sneppen.  Stability puzzles in phage $\lambda$.  Physical Review E, 65(5):051914, 2002.
  * 52. Erik Aurell and Kim Sneppen.  Epigenetics as a first exit problem.  Physical Review Letters, 88(4):048101, 2002.
  * 53. X.-M. Zhu, L. Yin, L. Hood, and P. Ao.  Robustness, stability and efficiency of phage lambda genetic switch: dynamical structure analysis.  Journal of Bioinformatics and Computational Biology, 2:785–817, 2004.
  * 54. X.-M. Zhu, L. Yin, L. Hood, and P. Ao.  Calculating biological behaviors of epigenetic states in the phage $\lambda$ life cycle.  Functional & Integrative Genomics, 4(3):188–195, 2004.
  * 55. Shuohao Liao, T Vejchodsky, and Radek Erban.  Tensor methods for parameter estimation and bifurcation analysis of stochastic reaction networks.  Journal of the Royal Society Interface, 12(108):20150233, 2015.
  * 56. F. Verstraete and J. I. Cirac.  Matrix product states represent ground states faithfully.  Phys. Rev. B, 73:094423, Mar 2006.
  * 57. Philip Nelson.  Physical Models of Living Systems.  Macmillan, 2015.
  * 58. Carl D Meyer.  Matrix analysis and applied linear algebra.  SIAM, 2000.
  * 59. M. Li, W. McClure, and M. Susskind.  Changing the mechanism of transcriptional activation by phage lambda repressor.  Proceedings of the National Academy of Sciences of the United States of America, 94(8):3691–3696, 1997.
  * 60. D. Hawley and W. McClure.  In vitro comparison of initiation properties of bacteriophage lambda wild-type PR and x3 mutant promoters.  Proceedings of the National Academy of Sciences of the United States of America, 77(11):6381–6385, 1980.
  * 61. D. Hawley and W. McClure.  Mechanism of activation of transcription initiation from the lambda PRM promoter.  Journal of Molecular Biology, 157(3):493–525, 1982.
  * 62. Madeline A. Shea and Gary K. Ackers.  The $OR$ control system of bacteriophage lambda a physical-chemical model for gene regulation.  Journal of Molecular Biology, 181(2):211–230, 1985.
  * 63. Céline Kuttler and Joachim Niehren.  Gene Regulation in the Pi Calculus: Simulating Cooperativity at the Lambda Switch.  Transactions on Computational Systems Biology VII, 4230:24–55, 2006\.

## Figure Legends

Figure 1:  An illustration for the boundary probabilities. The solid black line represents the true probability landscape on the exact state space. The dashed red line represents the probability landscape computed from the truncated state space with buffer capacity $N$. The gray shaded area represents the true error due to state space truncation with buffer capacity $N$. The probability of copy number $N$ on the true landscape is the true boundary probability, and the probability of $N$ on the computed landscape is the boundary probability on the truncated state space. In this study, we show that the computed boundary probability on the truncated state space can be used to bound the true error from the above.  Figure 2: The birth-death system associated with the aggregated rate matrix $B$. Each circle represents an aggregated state consisting of all microstates with the same copy number of elementary molecules in the MEG. These aggregated states are connected by aggregated birth and death reactions, with apparent synthesis rates $\alpha_{i}$ and degradation rates $\beta_{i+1}$ (see Lemma 1).  Figure 3: Error quantification and comparisons for the birth-death model. (A): The a priori estimated error (red solid line), the computed error (green line and circles), and the true error (blue line and crosses) of the steady state probability landscape. The inset shows the ratio of the true errors to the computed errors at different sizes of the MEG, and the grey straight line marks the ratio one. The computed errors are larger than the true errors when the black line is below the grey straight line. (B): The steady state probability landscapes of $X$ obtained with different truncations of net molecular number in the MEG. Note that probability distributions end at $X$ where truncation occurs. The probabilities in the landscapes are inflated when truncating the state space at smaller net molecular numbers of the MEG.  Figure 4: Error quantification and comparisons for the single gene expression model. (A) and (B): The a priori estimated error (red solid lines), the computed error (green lines and circles), and the true error (blue lines and crosses) of the steady state probability landscapes of $mRNA$ and $Protein$ at different sizes of truncations. The insets in (A) and (B) show the ratio of the true errors to the computed errors at different sizes of the MEG, and the grey straight line marks the ratio one. The computed errors are larger than the true errors when the black line is below the grey straight line. (C): The steady state probability landscapes of $Protein$ solved using different truncations of net molecular number in the MEG2. Note that probability distributions end at where truncation occurs. The probabilities in the landscapes are significantly inflated when truncating the state space at smaller net molecular numbers of the corresponding MEG. (D): The steady state probability landscapes of $mRNA$ solved using different truncations of net molecular number in the MEG2. In this cases, the probabilities in the landscapes are not affected by the truncation of the opposite MEG.  Figure 5: Error quantification and comparisons for the genetic toggle switch model. (A) and (B): The a priori estimated error (red solid lines), the computed error (green lines and circles), and the true error (blue lines and crosses) of the steady state probability landscapes of $A$ and $B$ at different sizes of truncations. The insets in (A) and (B) show the ratio of the true errors to the computed errors at different sizes of the MEG, and the grey straight line marks the ratio one. The computed errors are larger than the true errors when the black line is below the grey straight line. (C) and (D): The steady state probability landscapes of $A$ and $B$ solved using different truncations of net molecular number in the MEG1 and MEG2, respectively. Note that probability distributions end at where truncation occurs. The probabilities in the landscapes are significantly inflated when truncating the state space at smaller net molecular numbers of the corresponding MEG. (E) and (F): The steady state probability landscapes of $A$ and $B$ solved using different truncations of net molecular number in the MEG2 and MEG1, respectively. The probabilities in the landscapes are also significantly inflated when truncating the state space at smaller net molecular numbers of the opposite MEG.  Figure 6: Error quantification and comparisons for the phage lambda bistable epigenetic switch model. (A) and (B): The a priori estimated error (red solid lines), the computed error (green lines and circles), and the true error (blue lines and crosses) of the steady state probability landscapes of $CI$ and $Cro$ dimers at different sizes of truncations. The insets in (A) and (B) show the ratio of the true errors to the computed errors at different sizes of the MEG, and the grey straight line marks the ratio one. The computed errors are larger than the true errors when the black line is below the grey straight line. (C) and (D): The steady state probability landscapes of $CI$ and $Cro$ dimers solved using different truncations of net molecular number in the MEG1 and MEG2, respectively. Note that probability distributions end at where truncation occurs. The probabilities in the landscapes are significantly inflated when truncating the state space at smaller net molecular numbers of the corresponding MEG. (E) and (F): The steady state probability landscapes of $CI$ and $Cro$ dimers solved using different truncations of net molecular number in the MEG2 and MEG1, respectively. The probabilities in the landscapes are also significantly inflated when truncating the state space at smaller net molecular numbers of the opposite MEG.  Table 1: Reaction scheme and rate constants in phage lambda epigenetic switch mode. We use $CORn$ denotes $Cro_{2}$ bound operator site $ORn$, $RORn$ denotes $CI_{2}$ bound $ORn$, where $n$ can be 1, 2, and 3. Note that molecular species enclosed in parenthesis are those whose presence is required for the specific reactions to occur, but their copy numbers do not influence the transition rates between microstates.  Reactions | Rate constants
---|---
Synthesis reactions [50, 59, 60, 61]
$\emptyset+(OR3+OR2)\rightarrow CI_{2}+(OR3+OR2)$ | $k_{sCI_{2}}=0.0069/s$
$\emptyset+(OR3+COR2)\rightarrow CI_{2}+(OR3+COR2)$ | $k_{sCI_{2}}=0.0069/s$
$\emptyset+(OR3+ROR2)\rightarrow CI_{2}+(OR3+ROR2)$ | $k_{s1CI_{2}}=0.069/s$
$\emptyset+(OR1+OR2)\rightarrow Cro_{2}+(OR1+OR2)$ | $k_{sCro_{2}}=0.0929/s$
Degradation reactions [62, 50]
$CI_{2}\rightarrow\emptyset$ | $k_{dCI_{2}}=0.0026/s$
$Cro_{2}\rightarrow\emptyset$ | $k_{dCro_{2}}=0.0025/s$
Association rate of binding reactions [63]
$CI_{2}+OR1\rightarrow ROR1$ | $k_{bOR1CI_{2}}=0.021/s$
$CI_{2}+OR2\rightarrow ROR2$ | $k_{bOR2CI_{2}}=0.021/s$
$CI_{2}+OR3\rightarrow ROR3$ | $k_{bOR3CI_{2}}=0.021/s$
$Cro_{2}+OR1\rightarrow COR1$ | $k_{bOR1Cro_{2}}=0.021/s$
$Cro_{2}+OR2\rightarrow COR2$ | $k_{bOR2Cro_{2}}=0.021/s$
$Cro_{2}+OR3\rightarrow COR3$ | $k_{bOR3Cro_{2}}=0.021/s$
Dissociation reactions - $CI_{2}$ dissociation from OR1
$ROR1+(OR2)\rightarrow CI_{2}+OR1+(OR2)$ | $0.00898/s$
$ROR1+(ROR2+OR3)\rightarrow CI_{2}+OR1+(ROR2+OR3)$ | $0.00011/s$
$ROR1+(ROR2+ROR3)\rightarrow CI_{2}+OR1+(ROR2+ROR3)$ | $0.01242/s$
$ROR1+(ROR2+COR3)\rightarrow CI_{2}+OR1+(ROR2+COR3)$ | $0.00011/s$
$ROR1+(COR2)\rightarrow CI_{2}+OR1+(COR2)$ | $0.00898/s$
Dissociation reactions - $CI_{2}$ dissociation from OR2
$ROR2+(OR1+OR3)\rightarrow CI_{2}+OR2+(OR1+OR3)$ | $0.2297/s$
$ROR2+(ROR1+OR3)\rightarrow CI_{2}+OR2+(ROR1+OR3)$ | $0.0029/s$
$ROR2+(OR1+ROR3)\rightarrow CI_{2}+OR2+(OR1+ROR3)$ | $0.0021/s$
$ROR2+(ROR1+ROR3)\rightarrow CI_{2}+OR2+(ROR1+ROR3)$ | $0.0029/s$
$ROR2+(COR1+OR3)\rightarrow CI_{2}+OR2+(COR1+OR3)$ | $0.2297/s$
$ROR2+(OR1+COR3)\rightarrow CI_{2}+OR2+(OR1+COR3)$ | $0.2297/s$
$ROR2+(COR1+COR3)\rightarrow CI_{2}+OR2+(COR1+COR3)$ | $0.2297/s$
$ROR2+(ROR1+COR3)\rightarrow CI_{2}+OR2+(ROR1+COR3)$ | $0.0029/s$
$ROR2+(COR1+ROR3)\rightarrow CI_{2}+OR2+(COR1+ROR3)$ | $0.0021/s$
Dissociation reactions - CI dissociation from OR3
$ROR3+(OR2)\rightarrow CI_{2}+OR3+(OR2)$ | $1.13/s$
$ROR3+(ROR2+OR1)\rightarrow CI_{2}+OR3+(ROR2+OR1)$ | $0.0106/s$
$ROR3+(ROR2+ROR1)\rightarrow CI_{2}+OR3+(ROR2+ROR1)$ | $0.0106/s$
$ROR3+(ROR2+COR1)\rightarrow CI_{2}+OR3+(ROR2+COR1)$ | $0.0106/s$
$ROR3+(COR2)\rightarrow CI_{2}+OR3+(COR2)$ | $1.13/s$
Dissociation reactions - Cro dissociation from OR1
$COR1+(OR2)\rightarrow Cro_{2}+OR1+(OR2)$ | $0.0202/s$
$COR1+(ROR2)\rightarrow Cro_{2}+OR1+(ROR2)$ | $0.0202/s$
$COR1+(COR2+OR3)\rightarrow Cro_{2}+OR1+(COR2+OR3)$ | $0.0040/s$
$COR1+(COR2+ROR3)\rightarrow Cro_{2}+OR1+(COR2+ROR3)$ | $0.0040/s$
$COR1+(COR2+COR3)\rightarrow Cro_{2}+OR1+(COR2+COR3)$ | $0.0040/s$
Dissociation reactions - Cro dissociation from OR2
$COR2+(OR1+OR3)\rightarrow Cro_{2}+OR2+(OR1+OR3)$ | $0.1413/s$
$COR2+(ROR1+OR3)\rightarrow Cro_{2}+OR2+(ROR1+OR3)$ | $0.1413/s$
$COR2+(OR1+ROR3)\rightarrow Cro_{2}+OR2+(OR1+ROR3)$ | $0.1413/s$
$COR2+(ROR1+ROR3)\rightarrow Cro_{2}+OR2+(ROR1+ROR3)$ | $0.1413/s$
$COR2+(COR1+OR3)\rightarrow Cro_{2}+OR2+(COR1+OR3)$ | $0.0279/s$
$COR2+(OR1+COR3)\rightarrow Cro_{2}+OR2+(OR1+COR3)$ | $0.053/s$
$COR2+(COR1+COR3)\rightarrow Cro_{2}+OR2+(COR1+COR3)$ | $0.0328/s$
$COR2+(ROR1+COR3)\rightarrow Cro_{2}+OR2+(ROR1+COR3)$ | $0.053/s$
$COR2+(COR1+ROR3)\rightarrow Cro_{2}+OR2+(COR1+ROR3)$ | $0.0279/s$
Dissociation reactions - Cro dissociation from OR3
$COR3+(OR2)\rightarrow Cro_{2}+OR3+(OR2)$ | $0.0022/s$
$COR3+(ROR2)\rightarrow Cro_{2}+OR3+(ROR2)$ | $0.0022/s$
$COR3+(COR2+OR1)\rightarrow Cro_{2}+OR3+(COR2+OR1)$ | $0.0008/s$
$COR3+(COR2+ROR1)\rightarrow Cro_{2}+OR3+(COR2+ROR1)$ | $0.0008/s$
$COR3+(COR2+COR1)\rightarrow Cro_{2}+OR3+(COR2+COR1)$ | $0.003/s$