\section{Introduction}

The low-energy physics of interacting fermions in one dimension (1D) is
determined by collective spin and charge density excitations that define what is
known as Luttinger liquid (LL) behavior \cite{Haldane1981, Giamarchi2003,
Fradkin2013}. Soon after high-temperature superconductivity was discovered in
cuprate oxides \cite{Bednorz1986}, it was proposed that the charges added upon
doping a Mott insulator could end up distributed in stripes \cite{Tranquada1994,
Zaanen1989, Machida1989, Kato1990, Emery1999}. This led Anderson \emph{et al.}{} to
suggest that confined fermionic excitations in such presumed LL arrays (the
stripes) could explain the non-Fermi-liquid nature of cuprates' ``normal'' state
\cite{Strong1994}.
Since then, theoretical investigation has assessed whether LL behavior can
emerge in higher dimensions, especially in 2D \cite{Strong1994, Wen1990,
Emery2000, Vishwanath2001, Mukhopadhyay2001, Mukhopadhyay2001a}, the natural
route to that having been to study systems of coupled LLs in different guises.
It is now known, for example, that, in an array of parallel LLs, marginal
inter-wire density-density and current-current interactions lead to strong
transverse charge-density fluctuations at incommensurate wave vectors which can
frustrate electron crystallization and indeed stabilize a LL state, commonly
designated ``smectic metal'' or ``sliding Luttinger liquid'' state
\cite{Emery2000, Vishwanath2001, Mukhopadhyay2001, Mukhopadhyay2001a}.
However, previous work has been limited to exploring consequences of couplings
among either one or two perpendicularly crossed arrays, without ever considering
LLs interlinked in the form of a triangular net, possibly for lack of a
realistic representative system.

The recent discovery of strongly correlated physics in marginally twisted
bilayer graphene (MTBG) near the magic angle $\theta {\,\approx\,} 1.1^{\circ}$
\cite{Cao2018, Cao2018a} set off a flurry of interest in the origin of the
observed insulating and superconducting (SC) phases. (For the purposes of this
work, MTBG refers to bilayers twisted by ${\,\sim\,} 1^{\circ}$ or less, including
the first magic angle.)
It had been previously suggested that, at magic angles, the quasi-flatness of
the electronic bands closest to the undoped Fermi level could promote electronic
instabilities \cite{SuarezMorell2010, Li2010, Bistritzer2011,
LopesdosSantos2012}.
The development of effective tight-binding models for those bands
\cite{LopesdosSantos2007, Koshino2018, Po2018, Kang2018} enabled, on the one
hand, predictions of possible broken symmetries arising from weak-coupling
mechanisms, such as Fermi surface nesting or enhanced density of states
\cite{Liu2018, Isobe2018, You2019, Po2018, Laksono2018}; on the other hand, it
revealed trilobed Wannier functions centered at the AB/BA positions
\cite{Koshino2018, Po2018, Kang2018}, which has in turn motivated
strong-coupling perspectives based on extended and non-conventional Hubbard-type
interactions \cite{Xu2018, Xiao2018, Gu2019}. Electronic interactions are an
undisputed factor given that the ratio of the local Coulomb integral to
bandwidth is estimated in the range $U/w {\,\sim\,} 5\text{--}10$
\cite{Koshino2018}.
The extremely large Moir\'e unit cells involved (${\,\sim\,} 172$\,nm$^2$, about 100
times those of canonical Mott-insulators like cuprates) has also prompted the
suggestion that the insulating phase can be a Wigner crystal, consistently
with the extremely low densities, and the emergence of SC a result of its
melting \cite{Padhi2018, Padhi2019}.

Evidence accumulated from recent experiments and theoretical work motivates
a different perspective over the effective electronic model governing
correlations in MTBG, which we develop in this paper.
We note, first, that, by allowing more than the minimum two orbitals
\cite{Po2018} per Moir\'e unit cell, Carr \emph{et al.}{} have recently shown that the
weight of the Bloch states belonging to the flat-band sector is overwhelmingly
distributed among Wannier functions situated at the AA positions \emph{and} at
the AB/BA \emph{domain boundaries} \cite{Carr2019b}.
Second, it is well known that, in the presence of perpendicular electrical
fields, AB-BA domain boundaries host protected helical modes \cite{Martin2008,
Kindermann2012, Zhang2013, Vaezi2013, Koshino2013}. MTBG accommodates well
defined, intrinsic, and periodically alternating AB/BA regions
\cite{LopesdosSantos2012} whose network of boundaries was shown to likewise
support the propagation of such confined states \cite{San-Jose2013, Efimkin2018,
Pal2018a}. Moreover, since AB is favored against AA stacking, a considerable
atomic relaxation within the Moir\'e unit cell maximizes the AB/BA regions,
leaving sharply defined, atomic-scale domain boundaries \cite{Alden2013,
Nam2017, Anelkovic2017, Gargiulo2018, Lucignano2019, Walet2019a}. Crucially,
there is now unequivocal spectroscopic \cite{Alden2013, Ju2015, Yin2016,
Huang2018, Sunku2018} and transport \cite{Ju2015, Rickhaus2018, Yoo2018, Xu2019}
evidence of the reality of this network of 1D modes in MTBG, including in
single-gated devices.

There are additional hints that warrant a description in terms of such
``network of linked quantum wires'' to describe the observed correlated
behavior: the confinement of electrons to 1D naturally boosts correlations; a
phase diagram similar to that of magic-angle samples arises at other twist
angles under pressure \cite{Yankowitz2019}, in line with the expectation that an
inter-layer coupling enhanced by pressure would amplify the lattice relaxation,
in turn defining sharp domain boundaries and the emergence of the wire network
for more generic twists \cite{Rickhaus2018, Yoo2018, Xu2019}.

All the above aspects and observations call for an investigation of the
implications of a coupled-wire description of the low-energy physics of MTBG. Wu
\emph{et al.}{} have recently advanced arguments to justify the insulating and SC phases
in such a scenario, but only considering coupling at the wire intersections and
spin isotropic interactions within each wire \cite{Wu2019}.
By generalizing to MTBG the approach developed to study sliding LL phases, we
scrutinize not only the competition among SC and charge-density wave (CDW)
states, but also the emergence of Fermi liquid (FL) and sliding LL phases that,
stabilized by inter-wire interactions, might explain the experimental
progression of MTBG from an insulator to SC to a metal as density deviates
from
commensurate fillings. The SC order parameter acquires either $s$ or $d{\,\pm\,} id$
symmetry, depending on the Josephson coupling at the wire crossings. Moreover,
it suggests that these correlated phases could also happen in ``marginally''
($\theta {\,\ll\,} 1^\circ$) twisted bilayer graphene.
Finally, we emphasize that, independently of its direct relevance to MTBG,
this work reports the first detailed study of coupled quantum wires and sliding
LL phases in a triangular net geometry.

\begin{figure}
\centering
\includegraphics[width=0.47\textwidth]{Fig-1}
\caption{(a) Schematic of the local net spanned by the three coupled arrays of
quantum wires, labeled $\{1,2,3\}$. One \emph{array} consists of a set of
parallel, identically colored wires. The superposition of the three arrays
reproduces the net of AB-BA domain boundaries in MTBG.
(b) Domain of stability. The array's effective Luttinger parameter,
$\kappa(k_\perp)$, is positive in the shaded domain for all $k {\,\in\,}
(-\pi/d,\pi/d)$. The LL phase exists \emph{only} in the red-to-yellow region
surrounding the upper boundary. Dashed lines indicate the cuts chosen to
generate the phase diagrams in \Fref{fig:phase-diagram}.
}
\label{fig:schematic}
\end{figure}

\section{Coupled-wire model}

\subsection{Interconnected Luttinger liquid array}

Mirroring the experimental network of AB-BA boundaries, we consider three
families of quantum wires, each family consisting of an array of parallel
wires depicted by the same color in \Fref{fig:schematic}(a). Although each link
in MTBG supports two helical channels per spin along the forward and backward
direction, we study the simplified case of a single mode per wire, which should
capture the essential physics at play.
Each independent wire is a LL whose excitations are best described in the boson
formalism by a separated spin ($s$) and charge ($c$) Hamiltonian
(see Appendix for the convention of bosonization)
\cite{Haldane1981, Giamarchi2003, Fradkin2013},
\begin{align}\label{eq:LL}
  H = \sum_{\alpha = c,s} \int\!dx \ \frac{v_{\alpha}}{2}
  \left[ \frac{1}{K_{\alpha}} (\partial_x \theta_{\alpha})^2
  + K_{\alpha} (\partial_x \phi_{\alpha})^2 \right]\!,
\end{align}
where $\theta_\alpha(x)$ and $\phi_\alpha(x)$ are the conventional phase-field
operators, $K_{c/s}$ is the (inverse) charge/spin Luttinger parameter
($K_c\,\text{\small$\gtrless$}\,1$ for repulsive/attractive interactions), and
$v_\alpha$ defines the velocity of each excitation. Spin back-scattering adds
the term $2 g_s/(2\pi \alpha)^2\int\!dx \cos\left( 2\sqrt{2\pi} \phi_s
\right)$, where $\alpha$ is the (length) cutoff of the theory, to Hamiltonian
\eqref{eq:LL}, leading to a sine-Gordon-like action and spin couplings that
flow according to the equations \citep{Giamarchi2003}
\begin{equation}\label{eq:sine-Gordon}
  \frac{dg_s}{dl} = 2(1-\frac{1}{K_s})g_s,
  \qquad
  \frac{dK_s}{dl} = \frac{g_s^2}{2\pi^2 v_s^2}.
\end{equation}
As a result, when there is a spin gap, $K_s {\,\rightarrow\,} \infty$.

To describe each periodic array of parallel wires separated by $d$ as in
\Fref{fig:schematic}, one must include the long-wavelength (charge)
density-density and current-current interactions among wires in the fixed-point
Hamiltonian, as first noted by Emery \emph{et al.}\ \cite{Emery2000}. It can then be
shown
that the charge part of the action reads, in Fourier space,
\begin{multline}\label{eq:action-array}
  S = \sum_{q} \frac{v(k_{\perp}) k^2}{2} \left[
    \frac{|\theta_{c,q}|^2}{\kappa(k_{\perp})} + \kappa(k_{\perp})
    |\phi_{c,q}|^2 \right]
    \\
    + i \omega_n k \phi_{c,q}^* \theta_{c,q},
\end{multline}
where $q {\,\equiv\,} (\omega_n,k,k_{\perp})$ and $k/k_{\perp}$ is the momentum
along/perpendicular to the wires \cite{Emery2000, Vishwanath2001,
Mukhopadhyay2001}. The spin part, for $\theta_s$ and $\phi_s$, is obtained by
replacing $[ v(k_{\perp}), \kappa(k_{\perp}) ] \to [ v_s, K_s ]$.
Direct comparison with \Eqref{eq:LL} shows that an array of LLs effectively
behaves as a LL, the net effect of the inter-wire coupling being a Luttinger
parameter $\kappa$ that is now a $2\pi/d$-periodic function of $k_{\perp}$.
This periodicity justifies a Fourier expansion,
\begin{equation}
  \kappa(k_{\perp}) = K_0[
    1 + K_1 \cos(k_{\perp}d) + K_2 \cos(2k_{\perp}d) + \dots],
  \label{eq:kappa}
\end{equation}
which we shall use below with $K_{0,1,2}$ as free parameters
\cite{Vishwanath2001, Mukhopadhyay2001, Mukhopadhyay2001a}.

Each of the three LL arrays depicted in \Fref{fig:schematic} is assigned a
(superscript) label $j {\,\in\,} \{1,2,3\}$. Within each array $j$, we consider
the single-electron hopping ($t_\perp$) between nearest-neighboring wires, as
well as inter-wire CDW ($\mathcal{V}_n$) and SC ($\mathcal{J}_n$) singlet
interactions between $n$-th neighboring wires. These are described,
respectively, by the \emph{intra}-array, \emph{inter}-wire couplings
\begin{subequations}\label{eq:intra-array}
\begin{align}
&\mathcal{H}^j_{h} = t_{\perp} \sum_{l,\sigma} \sum_{\nu = \pm1}
    \psi^{j \dagger}_{l,\nu,\sigma} \psi^j_{l+1,\nu,\sigma} + \text{H.\,c.},  \\
&\mathcal{H}^j_{c,n} = \mathcal{V}_n \sum_{l,\sigma,\sigma',\nu}
    \psi^{j \dagger}_{l,\nu,\sigma}\psi^j_{l,-\nu,\sigma}
    \psi^{j \dagger}_{l+n,-\nu,\sigma'} \psi^j_{l+n,\nu,\sigma'}, \\
&\mathcal{H}^j_{sc,n} = \mathcal{J}_n \sum_{l,\mu,\nu}
    \psi^{j \dagger}_{l,\mu,\uparrow}\psi^{j \dagger}_{l,-\mu,\downarrow}
    \psi^j_{l+n,\nu,\downarrow}\psi^j_{l+n,-\nu,\uparrow} + \text{H.\,c.},
\end{align}
\end{subequations}
where $\psi^{i}_{l,\pm1,\sigma}$ is the field operator for a right/left-moving
electron of spin $\sigma$ in the $l$-th wire of array $i$. The wires are
also coupled at each intersection [white dots in \Fref{fig:schematic}(a)],
requiring us to consider the additional \emph{inter}-array hopping, CDW and SC
interactions:
\begin{subequations}\label{eq:inter-array}
\begin{align}
  &H^{i,j}_{h} = t \sum_{l,m,\sigma}\sum_{\mu,\nu}
    \psi^{i \dagger}_{l,\mu,\sigma} \psi^{j}_{m,\nu,\sigma} + \text{H.\,c.}, \\
  &H^{i,j}_{c} = \mathcal{V}_0 \sum_{l,m} \sum_{\sigma,\sigma'} \sum_{\mu,\nu}
    \psi^{i \dagger}_{l,\mu,\sigma}\psi^{i}_{l,\mu,\sigma}
    \psi^{j \dagger}_{m,\nu,\sigma'} \psi^{j}_{m,\nu,\sigma'}, \\
  &H^{i,j}_{sc} = \mathcal{J}_0 \! \sum_{l,m,\mu,\nu} \!
    \psi^{i \dagger}_{l,\mu,\uparrow}\psi^{i \dagger}_{l,-\mu,\downarrow}
    \psi^{j}_{m,\nu,\downarrow} \psi^{j}_{m,-\nu,\uparrow} + \text{H.\,c.}.
\end{align}
\end{subequations}

\subsection{Renormalization group equations}

Once all the couplings  in \Eqrefs{eq:intra-array} and \eqref{eq:inter-array}
are written in terms of the bosonic fields, we proceed by developing a
perturbative renormalization group (RG) analysis. To the lowest order, the flow
equations for the hopping ($t$, $t_\perp$), CDW ($\mathcal{V}_n$), and SC
($\mathcal{J}_n$) coupling parameters read:
\begin{subequations} \label{eq:RG-flow}
\begin{align}
  &\frac{d\mathcal{V}_n}{dl} = (2-\delta_{n,0}-\frac{1}{K_s}-\Delta_{C,n})
    \mathcal{V}_n, \\
  &\frac{d\mathcal{J}_n}{dl} = (2-\delta_{n,0}-\frac{1}{K_s}-\Delta_{S,n})
    \mathcal{J}_n, \\
  &\frac{dt_{\perp}}{dl} = \left[ 2-\frac{1}{4} (K_s+\frac{1}{K_s})
    - \frac{1}{4}(\Delta_{C,1}+\Delta_{S,1}) \right] t_{\perp}, \\
  &\frac{dt}{dl} = \left[ 1-\frac{1}{4} (K_s+\frac{1}{K_s})
    - \frac{1}{4} ( \Delta_{C,0}+\Delta_{S,0} ) \right] t,
    \label{eq:RG-flow-t}
\end{align}
\end{subequations}
where,
\begin{subequations}\begin{align}
  \Delta_{C,n} & \equiv \int_{-\pi}^{\pi} \frac{dk}{2\pi} [ 1 - (1-\delta_{n,0})
  \cos(nk) ] \frac{1}{\kappa(k/d)} ,
  \\
  \Delta_{S,n} & \equiv \int_{-\pi}^{\pi} \frac{dk}{2\pi} [ 1 - (1-\delta_{n,0})
  \cos(nk) ] \, \kappa(k/d).
\end{align}\end{subequations}
It is physically reasonable to expect the intra-array couplings to decay
rapidly so, henceforth, we only consider intra-array CDW and SC interactions up
to second-neighbors. As for $\kappa(k_{\perp})$, in line with Vishwanath and
Carpentier \cite{Vishwanath2001}, we truncate its Fourier expansion at the
second order. Furthermore, in order to have a stable theory, $\kappa(k_{\perp})$
must be positive for $k_{\perp} \in (-\pi/d, \pi/d)$, which constrains $K_0 {\,>\,}
0$ and $(K_1, K_2)$ to the shaded domain shown in \Fref{fig:schematic}(b).

\begin{table}
\centering
\begin{tabular}{| c || c | c | c | c | c |}
\hline
Coupling & $( 0, K_{-} ) $ & $( K_{-}, 1/2 )$ & $( 1/2, 1 )$ &
$( 1, K_{+} )$ & $( K_{+}, \infty )$ \\ [0.5ex]
\hline
\multicolumn{6}{c}{intra-array (among parallel wires)} \\
\hline
$\mathcal{V}_n$ (CDW) & -- & -- & \ding{51} & \ding{51} & \ding{51} \\
$\mathcal{J}_n$ (SC) & -- & -- & \ding{51} & \ding{51} & \ding{51} \\
$t_\perp$ (hop) & -- & \ding{51} & \ding{51} & \ding{51} & -- \\[0.5ex]
\hline
\multicolumn{6}{c}{inter-array (at wire crossings)} \\
\hline
$\mathcal{V}_0$ (CDW) & -- & -- & -- & \ding{51} & \ding{51} \\
$\mathcal{J}_0$ (SC) & -- & -- & -- & \ding{51} & \ding{51} \\
$t$ (hop) & -- & -- & -- & -- & -- \\[0.5ex]
\hline
\end{tabular}
\caption{
Relevance of each coupling for the different ranges of the spin Luttinger
parameter ($K_s$) specified in the first row. The symbol \ding{51}{} means that a
coupling \emph{may} be relevant while ``--'' indicates it is \emph{always}
irrelevant within that interval of $K_s$. $K_{-} {\,\equiv\,} 3 {\,-\,} 2\sqrt{2}
{\,\simeq\,} 0.17$, $K_{+} {\,\equiv\,} 3 {\,+\,} 2\sqrt{2} {\,\simeq\,} 5.83$.
}
\label{table:phases}
\end{table}

At this level of approximation, the RG equations \eqref{eq:RG-flow} are
independent. The relevancy of the different couplings can thus be immediately
established and is summarized in \Tabref{table:phases}.
Since $\kappa$ is strictly positive, $\Delta_{C,0}+\Delta_{S,0} {\,\ge\,} 2$ which,
according to \Eqref{eq:RG-flow-t}, implies that the single-electron hopping at
the wire intersections ($t$) is, at most, marginal if $K_s {\,=\,} 1$ and
$\kappa(k/d) {\,=\,} 1$ for all $k$; it is otherwise irrelevant in nearly the whole
phase space. This justifies considering $t$ globally irrelevant and,
accordingly, it will be ignored in the subsequent analysis.
Similarly, one can see that $\Delta_{C,1}+\Delta_{S,1} {\,\ge\,} 2$ so that the
intra-array hopping ($t_{\perp}$) may be relevant when $3-2\sqrt{2} {\,\le\,} K_s
{\,\le\,} 3+2\sqrt{2}$. The CDW and SC couplings are relevant only if $K_s {\,>\,} 1/2$
in the intra-array case ($\mathcal{V}_{1,2}$ and $\mathcal{J}_{1,2}$), while the
corresponding inter-array couplings ($\mathcal{V}_0$ and $\mathcal{J}_0$) are
relevant for $K_s {\,>\,} 1$.

Up to this point, the spin Luttinger parameter $K_s$ has been considered free;
\Tabref{table:phases} thus covers the most general scenario in relation to the
possible magnetic phases. However, addition of the spin backscattering term
mentioned earlier to \Eqref{eq:LL} makes $K_s$ a running coupling, governed by
the flow \Eqrefs{eq:sine-Gordon}. The solution where $K_s {\,\rightarrow\,} \infty$
corresponds to a spin-gapped state, in which case we find the single-electron
hoppings $t$ and $t_{\perp}$ to be irrelevant (last column of
\Tabref{table:phases}), in correspondence with previous calculations for a
single array of coupled quantum wires \cite{Emery2000}. In contrast, if $K_s
{\,\rightarrow\,} 0$ we have a spin gapless state and all the couplings considered here
are irrelevant\,---\,the system consists of decoupled LLs.

\begin{figure}
\centering
\includegraphics[width=0.45\textwidth]{Fig-2}
\caption{
Phase diagram along the horizontal cuts marked in \Fref{fig:schematic} and for
$K_s {\,=\,} 2$: (a) $K_2 {\,=\,} 0.5$, (b) $K_2 {\,=\,} 0$. The parameters $K_0$ and
$K_1$ (axes) are the Fourier coefficients defined in \Eqref{eq:kappa}.
In the region above the solid-blue line, one of $\mathcal{J}_1$ or
$\mathcal{J}_2$ is relevant (SC order). In the region below the solid-orange
line, at least one of  $\mathcal{V}_1$ and $\mathcal{V}_2$ is relevant (CDW
order). The intra-array hopping ($t_\perp$) is relevant in the gray domain
bounded by the dash-dotted line, implying that the system might be a Fermi
liquid in this region. To the left of the blue-dashed line, the inter-array SC
coupling is relevant, whereas the inter-array CDW coupling is relevant to the
right of the orange-dashed line. The green area indicates a regime where all the
couplings are irrelevant, corresponding to a 2D LL state.
The main difference between (a) and (b) is the absence of the LL phase in the
latter.
}
\label{fig:phase-diagram}
\end{figure}

\section{Phase diagram and analysis}

\subsection{Instability tendencies}

While one may explore any range of $K_s$, we will now focus on $K_s=2$.
\Tabref{table:phases} shows that this falls in the regime where all couplings
but $t$ are relevant and, therefore, it is representative of the physical
scenarios involving phase competition, as is the case of MTBG, either driven by
inter- or intra-array interactions (or both). Figure \ref{fig:phase-diagram}
shows the phase diagram in two representative scenarios, defined by different
magnitudes of the second harmonic in the Fourier expansion \eqref{eq:kappa}.
Although $K_0$ is not strictly the Luttinger (charge) parameter of an individual
wire, \Eqref{eq:action-array} implies it does represent the effective Luttinger
parameter of an array behaving collectively as a LL \cite{Emery2000}. Therefore,
$K_0 {\,<\,} 1$ signals an effectively attractive regime while $K_0 {\,>\,} 1$
describes repulsion. In this context, one qualitatively understands the fact
that the SC phase (blue region) dominates in the small-$K_0$ portion of the
phase diagram, while the CDW eventually becomes the only relevant phase for
large $K_0$.
In the crossover region $K_0 {\,\sim\,} 1$, the domains of relevancy for the CDW and
SC orders overlap; in addition, the intra-array hopping is relevant as well in
this case (gray region enclosed by the dot-dashed line) which, should the
hopping become dominant over the CDW and SC instabilities, implies the existence
of a 2D FL phase. This indicates that the transition between SC and CDW with
increasing repulsion (increasing $K_0$) can occur either directly or via an
intervening FL phase, depending on the magnitude of $K_1$ (which is a measure of
the nearest-neighbor inter-wire coupling within an array).
The precise outcome of this phase competition, or coexistence, depends on how
the running of one coupling constant affects the others, whose analysis requires
a perturbative RG calculation beyond first order, which is  not in the scope of
this paper.

When $K_2$ is finite, in addition to a FL, one finds that a LL phase is
stabilized between the SC and CDW region close to the parameter-space boundary
at $K_1 {\,\rightarrow\,} \sqrt{2}$. This is marked by the green area in
\Fref{fig:phase-diagram}(a) or the orange region in \Fref{fig:schematic}(b).
Physically, the appearance of a LL phase in this case arises from the fact that
$K_2$ promotes interaction between next-nearest-neighboring wires within an
array, which is detrimental to the stability of the CDW. As first pointed out by
Vishwanath and Carpentier \cite{Vishwanath2001}, when $K_2$ is included, near
the $(K_1,K_2)$ space boundary, the minimum of the $\kappa(k_\perp)$ is located
at some \emph{incommensurate} $k_\perp$ and the value is close to zero, which
indicate strong fluctuations of a transverse incommensurate CDW order [the
density correlation $\langle \phi^*_{c,k,k\perp} \phi_{c,k,k_\perp} \rangle
\propto 1/\kappa(k_\perp) $]. The \emph{incommensurate} CDW fluctuations then
destroy the crystallization so that all  couplings are irrelevant within that
region.
In contrast, for $K_2 {\,=\,} 0$, the most divergent transverse CDW is commensurate
as $K_1 {\,\rightarrow\,} 1$; i.e., $1/\kappa(k_{\perp}) {\,=\,} 1/K_0[1 {\,+\,}
\cos(k_{\perp} d)]$ diverges at $k_{\perp} {\,=\,} \pi/d$, so the
next-nearest-neighbor intra-array CDW coupling will crystallize the system and
there is no LL phase. Indeed, comparing \Frefs{fig:phase-diagram}(a) and
\ref{fig:phase-diagram}(b), we see that a LL phase is stabilized at the expense
of the CDW phase in the repulsive region ($K_0 {\,>\,} 1$), without much impact on
the region of SC stability. At large $K_0$, both intra-array (rightward off the
solid orange line) and inter-array (rightward off the dashed orange line) CDW
couplings are relevant\,---\,the electrons crystallize and an insulator ensues.

If $K_0 {\,<\,} 0.5$, the inter-array Josephson coupling is relevant. At each wire
crossing, the phases of the three SC order parameters couple via
\begin{equation}\label{eq:sc-coupling}
  \propto \mathcal{J}_0\sum_{i,j \in \{1,2,3\}}
    \!\!\!\!\! \cos(\varphi^i_{\mathbf{r}} - \varphi^j_{\mathbf{r}}).
\end{equation}
Assuming that the intra-array SC coupling promotes uniform SC within each array,
\Eqref{eq:sc-coupling} indicates that the global SC phase depends on the sign of
$\mathcal{J}_0$: if $\mathcal{J}_0 {\,<\,} 0$, the Josephson coupling favors
$s$-wave SC with all $\varphi_i$ equal; but, if $\mathcal{J}_0 {\,>\,} 0$, that
coupling is frustrated and will result in a $2\pi/3$ difference between the
phase of the SC order parameter of one array ($j$) with respect to the next
($j+1$). This originates a $d {\,\pm\,} id$ SC symmetry. A similar conclusion has
been drawn by Wu \emph{et al.}{} who have further considered triplet pairing and discuss
the additional possibility of $p{\,\pm\,} ip$ symmetry \cite{Wu2019}.

\subsection{Commensurability}

In a conventional (i.e., single) LL problem, the proximity to commensurate
electron densities is described by considering the Umklapp process within
each wire \cite{Giamarchi2003} which, in the notation of \Eqref{eq:LL}, has the
form
\begin{equation}
  \frac{g_{U}}{(2\pi \alpha)^2} \int dx \cos[2\sqrt{2\pi}\phi_c +(4k_F-G)x],
\end{equation}
where $G$ is a vector of the reciprocal superlattice. The couplings $g_{U}$ and
$K_c$ flow according to \Eqref{eq:sine-Gordon}, with the replacements $g_s {\,\to\,}
g_{U}$, $K_s {\,\to\,} K_c$ \cite{Giamarchi2003}. In the present case, however, the
effective Luttinger parameter $\kappa$ is a function of the transverse momentum
[cf. \Eqref{eq:action-array}] due to the marginal interactions between wires
within each array; this complicates the flow equations in the charge sector. We
proceed by assuming, as a first approximation, that the flow equations for
$g_{U}$ and $K_0$ behave analogously to those in \Eqref{eq:sine-Gordon}, in
which case we naturally obtain distinct behavior at and away from half-filling:
Our phase diagram in \Fref{fig:phase-diagram} indicates that, away from
half-filling, the system is a SC provided $K_0$ is not too large; at (or near)
half-filling, a large enough $g_{U}$ is able to drive the system to an
insulating state even for very small $K_0$. Such SC-to-insulator transition is a
general feature of the competing instabilities in a LL with commensurate
density, because the Umklapp terms provide a ``condensation'' energy gain that
ultimately makes the charge-gapped CDW state energetically favorable
\cite{Giamarchi1997,Giamarchi2003,Emery1979}.

This competition between SC and CDW insulating states bears directly on the
current experimental observations with magic-angle MTBG, which show the
ground-state to be either a FL at generic densities, a SC near commensurate
fillings, or an insulator at commensurability \cite{Cao2018, Cao2018a}\,---\,the
coupled LL scenario is consistent with such observation.
For quantitative comparisons in this regard, it is worth noting that the
electronic filling/density reported for the 2D experimental system needs to be
converted to 1D electronic densities by taking into account that, in the
coupled-wire picture, each Moir\'e unit cell contains three non-equivalent wires.
For example, increasing the electron density by one electron per Moir\'e unit
cell away from charge neutrality corresponds to adding $1/3$ electrons per
segment of each non-equivalent wire within that unit cell.

\section{Discussion}

\subsection{Landscape of correlated states}

The propagation of interacting electrons along the quantum channels provided
by the well defined AB-BA domain boundaries of small-angle MTBG
[\Fref{fig:schematic}(a)] provides a natural low-energy picture for the
emergence of competing SC and insulating states. At generic densities and
moderate Luttinger parameter ($|K_0| {\,\sim\,} 1$), we obtain SC and possibly
FL as the dominant phases, with SC stabilized even for \emph{repulsive} Coulomb
interactions ($K_0 {\,>\,} 1$), which is noteworthy; at commensurate densities,
the system is a charge-insulator. This holds both when $K_s {\,<\,} 1$
and $K_s {\,>\,} 1$, particularly in the spin-gapped regime ($K_s {\,\gg\,} 1$)
where the only qualitative difference is the possible loss of the FL phase at
very high $K_s$\,---\,this is significant for the model applicability to
MTBG, where a magnetic field has been seen to destroy the insulating state
\cite{Cao2018, Yankowitz2019}.

\subsection{Anomalous metallic behavior}

Most interestingly, we see that the interaction among parallel wires contributes
to stabilize both ``sliding'' and ``crossed sliding'' LL phases, thus extending
previous findings for square arrays \cite{Emery2000, Vishwanath2001,
Mukhopadhyay2001, Mukhopadhyay2001a} to this triangular geometry as well. These
phases are extremely interesting because, on the one hand, they define a regime
of metallic 2D transport underpinned entirely by Luttinger-liquid behavior and
interactions, with the consequence that physical observables scale anomalously
with temperature, size, and fields \cite{Mukhopadhyay2001}. On the other hand,
and as a result, these regimes of 2D transport are entirely different from that
of an effective circuit of independent 1D wires.
Perhaps most significantly for current experiments is the fact that charge
transport in these phases would have anisotropic fingerprints and an anomalous
temperature dependence, thus being a \emph{natural} candidate for the ``strange
metal'' behavior reported in magic-angle MTBG right above the temperatures
where the insulating and SC states disappear \cite{Cao2019}. Moreover, in this
picture, a symmetry breaking among the three equivalent wire arrays would
naturally impart both local and global electronic ``nematicity'', a feature that
has recently been inferred from high-resolution scanning tunneling
microscopy (STM) experiments \cite{Jiang2019, Kerelsky2019, Choi2019a, Xie2019}.

\subsection{The nature of the 1D wire net}

As the presence of 1D modes traveling along the AB/BA domain boundaries is a
decisive precondition for modeling MTBG in terms of the proposed network of
coupled Luttinger liquids, in the remainder we elaborate on their evidence so
far as well as on means through which they can be ensured.

It is now well established that MTBG undergoes considerable internal deformation
within each Moir\'e unit cell so as to maximize extension of the energetically
more favorable Bernal stacking at the expense of the AA-stacked regions. This
energetic tendency is constrained by frustration at the interface between AB and
BA regions and results in bilayers with uniform Bernal stacking essentially
everywhere, except at sharp AB/BA domain walls and AA vertices whose geometry is
depicted in \Fref{fig:schematic}(a).
Numerical calculations reveal this relaxation effect very clearly
\cite{Alden2013, Nam2017, Anelkovic2017, Gargiulo2018, Lucignano2019}: it
starts becoming prominent for twist angles below ${\,\sim\,} 2\text{--}3^\circ$
\cite{Wijk2015, Zhang2018b, Gargiulo2018} and is completely established for
$\theta{\,\lesssim\,} 1^\circ$, at which point the width of the domain
boundaries saturates at ${\,\sim\,} 6\text{--}9$\,nm (becomes independent of twist
amount for smaller angles) \cite{Alden2013, Yin2016, Gargiulo2018}.
This threshold is confirmed experimentally \cite{Yoo2018} and MTBG samples have
been reported with domain walls extending up to the micron scale while retaining
their atomic-scale width \cite{Alden2013}.

The electronic modes localized at the AB/BA domain boundaries of deliberately
\emph{biased} bilayer graphene are expected to behave as perfect 1D quantum
wires, so long as intervalley scattering remains unimportant \cite{Martin2008};
explicit bandstructure calculations have recently shown this to be indeed
realized in relaxed MTBL \cite{Walet2019a,Carr2019b}.
More importantly, this has been confirmed by measurements that probed
electric transport along isolated AB/BA boundaries extending over several
microns, which revealed the expected conductance quantization at ${\,\approx\,}
4e^2/h$ and orders of magnitude enhancement of the mean-free path associated
with these modes, when compared with that elsewhere in the sample \cite{Ju2015,
Li2016d}. Their confinement to the domain walls has been confirmed directly by
local STM and STS measurements \cite{Yin2016, Huang2018} as well as indirectly:
(i) by the observation of Aharanov-Bohm oscillations in magnetotransport with
spatial periods that correlate with closed paths along adjacent domain walls
\cite{Rickhaus2018, Xu2019}; (ii) by the saturation of resistance near
$h/(4e^2)$ \cite{Xu2019, Yoo2018} and metallic temperature dependence despite
increases in interlayer bias (i.e., resistivity saturation and metallic
temperature dependence with increasing bulk gap) \cite{Xu2019}; (iii) and by
the local enhancement of infrared optical conductivity at the AB/BA interfaces,
which is associated with the presence of the 1D modes \cite{Sunku2018}.

All of the above indicates that the coupled-wire model should provide an
adequate description of biased MTBG, where the tunable bulk gap and topological
character of the 1D states ensure the robustness of these electronic modes, in
the regime where the Fermi level remains within the bulk gap. In the case of
magic-angle MTBG, it would be interesting to experimentally investigate the fate
of the correlated insulator and SC states under a finite interlayer bias.

Whether these domain-wall-bound modes survive and remain influential at zero
interlayer bias can depend on the conditions of the substrate. For example:
these confined modes have been seen directly by STM/STS against a gapped Bernal
background of graphene bilayers deposited on graphite, without any electrostatic
bias \cite{Yin2016}; and it is know that, similarly to a graphene monolayer
\cite{Kindermann2012, Hunt2013, Chen2014b, Moon2014a, Jung2015}, the Moir\'e and
relaxation induced by boron nitride substrates generates a spectral gap for
certain crystallographic orientations in bilayer graphene as well
\cite{Mucha-Kruczynski2013, Moon2014a, Dean2013, Kim2018a}. Therefore, a bulk
gap that guarantees and stabilizes the 1D modes can be engineered with
appropriate substrate conditions.

Finally, we note that there is a strong pseudomagnetic field due to the lattice
relaxation, with magnitudes that might exceed 10\,T \cite{Wijk2015, Nam2017}.
As a result of the triangular shape of the AB and BA domains, that field is
quasi-uniform within the Bernal regions, but with opposite polarity\,---\,the
polarity sharply switches precisely along the domain walls. The combined effect
of large pseudomagnetic fields and abrupt polarity changes along the domain
walls is likely to efficiently confine snake-type chiral states
\cite{Jones:2017, Muller1992}. These would be chiral 1D modes of a different
kind, which do not require a bulk gap \cite{Jones:2017}.

\subsection{Conclusions}

Different experimental probes and theoretical work are persuasive enough of the
conclusion that the triangular array of coupled quantum wires illustrated in
\Fref{fig:schematic}(a) is the natural starting point to describe transport and
correlated states in biased MTBG. The accuracy of this picture increases with
larger bulk gaps and angles $\lesssim 1^{\circ}$, which ensures the sharpest
domain boundaries as well as sufficiently long channels between the AA vertices
for a valid LL description of each quantum wire (the inter-vertex distance is
14\,nm for $\theta {\,=\,} 1^\circ$). The facts that such modes have been equally
seen in unbiased experiments and that pseudomagnetic fields can themselves beget
additional 1D states of a different nature, suggest the relevance of this
description to unbiased devices as well. Indeed, the phenomenology of the
correlated states, which so far has been scrutinized only in the unbiased case,
tallies with the phase diagram arising from the coupled LL model, namely when it
comes to: the types of correlated states involved and their competition, the
sequence of phase transitions with doping, the association of CDW insulating
states with commensurate fillings, the existence of non-Fermi liquid metallic
states, and nematicity.

Seeing as the detailed mechanisms underpinning both the insulating and SC states
in MTBG remain an open problem, it is of utmost interest to experimentally
scrutinize the evolution of the correlated phase diagram in MTBG (at the magic
as well as smaller angles) \emph{as a function of the bulk gap} through
interlayer bias. This would place the system in the regime where our model
most reliably applies, while it would also assess its relevance to
the strictly unbiased case.

\emph{Note}\,---\,Recently, a preprint emerged with a
similar formulation  \cite{Chou2019}, but with more restricted applicability
directly to MTBG: it considers a wire net with $C_4$ symmetry rather than $C_6$,
assumes persistent local SC order in puddle regions that encompass the wire
intersections, and couplings are considered only at those intersections, without
intra-array interactions. In specific cases (parameter ranges) where the two
models can be compared, the conclusions agree.

\begin{acknowledgments}
CC thanks X.~Y.~Gu and J.~N.~Leaw for helpful discussions. This work was
supported by the National Research Foundation of Singapore under its
Medium-Sized Centre Programme.
\end{acknowledgments}

\bibliographystyle{apsrev4-1}


\section{Introduction} \label{introduction}
Human genome contains about 20,000 genes, and those genes undergo diverse expression dynamics in different tissues throughout a life from its development to ageing. In a given cell of any tissue, tens of thousands of various gene products interact with one another, forming an incredibly complicated network, to support certain living processes. Such a network is often referred to as the gene regulatory network (GRN), and it has been a hot research topic for decades in the fields of genomics and system biology.

Many computational approaches have been developed to model GRN, such as information theory models, boolean networks, differential equations models, Bayesian networks, Petri nets and so on \cite{DeJong2002a, Delgado2018, Steggles2007}. Among them, Petri nets present one of the most promising tools, as Petri nets not only excel in modelling concurrent dynamic systems, but also have a wide application community as well as a strong theoretical support \cite{Reisig1985, Murata1989, Steggles2007, Bordon2012}.

Standard Petri nets are able to specify clearly the structure and behaviour of a particular process. However, it's not convenient to compose larger nets from smaller ones as Petri nets suffer from a lack of compositionality and abstraction \cite{Ermel1996}. While complex systems such as GRN often need to abstract from internal transitions and focus on the communication behaviour between processes so that a holistic picture of the whole system could be captured. To address the question, putting Petri nets into the framework of category theory provides an optimal strategy. In fact, this line of research has been explored since 1990 and reached its climax in the recent years with the introduction and optimization of open Petri nets \cite{Meseguer1990, Ermel1996, BALDAN2005,Rathke2014, Baldan2015,Baez2017a, Baez2018}.

Standard Petri nets and their variants (except open Petri nets) have been employed to model GRN for a while \cite{Chaouiya2006, chaouiya2007petri, Chaouiya2008, Steggles2007, Ruths2008, Bordon2012, Liu2014, Liu2017, Hamed2018}. On the other hand, open Petri nets have been successfully applied in modelling electrical circuits \cite{Baez2015a}, Markov processes \cite{Baez2015}, chemical reaction networks \cite{Baez2017a}, etc. However, open Petri nets haven't been used for GRN so far. Therefore, in this work we attempt to explore the possibility and potential of an open Petri nets implementation for GRNs, with the motivation of bringing a better composability to GRN modelling.

\section{Gene regulatory networks} \label{grn}
  In order to support normal cellular functions of a living organism, genes interact with each other to carry out many crucial molecular processed. These interactions usually involve a regulation of the gene expression. If a change in the expression of gene X could induce a change in the expression of gene Y, we say that gene X regulates gene Y. This regulation can be either up-regulation or down-regulation, and we call it activation or inhibition, respectively. Together, those genes and their regulations form a complex network which is referred to as gene regulatory network (GRN) \cite{Chai2014}.

  GRN could be inferred from gene expression data through computational approaches such as ordinary differential equation, Boolean network, Bayesian network, etc. Each of those methods has its pro and con. Ordinary differential equation (ODE), as its name suggests, uses continuous variables to represent the concentrations of mRNA or protein molecules (the products of gene expression), and it uses differential equations to model the dynamics of concentration changes. Therefore, ODE suits best for an accurate analysis of non-linear features of GRN. However, it's relatively complex and the parameters for the equations are quite often unavailable in the biological reality. On the other hand, Boolean network uses boolean values to represent gene expression level as well as the relationships between genes. In a Boolean network, a gene is either on or off and the interaction between any pairs of genes is either active or inactive. Boolean network has the advantage of simplicity, while it suffers an inability to capture many important details of GRN. Probabilistic Boolean network is an extension of Boolean network which allows more than one transition Boolean functions and each of them is randomly selected during state transitions to update the target gene \cite{Chai2014}.

  A more sophisticated probabilistic GRN model is Bayesian network. Bayesian network combines probability and graph theory to model the qualitative properties of GRN. It depicts the GRN as a directed acyclic graph (DAG) G=(X,A), X contains a set of nodes $\{x_1, ..., x_n\}$ which represent the gene variables, and A contains the directed edges which represents the probabilistic dependence interactions among genes. Bayesian network carries out the inference of network structure and parameters through belief propagation, which takes both the prior biological knowledge and the expression data into account. Bayesian network has been successfully applied to improve the quality of predicted GRN \cite{Chai2014}.

\section{Petri net modelling of gene regulatory networks} \label{pnGrn}
Petri nets provide a complementary network modelling to the above mentioned traditional methods. Petri nets were introduced by C.A. Petri in 1962 to model distributed and concurrent systems \cite{petri1962kommunikation}. Petri nets combine a well defined mathematical theory with a graphical representation; the former allows a precise analysis of the system behaviour, while the latter provides an intuitive visualization of the structure and state changes of the system. Due to these advantages, Petri nets have been successfully applied in modelling many kinds of dynamic systems such as computer networks, communication systems, logistic networks, etc. \cite{Liu2014}

By definition, a Petri net is a bipartite directed graph containing places and transitions connected by directed arcs; arcs can only connect place to transition and vice versa. Figure \ref{pn_1} shows a simple Petri net with 2 places (P1, P2) and one transition (T1).

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_1"}
	\caption{A simple Petri net}
	\label{pn_1}
\end{figure}
\noindent

A place can hold tokens (in Figure \ref{pn_1}, P1 holds one token as denoted by a red dot inside the circle), and an arc has its capacity (1 by default, otherwise the capacity is marked on the arc), and transitions have neither capacity, nor tokens. The marking (or state) of a Petri net is its token assignments of places. When the number of tokens in each input place is equal or larger than the arc weight (capacity), a transition is able to fire. When a transition fires, the tokens in input places are transferred to the output places, resulting in a new marking of the net. Figure \ref{pn_2} illustrates the firing of a transition, which reveals that the firing of a transition removes tokens from its input places and add tokens to its output places according to the capacity of the corresponding arcs. Note that tokens may disappear or appear going through a transition \cite{petri1966communication}.
\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_2"}
	\caption{Firing of a transition}
	\label{pn_2}
\end{figure}

\noindent
As the theory and methodology of Petri nets have been developed for more than 50 years, using Petri nets to model GRN will enrich our ability to analyse the dynamic properties of GRN. Several attempts have been made to model GRN using Petri nets \cite{Chaouiya2004, Chaouiya2006, chaouiya2007petri, Chaouiya2008, Liu2014, Liu2017, Liu2017a, Liu2017b, Liu2017c, Liu2018, Liu2018a, Liu2018b}. Here we adopt a simple model based on \cite{Chaouiya2004, DraganBosnacki2004} and \cite{Ruths2008}. In this model we draw an arrow from gene g2 to gene g1 to represent that g2 activates g1 (Figure \ref{pn_3}a), and a blunt-end line from g2 to g1 will represent the inhibition of g1 by g2 (Figure \ref{pn_4}a). For these two elementary GRNs we may construct the corresponding Petri nets, where genes are modelled as places, and the interaction between two genes is modelled by a transition together with directed arcs that connect the places and transitions.

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_3"}
	\caption{Activation of gene g1 by g2}
	\label{pn_3}
\end{figure}

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_4"}
	\caption{Inhibition of gene g1 by g2}
	\label{pn_4}
\end{figure}

In both activation and inhibition cases, g2 may or may not be consumed. We use double sided arrow to indicate that g2 activates g1 without itself being consumed (Figure \ref{pn_3}c and \ref{pn_4}c), and one sided arrow otherwise (Figure \ref{pn_3}b and \ref{pn_4}b). For activation, we draw an arrow pointing towards g1 from transition ${T_{g1, g2}}$, while for inhibition, we draw the arrow the other way round. Figure \ref{pn_3} and Figure \ref{pn_4} also describe what happens when ${T_{g1, g2}}$ fires. For example, in Figure \ref{pn_3}b, before firing there is a token in g2 but none in g1, while after firing the token in g2 is consumed and g1 is activated as indicated by having a token in it.

Once we have these basic building blocks, we are able to convert a normal GRN into its corresponding Petri net. For example, if we have a GRN as shown in Figure \ref{pn_5}a, then the Petri net for it will look like that in Figure \ref{pn_5}b.
\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_5"}
	\caption{A simple GRN and its Petri net representation}
	\label{pn_5}
\end{figure}

\noindent
Note that in this GRN gene A activates gene B and C, gene B activates gene D, and gene C inhibits gene D while activates gene E. Except in the case of C activating E, where C is consumed, all other transitions don't consume the source gene.

Now that we obtain the Petri net modelling of a GRN, we shall be able to simulate and analyse the dynamics of the GRN \cite{Ruths2008}.

\section{A Categorical Treatment of Petri nets}
Unfortunately, standard Petri net lacks proper mechanisms for composition and abstraction, which makes it inadequate to model large or multi-scale GRN. The problem can be solved by viewing Petri net from the category theory perspective, i.e., treating Petri nets as monoids (\cite{Meseguer1990, Ermel1996}). In order to do that, we need to introduce two definitions, following Meseguer and Montanari's work:

(1) A standard Petri net is a triple $<S, T, F>$, where S is a set of places, T is a set of transitions, and $F: (S \times T)+(T \times S) \rightarrow \mathbb{N}$ is a multiset which describes the connections between places and transitions. Here $\mathbb{N}$ represents natural numbers, $\times$ is the Cartesian product of sets, and + denotes the disjoint union of sets.

(2) A graph G is a set T of arcs, a set V of nodes and two functions $\partial_0$ and $\partial_1$ from T to V ($\partial_0, \partial_1: T \rightarrow V $), called source and target, respectively. A graph morphism from G to G' is a pair of functions $<f, g>, f: T \rightarrow T', g: V \rightarrow V'$ such that:
$g \circ \partial_0 = \partial_0' \circ f$ and $g \circ \partial_1 = \partial_1' \circ f$.

In support of their motto "Petri Nets Are Monoids", Meseguer and Montanani argued that Petri nets can be viewed as ordinary, directed graphs equipped with parallel and sequential composition of transitions \cite{Meseguer1990}. Accordingly, they provided a new definition of Petri net:

A Petri net is a graph where the arcs in set are called transitions and where the set of nodes is the free commutative monoid $S^\oplus$ generated by the set of places S. The elements of $S^\oplus$ has the form $n_1a_1 \oplus n_2a_2 \oplus ...n_ma_m$, where $a_i \in S, n_i \in \mathbb{N}, i \in 1..m$, and m is the number of elements in S. The source and target functions go from T to $S^\oplus$ ($\partial_0, \partial_1: T \rightarrow S^\oplus $). And a Petri net morphism $<f, g>$ is similar to a graph morphism, except that g is now a monoid homomorphism, which respects the monoid structure of $S^\oplus$.

This definition allows a different representation of the same standard Petri net. For example, if we have a Petri net $<S, T, F>$ as shown in Figure \ref{pn_6} where $S = \lbrace g_1, g_2, g_3 \rbrace$, $T = \lbrace t \rbrace $, $F(g_1,t)=1, F(g_2,t) = 1, F(t, g_1)=1, F(t, g_3)=1$, then we can represent it as $S^\oplus = \lbrace g_1, g_2, g_3 \rbrace ^\oplus$, $T = \lbrace t \rbrace $ and  $\partial_0(t) = g_1+g_2$, $\partial_1(t) = g_1+g_3$.

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_6"}
	\caption{A simple GRN modelled in Petri net}
	\label{pn_6}
\end{figure}

Under such a representation, sequential and parallel compositions of Petri nets based on transitions could be obtained in a straightforward way. For example, a sequential composition of two Petri nets is illustrated in Figure \ref{pn_7}, where the Petri net in Figure \ref{pn_7}c is a composition of the Petri net in Figure \ref{pn_7}b after that in Figure \ref{pn_7}a.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.4]{"pn_7"}
\caption{Sequential composition of Petri nets}
\label{pn_7}
\end{figure}

On the other hand, a parallel composition of two Petri nets is shown in Figure \ref{pn_8}, in which we calculate the products of places and transitions from Figure \ref{pn_8}a and Figure \ref{pn_8}b, and got the result Petri net in Figure \ref{pn_8}c.

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.4]{"pn_8"}
	\caption{Parallel composition of Petri nets}
	\label{pn_8}
\end{figure}

\noindent
With these two compositions at hands, we will be able to combine small Petri nets into large ones for the modelling of large GRNs.

Alternatively, we could construct the "open" Petri net following Baez and Master's recent work (\cite{Baez2018}).

The basic idea underlies an Open Petri net is that certain places in the Petri net are designated as inputs or outputs via a cospan of sets, which allows tokens to flow in or out of the Petri net and therefore make it open. From the categorical theory point of view, the input and output sets are objects and the open Petri net is the morphism between them, together they form a symmetric monoidal category (\cite{Baez2017a, Baez2018}). Further, in order to better describe vertical composition between sets and horizontal composition between open Petri nets, a symmetric monoidal double category $\mathbb{O}pen(Petri)$ was introduced. $\mathbb{O}pen(Petri)$ has the following components:\\
\hspace*{1em}(1) objects: sets X, Y, Z, ... \\
\hspace*{1em}(2) vertical 1-morphisms: functions from set X to Y $(f:X \rightarrow Y)$ \\
\hspace*{1em}(3) horizontal 1-cells: open Petri nets $(P:X \nrightarrow Y)$, which are cospans in Petri of the form: \\

\[
\begin{tikzpicture}[commutative diagrams/every diagram]
\node (P0) at (90:2.3cm) {$P$};
\node (P1) at (90+72:2cm) {$LX$} ;
\node (P4) at (90+4*72:2cm) {$LY$};

\path[commutative diagrams/.cd, every arrow, every label]
(P1) edge node {$i$} (P0)
(P4) edge node[swap] {$o$} (P0);
\end{tikzpicture}
\]
here L is a functor from Set to Petri.

\hspace*{1em}(4) 2-morphisms: morphisms between open Petri nets. A 2-morphism in $\mathbb{O}pen(Petri)$ has the following form:

\[ \begin{tikzcd}
LX_1 \arrow{r}{i_1} \arrow[swap]{d}{Lf} &P_1 \arrow[swap]{d}{\alpha}  &LY_1 \arrow[swap]{l}{o_1} \arrow{d}{Lg}\\%
LX_2 \arrow{r}{i_2}& P_2 &LY_2 \arrow[swap]{l}{o_2}
\end{tikzcd}
\]

We can compose 2-morphisms vertically by usual composition of functions as well as horizontally via the pushout of cospans. In addition, we can also obtain the tensor products of two 2-morphisms \cite{Baez2018}. For example, if we have two 2-morphisms as follows:

\[ \begin{tikzcd}
LX_1 \arrow{r}{i_1} \arrow[swap]{d}{Lf} &P_1 \arrow[swap]{d}{\alpha}  &LY_1 \arrow[swap]{l}{o_1} \arrow{d}{Lg}\\%
LX_2 \arrow{r}{i_2}& P_2 &LY_2 \arrow[swap]{l}{o_2}
\end{tikzcd}
\]

\vspace*{1em}

\[ \begin{tikzcd}
LX_2 \arrow{r}{i_2} \arrow[swap]{d}{Lf'} &P_1 \arrow[swap]{d}{\beta}  &LY_2 \arrow[swap]{l}{o_2} \arrow{d}{Lg'}\\%
LX_3 \arrow{r}{i_3}& P_3 &LY_3 \arrow[swap]{l}{o_3}
\end{tikzcd}
\]

\noindent
we can compose them vertically to get the 2-morphism:
\[ \begin{tikzcd}
LX_1 \arrow{r}{i_1} \arrow[swap]{d}{L(f' \circ f)} &P_1 \arrow[swap]{d}{\beta \circ \alpha}  &LY_1 \arrow[swap]{l}{o_1} \arrow{d}{L(g' \circ g)}\\%
LX_3 \arrow{r}{i_3}& P_3 &LY_3 \arrow[swap]{l}{o_3}
\end{tikzcd}
\]

\noindent
And, for the following two 2-morphisms:\\

\hspace*{3em}
\begin{tikzcd}
LX_1 \arrow{r}{i_1} \arrow[swap]{d}{Lf} &P_1 \arrow[swap]{d}{\alpha}  &LY_1 \arrow[swap]{l}{o_1} \arrow{d}{Lg}  &LY_1 \arrow{r}{i_1'} \arrow[swap]{d}{Lg} &P_1' \arrow[swap]{d}{\alpha'}  &LZ_1 \arrow[swap]{l}{o_1'} \arrow{d}{Lh} \\%
LX_2 \arrow{r}{i_2}& P_2 &LY_2 \arrow[swap]{l}{o_2} &LY_2 \arrow{r}{i_2'}& P_2' &LZ_2 \arrow[swap]{l}{o_2'}
\end{tikzcd}\\

\noindent
we can compose them horizontally to get the following 2-morphism:
\[ \begin{tikzcd}
LX_1 \arrow{r}{i_1' \circ i_1} \arrow[swap]{d}{Lf} &{P_1+_{LY_1}P_1'} \arrow[swap]{d}{\alpha +_{Lg} \alpha'}  &LZ_1 \arrow[swap]{l}{o_1' \circ o_1} \arrow{d}{Lh}\\%
LX_2 \arrow{r}{i_2' \circ i_2} &{P_2+_{LY_2}P_2'} &LZ_2 \arrow[swap]{l}{o_2' \circ o_2}
\end{tikzcd}
\]

\noindent
Moreover, for the two 2-morphism like these:\\

\hspace*{3em}
\begin{tikzcd}
	LX_1 \arrow{r}{i_1} \arrow[swap]{d}{Lf} &P_1 \arrow[swap]{d}{\alpha}  &LY_1 \arrow[swap]{l}{o_1} \arrow{d}{Lg}  &LX_1' \arrow{r}{i_1'} \arrow[swap]{d}{Lf'} &P_1' \arrow[swap]{d}{\alpha'}  &LY_1' \arrow[swap]{l}{o_1'} \arrow{d}{Lg'} \\%
	LX_2 \arrow{r}{i_2}& P_2 &LY_2 \arrow[swap]{l}{o_2} &LY_2' \arrow{r}{i_2'}& P_2' &LY_2' \arrow[swap]{l}{o_2'}
\end{tikzcd}\\

\noindent
we can calculate their tensor product to get the following 2-morphism:
\[ \begin{tikzcd}
L(X_1+X_1') \arrow{r}{i_1+i_1'} \arrow[swap]{d}{L(f+f')} &{P_1 +P_1'} \arrow[swap]{d}{\alpha + \alpha'}  &L(Y_1+Y_1') \arrow[swap]{l}{o_1+o_1'} \arrow{d}{L(g+g')}\\%
L(X_2+X_2') \arrow{r}{i_2+i_2'}&{P_2+P_2'} &L(Y_2+Y_2') \arrow[swap]{l}{o_2+o_2'}
\end{tikzcd}
\]

\section{Open Petri nets for gene regulatory networks} \label{opnGrn}
Finally, we will illustrate with an example on how to compose two GRNs (in the form of open Petri nets). The example presented here is based on that created in \cite{Baez2018}. Suppose we have the following two GRNs:\\
\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.6]{"pn_9"}

	\label{pn_9}
\end{figure}

\noindent
To compose them horizontal via pushout, we could first stick them side by side through the set Y and get the following intermediate network: \\
\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.6]{"pn_10"}

	\label{pn_10}
\end{figure}

\noindent Then, we identify those places that are images of the same point in Y, group them together, and randomly choose one place to represent the group, removing set Y meanwhile \cite{Baez2018}. For example, point 2 in Y has images g3 and g6, so g3 and g6 are grouped together. Similarly, g4 and g6 are grouped together, and therefore we have g3, g4 and g6 all in the same group. Now we could merge those 3 places together, and we use g6 to represent this group. Also we merge g5 and g7 using g7 as representative place for them. We no longer need set Y and will obtain the composed open Petri net:
\\
\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.6]{"pn_11"}

	\label{pn_10}
\end{figure}

\section{Conclusion and future work}
In this preliminary work we attempt to model the gene regulatory networks (GRNs) with open Petri nets (OPNs) \cite{Baez2018}. OPN incorporates the power of both category theory and Petri nets, therefore it holds the great potential to become an ideal framework for systematic representation of GRNs. At this stage we only provide an example of how to compose two hypothetical GRNs written in the forms of OPNs. In the future we need to formalize on how to convert real GRNs into composable OPNs and further demonstrate that these OPN models facilitate the construction of large GRNs from smaller parts.

\bibliographystyle{apalike}