# Plaquette Models, Cellular Automata, and Measurement-induced Criticality

Hanchen Liu<EMAIL_ADDRESS>Xiao Chen<EMAIL_ADDRESS>Department of
Physics, Boston College, Chestnut Hill, Massachusetts 02467, USA

###### Abstract

We present a class of two-dimensional randomized plaquette models, where the
multi-spin interaction term, referred to as the plaquette term, is replaced by
a single-site spin term with a probability of $1-p$. By varying $p$, we
observe a ground state phase transition, or equivalently, a phase transition
of the symmetry operator. We find that as we vary $p$, the symmetry operator
changes from being extensive to being localized in space. These models can be
equivalently understood as 1+1D randomized cellular automaton dynamics,
allowing the 2D transition to be interpreted as a 1+1D dynamical absorbing
phase transition. In this paper, our primary focus is on the plaquette term
with three or five-body interactions, where we explore the universality
classes of the transitions. Specifically, for the model with five-body
interaction, we demonstrate that it belongs to the same universality class as
the measurement-induced entanglement phase transition observed in 1+1D
Clifford dynamics, as well as the boundary entanglement transition of the 2D
cluster state induced by random bulk Pauli measurements. This work establishes
a connection between transitions in classical spin models, cellular automata,
and hybrid random circuits.

††preprint: APS/123-QED

###### Contents

  1. I Introduction
  2. II Preliminaries
    1. II.1  Ground state physics
    2. II.2 Spin-flip symmetries
    3. II.3 Characteristics of symmetry operators
  3. III Phase transitions in RPMs
    1. III.1 Random triangular plaquette model (RTPM)
    2. III.2 Random X plaquette model (RXPM)
  4. IV Boundary physics and cellular automaton dynamics
    1. IV.1 Spin-flip symmetry and cellular automaton dynamics
      1. IV.1.1 Deterministic cellular automaton (DCA) dynamcis
      2. IV.1.2 Randomized cellular automata (RCA)
    2. IV.2 Dynamical phase transitions and the boundary physics
      1. IV.2.1  Free boundaries on both sides
      2. IV.2.2  One boundary free and the other fixed
    3. IV.3 Plaquette models, stabilizer state measurement, and boundary entanglement structure
      1. IV.3.1 Plaquette models and measuring the stabilizer state
      2. IV.3.2 Bulk measurement induced boundary entanglement
      3. IV.3.3 RXPM, MBQC boundary, and Clifford MIPT
  5. V Conclusion and discussion

##  I Introduction

Classical Ising spin models with multi-spin interaction have been extensively
studied in physics [1, 2, 3, 4, 5, 6]. These models represent a fundamental
tool for understanding non-trivial physical phenomena such as phase
transitions, critical behaviors, and glassy dynamics. The symmetries
associated with spin-flip operations play a crucial role in elucidating the
physics of these spin models.

The shape of the spin-flip symmetry varies among different Ising-type
classical spin models. In the standard Ising model with two-body interactions,
there exists a global spin-flip symmetry. This symmetry arises because the
energy remains invariant when all Ising spins in the model are flipped
simultaneously [1]. On the other end, models like the Ising gauge theory [6],
where the classical spins are located on the edges of the lattice, have local
spin-flip symmetries associated with each site. In this model, flipping spins
connected to a specific site does not alter the Hamiltonian. This form of
symmetry is also known as local gauge symmetry.

In between the “global” and “local” lies the subsystem symmetry and fractal-
shaped symmetry. The subsystem symmetry operates on low-dimensional subsystems
within the entire system. This symmetry is evident in models such as the
4-body plaquette Ising models [7, 8, 9], where it is possible to flip a line
or a plane of spins without altering the form of the Hamiltonian. This notion
of subsystem symmetry can be further generalized into fractal-shaped symmetry
[10], where the spin-flips form a fractal-shaped pattern. One famous example
is the Sierpinski triangular shape spin-flip symmetry in the Newman-Moore
model [3].

Motivated by previous works, this work introduces a class of two dimensional
$q$-body Ising models. We replace these multi-spin interactions randomly with
single-site terms. Specifically, we focus on the cases of $q=3$ and $q=5$, and
observe that in these random models, varying the probability $1-p$ of single-
site terms induces a phase transition in the ground state. We investigate this
transition by analyzing the ground state degeneracy and the associated spin
flip symmetry operators, which transform between different ground states. We
find that this transition can be characterized by a structural change in the
symmetry operator. When $p>p^{c}$, the symmetry operator becomes non-local and
can spread over the entire system. In contrast, when $p<p^{c}$, the symmetry
operator is localized within a finite region.

Furthermore, we demonstrate that these 2D random spin models can be
effectively treated as 1+1D randomized cellular automaton dynamics. Through
this approach, we show that the localization transition of the symmetry
operator can now be mapped to the absorbing phase transition. When $p>p^{c}$,
a symmetry operator starting from one boundary at time $t=1$ can reach the
other boundary even after a long time evolution. Conversely, when $p<p^{c}$,
this boundary symmetry operator quickly vanishes into the bulk. We find that
when $q=3$, this cellular automaton dynamics corresponds to a special limit of
the Domany-Kinzel cellular automaton(DKCA) [11, 12], while the dynamics with
$q=5$ is a second-order automaton [13] with some random constraints.
Interestingly, this $q=5$ model has a notable connection with the measurement-
induced entanglement phase transition (MIPT).

The MIPT has been extensively studied in various hybrid random circuits [14,
15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. These studies reveal a generic
phase transition from a highly entangled volume-law phase to a disentangled
area-law phase as the measurement strength or rate is varied. The universality
class of these transitions depends on the specific model. Due to the
randomness in these circuits, understanding these criticalities poses a
challenging problem. Numerically, the 1+1D hybrid random Clifford circuit has
been well studied due to its efficient simulation on classical computers.

Besides the $1+1$D hybrid circuits, investigations on the boundary phase
transition on holographic tensor networks [26], the random stabilizer tensor
networks [27], as well as those boundary criticalities induced by bulk
measurements on 2D quantum states [23, 28, 29], offer a fresh perspective on
comprehending the entanglement phase transitions.

With the aforementioned context in mind, in this paper, this paper rigorously
establishes the equivalence between the phase transition of the symmetry
operator in the 2D classical random spin model with $q=5$ and the boundary
entanglement phase transition induced by bulk Pauli measurements on 2D cluster
states. This boundary transition can be effectively mapped to a 1+1D hybrid
Clifford circuit. Numerical investigations confirm that the transitions in all
three models belong to the same universality class. Notably, the dynamics of
the symmetry operator offers a novel perspective for understanding the phase
transition in the 1+1D hybrid Clifford circuit.

The subsequent sections of this paper are structured as follows. In Section
II, we present the fundamental concepts of random plaquette models, including
symmetries and quantities designed to elucidate the structure of the symmetry
group, such as symmetry entropy and mutual information. Moving on to Section
III, we examine the ground state phase transitions in random plaquette models
without boundary, characterizing these transitions in terms of the structure
of the symmetry group. In Section IV.1, we discuss the boundary phase
transition of the random plaquette models defined on the cylinder, showing
that it can be treated as the dynamical phase transition in randomized
cellular automata. This section further establishes the equivalence between
Pauli measurement on a stabilizer state and the plaquette models. Moreover, we
argue that the phase transition of the random Ising model with $q=5$ falls in
the same universality class as the boundary phase transition induced by bulk
Pauli measurement on the 2D cluster state and the 1+1D measurement-induced
phase transition in random Clifford circuits. Section V devotes to the
conclusion and discussion.

##  II Preliminaries

We examine a class of $q$-spin Ising models defined on two-dimensional
lattices:

$H[\sigma]=-\sum_{m_{q}}J_{m_{q}}\prod_{i\in m_{q}}\sigma_{i},$ (1)

where $\sigma_{i}\in\\{-1,1\\}$ represents an Ising spin operator at $i$-th
site and $\prod_{i\in m_{q}}\sigma_{i}$ describes a group of $q$ neighboring
spins that are linked together, referred to as the $m_{q}$-th plaquette term.
The Hamiltonian sums over all these plaquette terms. We set the coupling
constant uniformly with $J_{m_{q}}=1$. It is worth noting that we could also
consider a random coupling strength with $J_{m_{q}}=\pm 1$, and the ground
state physics explored in this paper would remain the same. However, for
simplicity, we focus on $J_{m_{q}}=1$ in this work. One well-known instance is
the two-dimensional Ising model, where each plaquette term is a two-body
interaction, and its ground state is twofold degenerate. This model also
displays a finite-temperature phase transition.

In this work, we focus on two spin models with plaquette terms representing
three-body and five-body interactions, respectively, defined on a square
lattice. At each site, we associate a three/five-body interaction, as shown in
Fig. 1a. We then randomly replace some plaquette terms with single spin terms,
resulting in two random plaquette models (RPMs).

The first model is the random triangular plaquette model (RTPM), where each
plaquette term is drawn from a binary distribution:

$\mu(\prod_{i\in m_{q}}\sigma_{i})=\left\\{\begin{matrix}p_{3}&q=\Delta\\\
1-p_{3}&q=\bullet\\\ \end{matrix}\right.$ (2)

where $p_{3}$ is the probability of the three-body term denoted as
‘$q=\Delta$,’ and $1-p_{3}$ is the probability of taking the single site
$\sigma$, denoted as ‘$q=\bullet$,’ as illustrated in Fig. 1a. This model
converges to the Newman-Moore model in the limit $p_{3}\to 1$.

Another model we investigate is the random X plaquette model (RXPM), the
random couplings of which are taken from

$\mu(\prod_{i\in m_{q}}\sigma_{i})=\left\\{\begin{matrix}p_{5}&q=+\\\
1-p_{5}&q=\bullet\\\ \end{matrix}\right.$ (3)

where $p_{5}$ is the probability of the five-body coupling, denoted as
‘$q=+$,’ and $1-p_{5}$ is again the probability of taking the single site
$\sigma$, denoted as ‘$q=\bullet$,’ as illustrated in Fig. 1b.

(a) RTPM : $q=\Delta$ and $q=\bullet$

(b) RXPM: $q=+$ and $q=\bullet$

Figure 1:  Interaction pattern of the RPM models.

###  II.1  Ground state physics

Figure 2:  The phase diagram and snapshots of the symmetry operator of RTPM
and RXPM. Blue-colored sites host the spin-flip while the green background
represents zeros. The snapshots are taken at $p_{3/5}=1,0.9,0.5$, and the
spin-flip symmetry operator $g$ is randomly draw from the full symmetry group
$G$.

The ground states of these models must satisfy all of the $\prod_{i\in
m_{q}}\sigma_{i}=1$ defined in the Hamiltonian. There is one trivial ground
state where all spins $\sigma=1$. To find other non-trivial ground states, we
reformulate this problem in terms of a parity check problem. Through the
mapping $\sigma=(-1)^{x}$ (with $x\in\\{0,1\\}$), each term is mapped to a
parity check

$\sum_{i}P_{m_{q}i}x^{\color[rgb]{0,0,1}i}=0\ \mbox{mod}\ 2$ (4)

where $P_{m_{q}i}=1$ if site $i$ is in the $m_{q}$-th plauqette, otherwise
$P_{m_{q}i}=0$. The index $m_{q}/i$ runs from $1$ to $Q/N$, with $Q$ being the
number of plaquettes and $N$ being the number of sites.

This problem can then be compactly rewritten as

$\displaystyle PT=0\ \mbox{mod}\ 2,$ (5)

where $P$ is a $Q\times N$ parity check matrix and in RTPM/RXPM models, the
number of ones per row can either be 1 or 3/5. The solution space is spanned
by $T$, which is a $N\times M$ binary matrix defined as

$T=(\vec{x}_{1},\vec{x}_{2},\ldots,\vec{x}_{M}).$ (6)

where each column in the $T$ matrix represents an independent generator, and
$M$ is the total number of generators. Each ground state
$\sum_{i=1}^{M}\vec{x}_{i}^{\alpha_{i}}$ is labelled by a binary vector
$\vec{\alpha}=(\alpha_{1},\cdots,\alpha_{M})$ and there are $2^{M}$ such
states.

The $T$ matrix can be directly solved by the standard Gaussian elimination in
the following manner. By row reduction, properly relabeling the column indices
and removing the rows that are all zeros, we have

$P=(I_{r}\mid C),$ (7)

where $I_{r}$ denotes the $r\times r$ identity matrix and $C$ is a
$r\times(N-r)$ binary matrix, with $r=\rank P$. It is then obvious that

$T=\left(\begin{matrix}C\\\ I_{M}\end{matrix}\right),$ (8)

as

$(I_{r}~{}C)\left(\begin{matrix}C\\\
I_{M}\end{matrix}\right)=C+C=0~{}\mathrm{mod}~{}2,$ (9)

where $M=N-r$, and here and in the remaining parts we use $\rank$ to refer to
the binary rank. This algorithm can be realized in polynomial time. Similar
problems arise in classical error-correcting codes [30] and the XORSAT problem
[31]. In classical error-correcting codes, the $T$ matrix describes the logic
space of classical information, while in XORSAT, the solution space $T$
exhibits interesting phase transitions.

###  II.2 Spin-flip symmetries

We define the spin-flip symmetry $g\in G$ as a tensor product of local spin-
flips $X_{i}$

$g=\prod_{i}X_{i}$ (10)

that does not change the energy given by Eq. 1:

$H[\sigma]=H[g\sigma],$ (11)

where the spin-flip $X_{i}$ flips the classical spin on site $i$ as follows:

$X_{i}\sigma_{i}=-\sigma_{i}.$ (12)

If $\vec{x}_{m}$ is a solution of Eq. 4, then $g\vec{x}_{m}$ is also a
solution, since the spin-flip symmetry $g$ does not change the energy given by
Eq. 1 by definition. Thus, we can associate each $\vec{x}_{m}$ with a
$g_{m}\in G$ by the spin-flips that connect the trivial ground state
$\vec{x}_{0}$, where $\sigma_{i}=1$ for all $i$, and a non-trivial solution
$\vec{x}_{m}$ of Eq. 4:

$\vec{x}_{m}=g_{m}\vec{x}_{0},\quad g_{m}\in G.$ (13)

The spin-flip symmetry operator can now be written in a more specific form as

$g_{m}=\prod_{i}X_{i}^{x_{m}^{i}}$ (14)

where $x_{m}^{i}$ is the $i$-th entry of vector $\vec{x}_{m}$.

The ground state properties are now encoded in the structure of the spin-flip
symmetry group. The total number of group elements is now

$|G|=2^{M}$ (15)

as there are a total of $2^{M}$ ground states, and the number of independent
generators is obtained via

$\text{\\# of gen.}=\log|G|=M.$ (16)

Here and in the remaining parts of this work, we take $\log$ to be 2-based.
Another important aspect of $G$ is that it is abelian, which indicates that
for any subgroup $G_{A}\subset G$, we are able to find a group quotient
$G/G_{A}$, the number of elements of which satisfies

$\displaystyle\log|G/G_{A}|$ $\displaystyle=\log(|G|/|G_{A}|)$ (17)
$\displaystyle=\log|G|-\log|G_{A}|.$

and for any subregion $\mathcal{D}$

$\log|G|-\log|G_{\mathcal{D}}|=\rank T_{\overline{\mathcal{D}}}$ (18)

where $T_{\overline{\mathcal{D}}}$ is the submatrix of $T$ taking only entries
associated with domain $\overline{\mathcal{D}}$. These characteristics of
group $G$ will be extensively used in the following sections.

###  II.3 Characteristics of symmetry operators

In this section, we introduce several quantities to further characterize the
ground state spin configurations, or equivalently the structure of the
symmetry operators of the RPMs.

Firstly, we define the configuration entropy as the log of the dimension of
the full solution space $\vec{x}$ of Eq. 4, which is equivalent to the rank of
the matrix $T$:

$S_{\text{cf}}\equiv\rank T=\log|G|.$ (19)

This quantity also refers to the number of independent generators of the
symmetry group $G$.

To capture the non-locality of the symmetry operator, we introduce the
symmetry entropy with respect to a subregion $A$:

$\displaystyle\mathrm{sym}S_{A}$ $\displaystyle\equiv\log|G/(G_{A}\cdot
G_{\overline{A}})|$ (20)
$\displaystyle=\log|G|-\log|G_{A}|-\log|G_{\overline{A}}|$
$\displaystyle=\rank T_{A}+\rank T_{\overline{A}}-\rank T,$

where $G_{A/\overline{A}}$ is a subgroup of $G$ containing only elements
located inside domain $A$ or its complement $\overline{A}$, and
$T_{A/\overline{A}}$ is submatrix of $T$ containing entries associated only
with sites in region $A/\overline{A}$. The symmetry entropy
$\mathrm{sym}S_{A}$ thus counts the number of symmetry operators shared by
both $A$ and $\overline{A}$.

Similar to $\mathrm{sym}S_{A}$, we can define another quantity to characterize
the spatial structure of the symmetry operator. This quantity is defined as
follows:

$\displaystyle\mathrm{sym}I_{AB|C}$
$\displaystyle\equiv\log\left|G^{C}/(G_{A}^{C}\cdot G_{B}^{C})\right|$ (21)
$\displaystyle=\log|G_{C}|-\log|G_{AC}|-\log|G_{BC}|+\log|G|$
$\displaystyle=\rank T_{A}+\rank T_{B}-\rank T_{AB}$

where we divide the system into three subsystems A, B and C with
$G^{C}=G/G_{C}$, $G_{A}^{C}=G_{AC}/G_{C}$, and $G_{B}^{C}=G_{BC}/G_{C}$. Here,
in the second line, we use Eq. 17, and in the third line, we use Eq. 18
together with the fact that $\rank T=\log|G|$. This quantity
$\mathrm{sym}I_{AB|C}$ computes the number of symmetry operators that have
non-trivial support on both subsystems $A$ and $B$.

##  III Phase transitions in RPMs

(a) RTPM operator size

(b) RXPM operator size

Figure 3:  Typical operator size of the spin-flip symmetry operators on a
$L\times L$ torus. The operator size $w(G)$ of group $G$ is defined by the
Hamming weight of the representation vector $\vec{x}_{m}$ of a randomly chosen
$g_{m}\in G$, and $\langle\ldots\rangle$ denotes the average of the
interaction patterns of RPM models.

We place RTPM and RXPM models on the torus by implementing periodic boundary
conditions. We observe that both models exhibit phase transitions as the
probabilities $p_{3}/p_{5}$ are varied, which are characterized by the
structure of the symmetry operator. When $p_{3}=p_{5}=1$, both models possess
fractal-shaped symmetry operators, with one example shown in Fig. 2. As we
depart from this fractal point, the symmetry operator becomes more extensive,
occupying $O(L^{2})$ number of sites, as present in Fig. 3. Further reduction
of $p_{3}/p_{5}$ results in another phase transition. Specifically, for the
RTPM model, when $p_{3}<p_{3}^{c}$, the symmetry operator disappears. In
contrast, for the RXPM model, symmetry operators persist but only occupy a
finite number of sites. We proceed to analyze these phase transitions in both
models quantitatively using the quantities introduced in the previous section.

(a) RTPM $S_{\text{cf}}$

(b) RTPM $S_{\text{cf}}$ data collapse

(c) RXPM $S_{\text{cf}}$

(d) RXPM $S_{\text{cf}}$ data collapse

Figure 4:  The data analysis for configuration entropy $S_{\text{cf}}$ of RTPM
defined on a $L\times L^{z}$ torus with $z=1.69$ and RXPM defined on a
$L\times L$ torus.

###  III.1 Random triangular plaquette model (RTPM)

We perform numerical computations to determine the configuration entropy of
RTPM as a function of $p_{3}$. This quantity represents the logarithm of the
ground state degeneracy or the number of symmetry operators. We find that when
$p_{3}$ is large, the configuration entropy $S_{\text{cf}}$ remains an $O(1)$
constant and decreases to 0 as $p_{3}$ decreases. The critical point
$p_{3}=0.81$ is marked by crossing point of $S_{\text{cf}}$ calculated in
different $L$ present in Fig. 4a.

As shown in Fig. 4b, around the critical point, we observe that
$S_{\text{cf}}$ for different system sizes collapse onto a single function
according to the following scaling form:

$S_{\text{cf}}\sim f(L^{1/\nu_{3}}(p-p_{c}))$ (22)

where $p_{3}^{c}\sim 0.81$ and $\nu_{3}\sim 1.21$. Here we take the torus to
be of size $L\times L^{z}$ with the geometry shown in Fig. 5a. The parameter
$z=1.697$ is used to capture the anisotropy of RTPM in these two directions.
The physical meaning of $z$ and the numerical method to calculate this
quantity will be further discussed in Sec. IV.1. The successful collapse of
the data implies that in the limit $L\to\infty$,

$S_{\text{cf}}\sim\left\\{\begin{aligned} \text{const.}&\quad
p_{3}>p_{3}^{c}\\\ 0&\quad p_{3}<p_{3}^{c}\end{aligned}\right..$ (23)

The ground state has finite degeneracy when $p_{3}>p_{3}^{c}$ and has no
degeneracy when $p_{3}<p_{3}^{c}$.

Moreover, we explore the spatial structure of the symmetry operators by
computing $\mathrm{sym}I_{AB|C}$. This quantity is defined in Eq. 21 and
counts the number of symmetry operators which have non-trivial support on both
$A$ and $B$. Specifically, we take the geometry as shown in Fig. 5a, where $A$
and $B$ are two antipodal line-shaped domains. When $p_{3}>p_{3}^{c}$, all
generators of the spin-flip symmetry are non-local, and there are $O(1)$ of
them, resulting in $\mathrm{sym}I_{AB|C}=S_{\text{cf}}\sim O(1)$. Conversely,
for $p_{3}<p_{3}^{c}$, there is no ground state degeneracy, and hence no
symmetry operator is shared by $A$ and $B$, leading to
$\mathrm{sym}I_{AB|C}=0$.

Through the above analysis, we observe that $S_{\text{cf}}$ and
$\mathrm{sym}I_{AB|C}$ are identical in both phases. This is confirmed by the
result in Fig. 5b, where we take the same data collapse for
$\mathrm{sym}I_{AB|C}$. The non-local structure of the symmetry operator is
also confirmed by the independence of $\mathrm{sym}S_{A}$ as a function of the
subsystem size, as demonstrated in Fig. 6d.

(a) Antipodal line-shaped domains $A,B$, with $L_{A}=L_{B}=1$ for RTPM defined
on a $L\times L^{z}$ torus with $z=1.697$, where boundaries with the same
arrow are glued together

(b) $\mathrm{sym}I_{AB|C}$ data collapse of RTPM, with $p_{3}^{c}=0.81$,
$\nu_{3}=1.21$.

(c) Antipodal line-shaped domains $A,B$, with $L_{A}=L_{B}=2$ for RXPM defined
on a $L\times L$ torus, where boundaries with the same arrow are glued
together

(d) $\mathrm{sym}I_{AB|C}$ data collapse of RXPM, with $p_{5}^{c}=0.743$,
$\nu_{5}=1.3$.

Figure 5:  $\mathrm{sym}I_{AB|C}$ of strip domain $A,B$ and $C=\overline{AB}$.

###  III.2 Random X plaquette model (RXPM)

For RXPM model, which is isotropic in both the horizontal and vertical
directions, we place this model on an $L\times L$ square torus. As we tune
$p_{5}$ from $1$ to $0$, a phase transition occurs, as indicated by the
configuration entropy density $S_{\text{cf}}/L^{2}$ shown in Fig. 4d, which
suggests that

$S_{\text{cf}}\sim\left\\{\begin{matrix}O(1)&\quad p_{5}>p_{5}^{c}\\\
O(L^{2})&\quad p_{5}<p_{5}^{c}\end{matrix}\right..$ (24)

When $0<p_{5}<p_{5}^{c}$, the extensive configuration entropy suggests
exponential ground state degeneracy in this phase. The phase transition is
also captured by the scaling of the symmetry entropy of a subregion $A$:

$\mathrm{sym}S_{A}\sim\left\\{\begin{matrix}O(1)&\quad p_{5}>p_{5}^{c}\\\
O(l_{A})&\quad p_{5}<p_{5}^{c}\end{matrix}\right.,$ (25)

where $l_{A}$ denotes the boundary length of subregion $A$, as present in Fig.
6e. When $p>p_{5}$, the number of symmetry operators $|G|\sim O(1)$, and the
symmetry operator $g\in G$ is non-local, contributing to a constant symmetry
entropy of any subregion $A$. As for $p<p_{5}$, the symmetry operators become
local and there are exponential number of them, resulting in an area-law
scaling of $\mathrm{sym}S_{A}$. At criticality $p_{5}\sim 0.743$, we observe
that

$\mathrm{sym}S_{A}/L\sim\tilde{c}\log\frac{L}{\pi}\sin\left(\frac{\pi
L_{A}}{L}\right)+b$ (26)

with $\tilde{c}\sim 3.6\times 10^{-3}$ and $b$ being some non-universal
constant,as shown in Fig. 6.

(a) Shape of subregion $A$ of RXPM on $L\times L$ torus

(b) Shape of subregion $A$ of RTPM on $L\times L$ torus

(c) RXPM critical symmetry entropy scaling $\mathrm{sym}S_{A}/L$, with
$\tilde{c}=6.9\times 10^{-3}$

(d) Trivial scaling of RTPM symmetry entropy

(e) Half-system symmetry entropy ($L_{A}=L/2$) of RXPM.

Figure 6:  Data analysis of symmetry entropy $\mathrm{sym}S_{A}$ of subregion
$A$.

To further investigate the structure of the symmetry operators, we compute the
symmetry mutual information $\mathrm{sym}I_{AB}$ of two disjoint antipodal
regions shown in Fig. 5a. As shown in Fig. 5d, $\mathrm{sym}I_{AB}$ collapses
to

$\mathrm{sym}I_{AB}\sim g(L^{1/\nu_{5}}(p_{5}-p_{5}^{c}))$ (27)

with $\nu_{5}=1.3$ and $p_{5}^{c}\sim 0.743$. This finding suggests that when
$p_{5}<p_{5}^{c}$, $\mathrm{sym}I_{AB}$ approaches zero in the limit
$L\to\infty$, indicating that the symmetry operators are local.

From the above analysis, it is evident that although both the RTPM and RXPM
models exhibit phase transitions, they differ significantly, especially in the
phase where $p_{3/5}<p^{c}_{3/5}$ and at criticality. In the following
sections, we will deepen our understanding of the transitions in both models
by introducing various boundary conditions.

##  IV Boundary physics and cellular automaton dynamics

In the previous section, we examine the localization phase transition of bulk
symmetry operators in RPMs placed on a torus. Here, we shift our focus to
models placed on geometries with boundaries. Specifically, we consider RPMs on
cylinders with two boundaries, subject to various boundary conditions. Our
objective is to study the behavior of symmetry operators originating from the
boundary and how they can be used to characterize phase transitions. To
achieve this, we introduce a dynamical approach that maps this 2D model to 1D
random cellular automaton dynamics. Through this mapping, we demonstrate that
the localization phase transition can be equivalently described as a dynamical
absorbing phase transition.

We consider two types of boundaries: fixed and free. In the fixed boundary,
all boundary spins $\sigma$ must be $+1$, while in the free boundary, no
constraints are imposed on the boundary spins. The free boundary condition
effectively removes $L_{f}$ constraints from the parity check matrix $P$ in
Eq. 4, resulting in $O(2^{L_{f}})$ spin-flip symmetries originating from the
boundary, where $L_{f}$ is the length of the free boundary. We consider two
combinations of boundary conditions, as illustrated in Fig. 7: one with only
one boundary free and the other fixed, and another with both boundaries free.

(a)  free boundaries (b)  Half-free boundary

Figure 7:  RPM defined on cylinder with two boundaries. Black vertical lines
with the arrows are glued together. Blue/Red line refers to free/fix boundary
respectively. Along the fixed boundary, spins are required to be $\sigma=+1$.

###  IV.1 Spin-flip symmetry and cellular automaton dynamics

Prior studies have demonstrated that spin-flip symmetry operators of plaquette
models can be derived through specific one dimensional classical cellular
automata[10, 32, 33]. In this part, we first revisit the relationship between
deterministic cellular automata (DCA) and plaquette models without randomness.
Subsequently, we map the RTPM and RXPM models to two distinct types of
randomized cellular automata: the probabilistic cellular automaton (PCA) and
cellular automaton with random impurity (CAwRI).

####  IV.1.1 Deterministic cellular automaton (DCA) dynamcis

Previously, we have shown in Eq. 14 that each spin-flip symmetry operator in
$G$

$G=\bigg{\\{}g_{k}=\prod_{i,\tau}X_{i,\tau}^{x_{k,i}^{\tau}}\bigg{\\}}$ (28)

can be characterized by a binary vector $\vec{x}_{k}$ with $k$ running from
$1$ to $|G|$. Here $(i,\tau)$ labels a site in the 2D lattice. This can be
alternatively understood as a spacetime configuration generated by a cellular
automaton at time $\tau$ on lattice site $i$ .

We first start with the Newman-Moore model with 3-body interaction. The
constraint imposed in the symmetry operator can be expressed as

$x_{k,i}^{\tau+1}+x_{k,i}^{\tau}+x_{k,i+1}^{\tau}=0,~{}\forall k.$ (29)

This can be treated as an automata updating rule $x^{\tau}\to x^{\tau+1}$

$x_{k,i}^{\tau+1}=x_{k,i}^{\tau}+x_{k,i+1}^{\tau},$ (30)

which is the rule-18 automata [34] that generates the Sierpinski triangle.

Following similar logic, for the X plaquette model, we have the automata
updating rule

$x_{k,i}^{\tau+1}=x_{k,i-1}^{\tau}+x_{k,i}^{\tau}+x_{k,i+1}^{\tau}+x_{k,i}^{\tau-1}.$
(31)

In this dynamics, the value of $x^{\tau+1}$ is determined by the
configurations of the two preceding layers $x^{\tau}$ and $x^{\tau-1}$. This
automaton is the reversible second-order automaton shown in [13].

####  IV.1.2 Randomized cellular automata (RCA)

We now discuss the randomized cellular automata dynamics corresponding to RTPM
and RXPM. Here we take the spatial/temporal direction shown in Fig 8.

(a) RTPM update rule. 3-body interaction:
$x_{k,i}^{\tau+1}=x_{k,i}^{\tau}+x_{k,i+1}^{\tau}$; single-site term:
$x_{k,i}^{\tau+1}=0$ (b) RXPM update rule. 5-body X interaction:
$x_{k,i}^{\tau+1}=x_{k,i}^{\tau}+x_{k,i+1}^{\tau}+x_{k,i-1}^{\tau}+x_{k,i+1}^{\tau}$;
single-site term: $x_{k,i}^{\tau+1}=0$

Figure 8:  Update rules for RPMs.

##### RTPM and Probabilistic cellular automaton (PCA) dynamics

For RTPM, if we treat the vertical direction as the time direction shown in
Fig. 8a, it can be mapped to a 1+1 dimensional classical cellular automata
dynamics. The update rule is defined as

$\mu(1|00)=0$ | $\mu(0|00)=1$
---|---
$\mu(1|10)=p_{3}$ | $\mu(0|10)=1-p_{3}$
$\mu(1|01)=p_{3}$ | $\mu(0|01)=1-p_{3}$
$\mu(1|11)=0$ | $\mu(0|11)=1$
(32)

where $\mu(x_{i}^{\tau+1}|x_{i}^{\tau}x_{i+1}^{\tau})$ denotes the probability
of site $i$ taking value $T_{i}^{\tau+1}$ at time $\tau+1$ under the condition
that $x_{i}^{\tau}$ and $x_{i+1}^{\tau}$ takes some certain value at time
$\tau$. This update rule is the same as the one in Wolfram rule 18 and is a
special limit of the Domany-Kinzel model [11, 12, 35]. We notice that
$\mu(0|00)=1$, effectively inhibiting the emergence of the symmetry operator
originating from the bulk.

##### RXPM and cellular automaton with random impurity (CAwRI)

For RXPM, we again take the vertical direction as the time direction shown in
Fig. 8b. When the random interaction takes the five-body form, the update rule
is

$x_{k,i}^{\tau+1}=x_{k,i}^{\tau}+x_{k,i+1}^{\tau}+x_{k,i-1}^{\tau}+x_{k,i}^{\tau-1}.$
(33)

The single-site term can be viewed as an impurity in the automaton dynamics,
imposing a constraint on the $T$ matrix, i.e.,

$x^{\tau}_{k,i}=0,~{}\forall(i,\tau)\in\Lambda_{I},$ (34)

where $\Lambda_{I}$ are these impurity sites. Furthermore, this constraint
leads to the generation of a new symmetry operator, starting at $(i,\tau+1)$.

To find the symmetry operators satisfying the above constraint, we propose an
algorithm to compute all generators of the symmetry group. Here’s the
approach: Assuming we have already obtained the symmetry group $G^{\tau}$ at
$t=\tau$ generated by $g_{k}^{\tau}$. We now want to add one more layer as
shown in Fig. 9 and compute the generators $g_{k}^{\tau+1}$ of $G^{\tau+1}$.
For the five-body interaction centered at $t=\tau$, $x_{k,i}^{\tau+1}$ follows
the update rules defined in Eq.(33). However, if we apply an on-site
$\sigma_{i,\tau}$, one requires $x^{\tau}_{k,i}=0$.

We denote generators that satisfy this constraint as
$\Omega^{1}=\\{g^{1}_{m}\\}$ and those that fail to meet this constraint as
$\Omega^{2}=\\{g^{2}_{n}\\}$. For the generators that do not meet this
constraint, we select the first one $g_{1}^{2}\in\Omega^{2}$ from the
unsatisfied set and modify the remaining generators to

$\tilde{g}^{2}_{n}=g_{1}^{2}g_{n}^{2}$ (35)

so that the associated $\tilde{x}^{\tau}_{n,i}=0$. We then remove $g_{1}^{2}$
from the $G^{\tau}$. Moreover, the on-site interaction induces a free spin at
$t=\tau+1$, introducing an extra spin flip symmetry $X_{i}$. The group
generators for the $\tau+1$ layer system are now

$G^{\tau+1}=\langle g^{1}_{m},~{}\tilde{g}^{2}_{n},~{}X_{i,\tau+1}\rangle.$
(36)

Note that the generators constructed in this manner may not be independent. To
resolve this, at the end of the evolution, it is necessary to apply Gaussian
elimination to eliminate any redundancies. In this algorithm, the initial
condition at $t=1$ is chosen as

$\displaystyle G^{1}=\langle X_{i,1},X_{i,2}\rangle$ (37)

with $i=1,\cdots L$. Here two layers of free spins are needed to evaluate the
third layer.

(a)  The dynamics up to the $(\tau-1)$th layer (b)  The dynamics at the
$\tau$th layer. The on-site term, denoted by the green dot on the $\tau$th
layer, induces a single site spin-flip $X$ on the $(\tau+1)$th layer denoted
as the yellow dot.

Figure 9:  The cellular automaton dynamics of the RXPM model.

###  IV.2 Dynamical phase transitions and the boundary physics

Due to the existence of the boundary, the symmetries $g\in G$ are now
categorized into two types: those with nontrivial support on the boundary and
those supporting trivially on the boundary. We define symmetries $g\in G$
supporting trivially on the boundary as the bulk symmetry $g^{bk}$. The bulk
symmetries form a group $G^{bk}\subset G$. In Sec. III, we have analyzed the
transitions in these bulk operators. In this section, we are interested in
these operators with non-trivial support on the boundary and we study the
transitions of these operators. They form a group and can be obtained as the
quotient

$G^{bd}\equiv G/G^{bk}.$ (38)

The size of $G^{bd}$ can be obtained by

$\mathrm{rank}T^{bd}=\mathrm{rank}T-\log|G^{bk}|=\log|G^{bd}|,$ (39)

following the same logic as in Eq. 18. One may regard $G^{bk}$ as the
disconnected contribution to the boundary. Here the submatrix $T^{bd}$ is
obtained by directly truncating the full tableau $T$ in Eq. 6, taking only
entries associated with the boundary.

For the cylinders with free boundary condition on both ends, we are interested
in these boundary operators which can have non-trivial support on both
boundaries. This form a subgroup of $G^{bd}$ and the number of symmetry
generators of this group is now given by

$\displaystyle\mathrm{sym}I$
$\displaystyle=\log\left|G^{bd}/\left(G_{\text{top}}^{bd}G_{\text{bot}}^{bd}\right)\right|$
(40)
$\displaystyle=\mathrm{rank}T_{\text{top}}^{bd}+\mathrm{rank}T_{\text{bot}}^{bd}-\mathrm{rank}T^{bd},$

where $T_{\mathrm{top/bot}}$ represents the submatrix of $T$ taking only
entries related to the top or bottom boundaries. From the viewpoint of
cellular automaton dynamics, the matrices $T_{\mathrm{top/bot}}$ are now the
initial/final state tableau $T_{1/\tau}$ of the cellular automaton dynamics.

The definition of boundary symmetry $G^{bd}$ also extends to the cylinder with
one fixed boundary shown in Fig. 7b, where the boundary symmetry operators
only start on the free boundary. These boundary symmetry operators can also be
used to characterize the phase transition. Analogous to the bulk symmetry
operator, we now define the boundary symmetry entropy as

$\displaystyle\mathrm{sym}S_{A}^{bd}\equiv\log\left|G^{bd}/\left(G_{A}^{bd}G_{\overline{A}}^{bd}\right)\right|$
(41) $\displaystyle=$
$\displaystyle\mathrm{rank}T_{A}^{bd}+\mathrm{rank}T_{\overline{A}}^{bd}-\mathrm{rank}T^{bd},$

which quantifies the number of group generators that have both nontrivial
support on the boundary domain $A$ and its boundary complement $\overline{A}$.
We can further define the boundary mutual information as

$\mathrm{sym}I_{AB}^{bd}=\mathrm{sym}S_{A}^{bd}+\mathrm{sym}S_{B}^{bd}-\mathrm{sym}S_{AB}^{bd}$
(42)

for two subsystems $A$ and $B$ on the boundary.

In what follows, we investigate the boundary phase transition using the
symmetry measures defined here. Through numerical analysis, we demonstrate a
transition from an active phase where boundary operators can penetrate deeply
into the bulk and exhibit extensive support there, to an absorbing phase where
boundary operators become localized near the boundary, as we vary the
parameter $p_{3/5}$ in the RPMs. Remarkably, the critical point $p_{3/5}^{bc}$
closely aligns with the one identified in the bulk localization transition,
suggesting that these transitions are indeed identical.

####  IV.2.1  Free boundaries on both sides

We begin by initializing a symmetry operator at the top boundary and analyze
its evolution over time.

As we increase the probability of the single-site term, the fraction of zeros
in this operator becomes more prevalent. Eventually, when $p<p^{c}$, the
entire operator rapidly decays to zero within the bulk of the system, and
becomes localized near the top boundary. We refer to this as the absorbing
phase.

In both models, we observe that for $p>p^{c}$, the operators can rapidly
spread throughout the entire system due to the update rule determined by the
three-body or five-body interaction. The occasional single-site term, which
enforces the corresponding site to remain zero, is unable to hinder the
operator’s growth. Even after a prolonged time evolution, the operator can
reach the bottom boundary. We refer to this as the active phase.

From this perspective, it is natural to introduce $\mathrm{sym}I$, as defined
in Eq. (40), between the top and bottom boundaries to characterize this
absorbing phase transition. Our numerical findings indicate that while both
models undergo phase transitions, they belong to distinct universality
classes.

(a)  $p_{3}=0.9$ (b)  $p_{3}=0.81$

(c)  $p_{3}=0.75$

Figure 10:  Space-time configuration of the RTPM-PCA evolution. Time flows
from top to bottom. Each blue site hosts a spin-flip operator. Gray sites
label the support of the full boundary group $G^{bd}$.

##### RTPM effective dynamical transition

The scaling of $\mathrm{sym}I$ with the effective time $L_{\tau}$ and boundary
length $L$ captures different phases of the RTPM. As shown in Fig. 12a,

$\mathrm{sym}I\sim\begin{cases}L&p_{3}>p_{3}^{c}\\\
L\exp(-dL_{\tau})&p_{3}<p_{3}^{c}\end{cases}$ (43)

where $d>0$ is some non-universal number, and $p_{3}^{c}=0.81$. When
$p_{3}>p_{3}^{c}$, the update rule corresponding to the three-body interaction
dominates. In this regime, symmetry operators originating from one boundary
can propagate to the other boundary even when $L_{\tau}\gg L$ (Fig. 10a). This
behavior leads to a plateau observed in Fig. 12a. On the other hand, for
$p_{3}<p_{3}^{c}$, the frequent single-site forces the symmetry operator to
take trivial values and therefore the operator starting from the boundary
quickly vanishes as shown in Fig. 10c. This results in an exponential decay in
$\mathrm{sym}I$, and thus the two boundaries are not connected by the boundary
symmetry operators when $L_{\tau}\gg L$.

The critical behavior of $\mathrm{sym}I$ is presented in Fig. 12b, where
$\mathrm{sym}I$ collapses to a universal scaling function:

$\mathrm{sym}I\sim h(L_{\tau}/L^{z}),$ (44)

where the scaling function takes the form

$h(\theta_{\tau})\sim\left\\{\begin{matrix}\pi
h^{0}\theta_{\tau}^{-1/z}&\theta_{\tau}\ll 1\\\
\exp(-h^{1}2\pi\theta_{\tau})&\theta_{\tau}\gg 1\end{matrix}\right.$ (45)

with dynamical exponent $z=1.697$, early time exponent $h^{0}=0.692$, and late
time exponent $h^{1}=0.075$. The dynamical exponent $z$ is very close to the
one computed in [36].

Figure 11:  RTPM with free boundaries (dark blue dots) on top and bottom
layers. The bulk three-body couplings are illustrated as green triangles and
the single-site terms are illustrated as the green dots.

(a) $\mathrm{sym}I$ of RTPM.

(b) critical scaling of $\mathrm{sym}I$ of RTPM.

Figure 12:  $\mathrm{sym}I$ between top and bottom boundaries.

(a)  $p_{5}=0.8$ (b)  $p_{5}=0.743$

(c)  $p_{5}=0.7$

Figure 13:  Space-time configuration of the RXPM-CAwRI evolution. Time flows
from top to bottom. Blue sites host spin-flips symmtries associated with the
boundary, the orange sites carry spin-flips related to the bulk, and the grey
sites support the full symmetry group.

##### RXPM effective dynamical transition

For RXPM, when the five-body coupling terms dominate, the symmetries are non-
local, and the top and bottom boundaries are connected by non-trivial symmetry
generators as shown in Fig. 13a. In contrast, in the area-law phase, the
symmetry operators are local, resulting in an exponential decay in
$\mathrm{sym}I$. The scaling of mutual information $\mathrm{sym}I$ between the
top and bottom boundary is again

$\mathrm{sym}I\sim\left\\{\begin{matrix}L&p_{5}>p_{5}^{c}\\\
L\exp(-dL_{\tau})&p_{5}<p_{5}^{c}\end{matrix}\right..$ (46)

with $d>0$ being some non-universal exponent. As shown in Fig. 15, at
criticality, $\mathrm{sym}I$ collapse to a universal scaling function

$\mathrm{sym}I\sim h(L_{\tau}/L),$ (47)

and the critical scaling function takes the form

$h(\theta_{\tau})\sim\left\\{\begin{matrix}\pi
h^{0}\theta_{\tau}^{-1}&\theta_{\tau}\ll 1\\\
\exp(-h^{1}2\pi\theta_{\tau})&\theta_{\tau}\gg 1\end{matrix}\right.$ (48)

Here, we have early time exponent $h^{0}=1.53$, and late time exponent
$h^{1}=0.125$.

Figure 14: RXPM with free boundaries (dark blue dots) on top and bottom
layers. The bulk five-body interactions are illustrated as green squares and
the single-site terms are illustrated as green dots.

(a) $\mathrm{sym}I$ of RXPM

(b) Early time critical scaling of $\mathrm{sym}I$ of RXPM. Red line:
$y=1.53(L_{\tau}/L)^{-1}$

(c) Late time critical scaling of $\mathrm{sym}I$ of RXPM. Red line:
$y\propto\exp(-0.125\times 2\pi L_{\tau}/L)$

Figure 15:  Symmetry mutual information between the top boundary and the
bottom boundaries.

####  IV.2.2  One boundary free and the other fixed

The boundary phase transition is also captured by the boundary symmetry
operator of a symmetry with half-free boundary condition, as shown in Fig. 7b.
In this case, the boundary symmetry only relates to the free boundary.

##### RTPM boundary transition

We investigate the boundary phase transitions in RTPM with half-free boundary
as shown in Fig. 16. The symmetry operators are evaluated via the randomized
updating rule presented in Eq. 32. The fixed boundary condition is realized by
imposing constraints on the state-vectors in the following manner.

Figure 16: RTPM with free boundary on top (dark blue dots) and fixed boundary
on bottom (red dots). The bulk three-body interactions are illustrated as
green triangles and the single-site terms are illustrated as green dots.

We encode the initial conditions of the RTPM automaton into a $L\times L$
matrix $T^{0}$, where each row vector in matrix $T^{0}$ represents an initial
condition. For simplicity, one can take $T^{0}=I_{L}$, the identity matrix,
where each initial condition only contains one non-trivial entry. Evolving
$T^{0}$ $L_{\tau}$ times, we obtain the final state matrix $T^{\tau}$, where
each row vector is related to the initial state vector. The full boundary
tableau

$T^{bd}=\left(T^{0}\mid T^{\tau}\right)$ (49)

comprises two parts, the initial state tableau $T^{0}$ and the final state
tableau $T^{\tau}$.

The fixed boundary condition requires that the generator must satisfy
$x^{\tau}=0$ at the final time. When $p_{3}<p_{3}^{bd}$, if $L_{\tau}$ is
large, the generator rapidly decays to zero in the bulk, automatically
satisfying $T^{\tau}=0$. However, when $p_{3}>p_{3}^{bd}$, the symmetry
operator starting from one boundary can reach the other boundary even after
the long time evolution. To identify the subspace of symmetry operators for
which $x^{\tau}=0$, we use the standard row reduction techniques. After the
row reduction, the state matrix takes the form

$T^{bd}=\left(\begin{matrix}Q^{0}&Q^{\tau}\\\ Q^{1}&0\end{matrix}\right),$
(50)

(a)  $p_{3}=0.9$ (b)  $p_{3}=0.81$

(c)  $p_{3}=0.75$

Figure 17:  Space-time configuration of the RTPM-PCA evolution. Time flows
from top to bottom. The bottom boundary is fixed, labeled by red. Each blue
site hosts a spin-flip operator, and the initial condition is a single spin-
flip $X$ reside in the middle. Grey sites label the support of the full
boundary group $G^{bd}$.

This phase transition is reflected in the scaling of $\mathrm{sym}I_{AB}$ and
entropy $\mathrm{sym}S_{A}$ of the boundary symmetries associated with the
free boundary. We find that the boundary symmetry $\mathrm{sym}I_{AB}$ of two
antipodal boundary domains $A,B$ with $L_{A}=L_{B}=L/8$ collapses to a
universal function

$\mathrm{sym}I_{AB}^{bd}\sim g((p_{3}-p_{3}^{c})L^{1/\nu_{3}^{\prime}})$ (51)

with $p_{3}^{c}\sim 0.81$ and $\nu_{3}^{\prime}\sim 2.43$, as shown in Fig.
18a. This collapse indicates a boundary phase transition as we tune the
coupling parameter $p_{3}$. Here the exponent $\nu_{3}^{\prime}$ is related to
$\nu_{3}\sim 1.21$ obtained previously in the following manner

$\nu_{3}^{\prime}=z\times\nu_{3}$ (52)

with $z=1.697$ being the effective dynamical exponent extracted in the
previous section.

When $p_{3}>p_{3}^{c}$, we have $\mathrm{sym}I_{AB}^{bd}>0$, indicating that
the boundary is long-range correlated and the boundary symmetry operators are
non-local. For $p_{3}<p_{3}^{bd}$, $\mathrm{sym}I_{AB}^{bd}=0$, meaning that
the boundary symmetry operators are local and the boundary is short-range
correlated. The boundary transition is further reflected in the scaling of
half-boundary symmetry entropy

$\mathrm{sym}S_{L/2}^{bd}\sim\left\\{\begin{matrix}L^{a}&p_{3}>p_{3}^{c}\\\
s_{0}&p_{3}=p_{3}^{c}\\\ \exp(-bL)&p_{3}<p_{3}^{c}\end{matrix}\right.$ (53)

with $s_{0}\sim 3$ and $a,b>0$, as shown in Fig. 18c. For $p_{3}>p_{3}^{c}$,
the non-local symmetry operators contribute to the power-law scaling of the
symmetry entropy, while at $p_{3}<p_{3}^{c}$, the symmetry entropy is
exponentially suppressed by system size $L$. In the active phase where
$p_{3}>p_{3}^{bd}$, $\mathrm{sym}S_{L/2}^{bd}$ exhibits a sublinear power-law
scaling with $\alpha\sim 0.5$. This finding is consistent with the scaling
behavior of the bulk symmetry operators on the cylinder geometry, indicating
that while they are non-local, they are not extensive. This is different from
the scaling behavior on the torus geometry as shown in Fig. 3b and the physics
of such a difference is left for the future work.

(a) RTPM boundary symmetry mutual information data collapse

(b) Geometry of the boundary symmetry mutual information $A$ and $B$, with
$L_{A}=L_{B}=L/8$

(c) RTPM boundary symmetry entropy scaling

(d) Geometry of half entropy $\mathrm{sym}S_{L/2}$. Here $L_{A}=L/2$

Figure 18:  Boundary phase transition of RTPM.

##### RPXM boundary transition

For RXPM, the free boundary is selected as shown in Fig. 19. Different from
that of RTPM, it contains two layers of spins. Again, we fix the last layer of
boundary spins to be $\sigma=+1$.

Figure 19: RXPM with free boundary on top (dark blue dots) and fixed boundary
on bottom. The bulk five-body interactions are illustrated as green squares
and the single-site terms are illustrated as green dots.

(a)  $p_{5}=0.8$ (b)  $p_{5}=0.743$

(c)  $p_{5}=0.7$

Figure 20:  Space-time configuration of the RXPM-CAwRI evolution. Time flows
from top to bottom. Blue sites host spin-flip symmetries associated with the
boundary, the orange sites carry spin-flips related to the bulk, and the grey
sites support the full symmetry group.

We observed a boundary phase transition as we tune the bulk coupling parameter
$p_{5}$, as present in Fig. 21b. When the five-body bulk couplings dominate,
$p_{5}>p_{5}^{c}$, the symmetries in the system are extensive and non-local,
resulting in a volume law scaling in the boundary symmetry entropy
$\mathrm{sym}S_{L/2}^{bd}\sim L$. On the other hand, when $p_{5}<p_{5}^{c}$,
the symmetry operator breaks into small clusters with finite length, resulting
in an area-law scaling of the boundary symmetry entropy
$\mathrm{sym}S_{L/2}^{bd}\sim O(1)$.

(a) Geometry of half entropy $\mathrm{sym}S_{L/2}^{bd}$

(b) $\mathrm{sym}S_{L/2}^{bd}$ of the RXPM

Figure 21: Boundary symmetry entropy on the top boundary of RXPM.

The critical point is further characterized by the peak of
$\mathrm{sym}S_{L/2}^{bd}$ as we tune the coupling parameter $p_{5}$. As shown
in Fig. 22b, $\mathrm{sym}S_{L/2}^{bd}$ collapses to a universal function

$\mathrm{sym}I_{AB}^{bd}\sim h((p_{5}-p_{5}^{c})L^{1/\nu_{5}})$ (54)

with $p_{5}^{c}=0.743$ and $\nu_{5}=1.3$.

(a) Geometry of anti-podal boundary domain $A,B$ with $L_{A}=L_{B}=L/8$

(b) Data collapse of $\mathrm{sym}I_{AB}^{bd}$ of RXPM

Figure 22: RXPM $\mathrm{sym}I_{AB}^{bd}$ data collapse.

At criticality $p_{5}=p_{5}^{c}=0.743$, the symmetry entropy scales
logarithmically

$\mathrm{sym}S_{A}^{bd}\sim c\log\sin\left(\frac{\pi L_{A}}{L}\right),$ (55)

with $c\sim 3.03$ as shown in Fig. 23a. The symmetry mutual information
$\mathrm{sym}I_{AB}^{bd}$ scales as

$\mathrm{sym}I_{AB}^{bd}\sim\chi_{AB}^{\Delta}$ (56)

with $\Delta=2$, and $\chi_{AB}$ being the cross-ratio

$\chi_{AB}=\frac{x_{12}x_{34}}{x_{13}x_{24}},\quad\mbox{with}~{}x_{ij}=\frac{L}{\pi}\sin\left(\frac{\pi}{L}|x_{i}-x_{j}|\right),$
(57)

where $x_{i=1,2,3,4}$ are endpoints of non-overlapping subregions
$A=[x_{1},x_{2}]$ and $B=[x_{3},x_{4}]$ present in Fig. 22a.

(a) Critical $\mathrm{sym}S_{A}^{bd}$ scaling of RXPM

(b) Critical $\mathrm{sym}I_{AB}^{bd}$ scaling of RXPM

Figure 23: RXPM boundary criticality data ($p_{5}=0.743$).

###  IV.3 Plaquette models, stabilizer state measurement, and boundary
entanglement structure

This part discusses the connection between RXPM, stabilizer state measurement,
and the measurement-induced phase transition in hybrid Clifford dynamics. We
show that Pauli measurements on a stabilizer state can be understood in terms
of the Plaquette models, and its spin-flip symmetry characterizes the
entanglement structure induced by the measurement. We use the cluster state
defined on the square lattice as an example to illustrate this connection
further. As previously shown in Refs. [23, 29], measuring the bulk qubits on
the square lattice cluster state will induce an interesting entanglement
structure on the boundary. In this section, we further show that the boundary
spin-flip symmetry group of RXPM can fully characterize this entanglement
structure.

####  IV.3.1 Plaquette models and measuring the stabilizer state

The stabilizer state is a special quantum state where all information is
conveniently stored in a stabilizer tableau. It plays a critical role in
simulating Clifford circuits, and an efficient classical algorithm for this is
outlined in [37].

A stabilizer state $\ket{\phi}_{S}$ adheres to the condition:

$\ket{\phi}_{S}=\mathbf{g}\ket{\phi}_{S}\quad\mathbf{g}\in\mathcal{S}\subset\mathcal{P}_{N}.$
(58)

Here, $\mathbf{g}$ is the group element of the stabilizer group $\mathcal{S}$,
and $\mathcal{P}_{N}$ denotes the $N$-qubit Pauli group. Performing a Pauli
measurement on such a quantum state induces a change in the stabilizer group
$\mathcal{S}\to\mathcal{S}^{\prime}$. To elucidate this change, we introduce
the commutator matrix $P$:

$O_{i}\mathbf{g}_{j}=(-1)^{P_{ij}}\mathbf{g}_{j}O_{i}.$ (59)

This $P$ matrix encapsulates the commutation relation between the pre-
measurement stabilizer group and the measurement observable. The post-
measurement stabilizer group $\mathcal{S}^{\prime}$ comprises combinations of
generators from the pre-measurement stabilizer group $\mathcal{S}$ and the
measurement observables $\\{O_{i}\\}$. Formally:

$S^{\prime}=\big{\langle}O_{i},~{}\mathbf{g}^{\prime}_{k}\equiv\prod_{j}\mathbf{g}_{j}^{x_{k}^{j}},~{}\mathbf{g}_{m}\big{\rangle},$
(60)

where $x_{k}^{j}=0,1$ characterizes the stabilizer generators induced by the
measurement and $\mathbf{g}_{m}$ denotes pre-measurement stabilizers that
commute with all the observables, $[\mathbf{g}_{m},O_{i}]=0$ for all $i,m$.
From here and in the remaining part of this section, we use $\prime$ to label
the post-measurement stabilizer group and its elements.

To maintain the abelian nature of the stabilizer group, it is imperative that

$[O_{i},\mathbf{g}^{\prime}_{k}]=0,$ (61)

which can be further written as

$\big{[}1-(-1)^{\sum_{j}P_{ij}x_{k}^{j}}\big{]}O_{i}\mathbf{g}^{\prime}_{k}=0,$
(62)

signifying that

$\sum_{j}P_{ij}x_{k}^{j}=0~{}\operatorname{mod}~{}2\quad\forall i,k.$ (63)

This equation is the same as Eq. 4, which is used to capture the ground state
physics of the plaquette models. Taking advantage of this equation, we may
directly write down the corresponding plaquette model Hamiltonian

$H[\sigma]=-\sum_{q}J_{q}\prod_{i:P_{qi}=1}\sigma_{i},$ (64)

where the product $\prod_{i:P_{qi}=1}$ is over non-trivial entries $i$ of the
$q$-th row of the matrix $P$. The Pauli measurement process is thus related to
the Plaquette models. Each of the post-measurement stabilizer

$\mathbf{g}^{\prime}_{k}\equiv\prod_{j}\mathbf{g}_{j}^{x_{k}^{j}}$ (65)

is thus related to the spin-flip symmetry operator in the plaquette model via
the vector $\vec{x}_{k}=(x_{k}^{1},\ldots,x_{k}^{N})$.

####  IV.3.2 Bulk measurement induced boundary entanglement

We now consider a specific setup where a $d$-dimensional quantum stabilizer
state is obtained via measuring the bulk of a $d+1$-dimensional stabilizer
state. An example for $d=1$ is presented in Fig. 24. This type of system has
been previously studied in [23, 29, 38, 39]. These systems are also referred
to as measurement-based quantum computers (MBQC) [40, 41], as the entanglement
structure of the post-measurement state is controlled by the bulk measurement
direction. We now show that the entanglement structure of such a quantum
system can be understood via the structure of the boundary spin-flip group of
the corresponding plaquette models.

Figure 24: Measuring the bulk of a 2D quantum state on a 2D manifold results
in a 1D quantum state living on the manifold boundary. The measurement-induced
stabilizers separates into two types: 1. Bulk stabilizers
${\mathbf{g}_{k}^{\prime}}^{1}$. 2. Boundary stabilizers
${\mathbf{g}_{k}^{\prime}}^{2}$. The entanglement strcutre of the boundary
state is captured by the boundary stabilizers ${\mathbf{g}_{k}^{\prime}}^{2}$.

As shown in Eq. 60, the post-measurement stabilizer group consists of three
parts: the measurement observables $\\{O_{i}\\}$, the measurement-induced
stabilizer $\mathbf{g}_{k}^{\prime}$. In the MBQC setup, if we further assume
that the bulk measurements are single-site Paulis $O_{i}=X_{i},Y_{i},Z_{i}$,
and no stabilizer lives solely on the boundary, the corresponding stabilizer
group is then

$\mathcal{S}^{\prime}=\langle\mathbf{g}^{\prime}_{k},O_{i}\rangle.$ (66)

The non-trivial entanglement structure of the unmeasured boundary is thus
captured by $\mathbf{g}^{\prime}_{k}$’s. The $\\{\mathbf{g}^{\prime}_{k}\\}$
are naturally separated into two types as present in Fig. 24: 1) the bulk
stabilizers, those support trivially on the boundary labeled by
$\\{{\mathbf{g}^{\prime}_{k}}^{1}\\}$ and 2) the boundary stabilizers, those
have non-trivial support on the boundary labeled by
$\\{{\mathbf{g}^{\prime}_{k}}^{2}\\}$. Obviously, the first type can be
directly decomposed into the bulk measurement observables

${\mathbf{g}^{\prime}_{k}}^{1}=\prod O_{i}.$ (67)

The post-measurement stabilizer group is then

$\mathcal{S}^{\prime}=\langle{\mathbf{g}^{\prime}_{k}}^{2},O_{i}\rangle,$ (68)

and non-trivial boundary entanglement structure is then fully captured by the
boundary part ${\mathbf{g}^{\prime}_{k}}^{2}$.

Formally, the bulk stabilizers form a subgroup of the full stabilizer group
${\mathcal{S}^{1}}^{\prime}=\\{{\mathbf{g}^{\prime}_{k}}^{1}\\}\subset{\mathcal{S}}^{\prime}$,
and the boundary part
${\mathcal{S}^{2}}^{\prime}=\\{{\mathbf{g}^{\prime}_{k}}^{2}\\}$ is the group
quotient

${\mathcal{S}^{2}}^{\prime}=\mathcal{S}^{\prime}/{\mathcal{S}^{1}}^{\prime}.$
(69)

Taking advantage of the relation between the post-measurement stabilizers and
the spin-flip symmetry obtained in Eq. 65, and the bulk-boundary
correspondence of symmetries given in Eq. 38

$G^{bd}=G/G^{bk},$

it is now obvious that the bulk stabilizers $\mathcal{S}_{1}^{\prime}$
correspond to the bulk spin-flip symmetry $G^{bk}$ and the boundary
stabilizers $\mathcal{S}_{2}^{\prime}$ then correspond to the boundary spin-
flip symmetry $G^{bd}$.

Recall that the entanglement entropy of the stabilizer state [42]

$S_{A}=\frac{1}{2}\log|\mathcal{S}_{A\overline{A}}|$ (70)

where $\log|\mathcal{S}_{A\overline{A}}|$ counts the rank group quotient

$\mathcal{S}_{A\overline{A}}\equiv\mathcal{S}/(\mathcal{S}_{A}\mathcal{S}_{\overline{A}})$
(71)

where $\mathcal{S}_{A/\overline{A}}\subset\mathcal{S}$ are the subgroup of
stabilizers that operate only on domain $A/\overline{A}$.

(a) 2D domain: domain $A,A_{0}$ is enclosed by red/green line

(b) 1D boundary domain: $A,A_{0}$ is enclosed by the red/green dashed line

Figure 25:  Relation between domain $A$ and $A_{0}$.

As shown in Eq. 64, the stabilizer group $\mathcal{S}^{\prime}$ and the
plaquette model spin-flip group $G$ are related via the binary vector set
$\\{x_{k}\\}$. As by assumption $g_{i}$’s are local operators, the
entanglement entropy of the boundary state $\ket{\psi}_{\partial M_{2}}$ is
then related to the boundary symmetry entropy $\mathrm{sym}$ up to an area-law
deviation in the following manner

$S_{A}=\frac{1}{2}\text{sym}S_{A}^{bd}+O(l_{A}),$ (72)

where $l_{A}$ is the size of the boundary of domain $A$. For a 1D domain, its
boundary $l_{A}\sim 1$. This deviation $O(l_{A})$ originates from the
difference in the local subgroup size $\log|G_{A}|$ and
$\log|\mathcal{S}^{\prime}_{A}|$, which can be directly calculated as

$\displaystyle\delta_{A}$
$\displaystyle\equiv\log|G_{A}|-\log|\mathcal{S}^{\prime}_{A}|$ (73)
$\displaystyle=\log|G_{A}|-\log|G_{A_{0}}|$ $\displaystyle=\rank
T_{\overline{A_{0}}}-\rank T_{\overline{A}}$ $\displaystyle\leq\rank
T_{A\backslash A_{0}}\leq\min\\{|G|,l_{A}\\}$ $\displaystyle\leq l_{A}$

where in the second line $G_{A_{0}}$ is isomorphic to the local stabilizer
group $\mathcal{S}_{A}^{\prime}$ as generators $g_{i}$’s are local we have
$A_{0}\subset A$ and $|A\backslash A_{0}|\sim l_{A}$, and the geometry of
domain $A$ and $A_{0}$ is presented in Fig. 25. In the third line, we used the
fact that

$\log|G_{\mathcal{D}}|=\rank T-\rank T_{\overline{\mathcal{D}}}$ (74)

for any domain $\mathcal{D}$ with $N$ being the total number of sites. In the
fourth line, we used the rank inequality that for any matrices $\mathbf{A}$
and $\mathbf{B}$

$\rank{\mathbf{A}}+\rank{\mathbf{B}}\geq\rank[\mathbf{A}\mid\mathbf{B}],$ (75)

and in the fifth line for any $a\times b$ matrix $\mathbf{A}$

$\rank\mathbf{A}\leq\min\\{a,b\\}.$ (76)

####  IV.3.3 RXPM, MBQC boundary, and Clifford MIPT

Building upon the mappings described above, we are now ready to show the
equivalence between the boundary symmetry phase transition observed in RXPM
and the measurement induced boundary phase transition on a two dimensional
cluster state.

Figure 26: $Y$ measurement induces $5$ non-zero entries in the corresponding
row of the matrix $P$, corresponding to site $i$ and its four neighboring
sites (highlighted in Green). $Z$ measurement only induces one non-zero entry
result in in the corresponding row of the matrix $P$ (green dot enclosed by
the dashed line).

For the square lattice graph state, the stabilizer generators are given by

$g_{i}=X_{i}\prod_{j\in N_{i}}Z_{j},$ (77)

where $N_{i}$ is the set of four neighboring sites of site $i$. We measure all
of the bulk qubits in $Y$ or $Z$ directions. As present in Fig. 26, the single
site $Y$ measurement anti-commutes with five $g_{i}$ and corresponds to the
five-body interaction. On the other hand, the $Z$ measurement corresponds to
the single-site term, as this single $Z$ only anti-commutes with one $g_{i}$.

The generators of the post-measurement stabilizer state corresponds to the
boundary symmetry operators of RXPM. In addition, the entanglement entropy of
the boundary stabilizer state is equivalent to half of the boundary symmetry
entropy $\mathrm{sym}S^{bd}_{A}$ with respect to subregion $A$ up to a finite
constant.

The critical exponents extracted in the RXPM model (this work), the Y/Z bulk
measurement-induced boundary criticality on a 2D square lattice cluster state
[29], and the measurement-induced phase transition in random Clifford circuits
[15, 21, 43] are shown in Tab. 1. We propose that these criticalities belong
to the same universality class, supported by the proximity of exponents. Here,
$c$ and $h^{0}$ for RXPM and Clifford MIPT/MBQC (square) differ by a prefactor
$2$. The reason lies in the definition of stabilizer entanglement entropy and
symmetry entropy shown in Eq. 72.

Model |  $z$ |  $h^{0}$ |  $h^{1}$ |  $\nu$ |  $c$ |  $\Delta$
---|---|---|---|---|---|---
RXPM |  $1$ |  $1.53$ |  $0.125$ |  $1.3$ |  $3.05$ |  $2$
MBQC(Square) |  $1$ |  - |  - |  $1.3$ |  $1.6$ |  $2$
Cliff. MIPT |  $1$ |  $0.76$ |  $0.125$ |  $1.3$ |  $1.6$ |  $2$
Table 1:  Critical exponents of RXPM boundary, random Y/Z measurement MBQC on
the square-lattice boundary [29] and Clifford MIPT [15, 21, 43]. Here, we drop
the subscript of $\nu_{5}$ for the RXPM model.

##  V Conclusion and discussion

In this work, we introduce a class of classical spin models with $q$-spin
Ising interactions defined on the two dimensional lattice. These interactions
are randomly replaced with single-site terms with probability $1-p$, leading
to what we term random plaquette models (RPMs). We investigate the ground
state phase transitions of these models by analyzing the behavior of the spin-
flip symmetry operators. Specifically, we concentrate on two instances: the
random triangular plaquette model (RTPM) with $q=3$ and the random X-plaquette
model (RXPM) with $q=5$. Our analysis reveals that for $p>p_{c}$, the symmetry
operators become non-local, spanning the entire lattice. Conversely, for
$p<p_{c}$, the symmetry operator localizes within a finite region. To
characterize these localization phase transitions, we develop various tools
tailored to the behavior of these symmetry operators.

We also present a dynamical perspective to elucidate the phase transitions in
the localization of the symmetry operator. We demonstrate that the symmetry
operator, when initiated from a one-dimensional boundary, evolves according to
an update rule determined by the specific interaction form. By employing this
approach, we construct classical 1+1D random cellular automata and illustrate
that the propagation of the symmetry operators, starting from one boundary and
moving into the bulk, can undergo an absorbing phase transition. In
particular, we establish that the RTPM corresponds to a special limit of the
Domany-Kinzel model. Regarding the RXPM, we determine that it possesses a
dynamical exponent $z=1$ at criticality and exhibits an intriguing connection
with the recently discovered measurement-induced entanglement phase transition
(MIPT).

The connection between them is twofold. First, our model offers a classical
counterpart to the MIPT, akin to other proposals discussed in [44, 45, 46,
47]. Here, the single-site constraint plays the role of the measurement in our
model. Second, we establish an equivalence between the boundary symmetry phase
transition of the RXPM and the boundary entanglement phase transitions induced
by the bulk measurement in the cluster state defined on the square lattice.
The latter can be effectively treated as the entanglement phase transition in
a 1+1D hybrid Clifford circuit [15], as previously demonstrated in [29]. By
numerically comparing the critical exponents in these three transitions, we
demonstrate that the phase transitions in these three models belong to the
same universality class, as presented in Tab. 1.

We conclude this work by highlighting several promising directions for further
research. Firstly, the physics explored in this paper can be generalized to
other 2D models with multi-spin interactions or even some non-local spin
models. It would be interesting to explore the universality classes of
potential phase transitions in these spin models.

Secondly, the symmetry operator defined in our classical spin model is
equivalent to the logical operator in classical LDPC codes. Exploring our
transition from the perspective of classical error correction and even
promoting it to the quantum version could be an interesting direction.

Thirdly, our model offers a new theoretical framework for understanding the
MIPT in hybrid Clifford circuits. Generalizing our approach to more complex
classical spin models could deepen our understanding of the MIPT in Haar
random circuits, which is currently only comprehended in specific limits or at
a mean-field level [21, 19, 48].

Fourthly, it is known that the Newman-Moore model with three-body interaction
does not have a thermal phase transition but does exhibit some glassy dynamics
at finite temperature [3]. We expect that substituting the three-body
interaction terms with single-site terms would suppress the glassy behavior at
finite temperature. Therefore it would be interesting to explore the finite
temperature physics of the models discussed in this paper and explore the
potential glass phase transitions in them.

Lastly, our model bears some resemblance to the XOR-SAT problem, where two
transitions have been identified by varying the number of random non-local
constraints [49, 50]. Our model appears to be a local version of XOR-SAT, and
exploring the connections between these two models could be interesting.

###### Acknowledgements.

We gratefully acknowledge computing resources from Research Services at Boston
College and the assistance provided by Wei Qiu. H.L. thanks Ethan Lake,
Shengqi Sang, Yabo Li and Tianci Zhou for insightful discussions. This
research is supported in part by the National Science Foundation under Grant
No. DMR-2219735.

## References

  * McCoy and Wu [1973] B. M. McCoy and T. T. Wu,  _The two-dimensional Ising model_ (Harvard University Press, 1973).
  * Lipowski [1997] A. Lipowski, Glassy behaviour and semi-local invariance in ising model with four-spin interaction, Journal of Physics A: Mathematical and General 30, 7365 (1997).
  * Newman and Moore [1999] M. E. J. Newman and C. Moore, Glassy dynamics and aging in an exactly solvable spin model, Phys. Rev. E 60, 5068 (1999).
  * Nishimori [1979] H. Nishimori, Transition point of the random ising model with bond dilution, Journal of Physics C: Solid State Physics 12, L641 (1979).
  * Barrat [1997] A. Barrat, The p-spin spherical spin glass model, arXiv preprint cond-mat/9701031 (1997).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Rev. Mod. Phys. 51, 659 (1979).
  * Suzuki [1972] M. Suzuki, Solution and critical behavior of some ”three-dimensional” ising models with a four-spin interaction, Phys. Rev. Lett. 28, 507 (1972).
  * Mueller  _et al._ [2017] M. Mueller, D. A. Johnston, and W. Janke, Exact solutions to plaquette ising models with free and periodic boundaries, Nuclear Physics B 914, 388 (2017).
  * You  _et al._ [2018] Y. You, T. Devakul, F. J. Burnell, and S. L. Sondhi, Subsystem symmetry protected topological order, Physical Review B 98, 035112 (2018).
  * Devakul  _et al._ [2019] T. Devakul, Y. You, F. Burnell, and S. L. Sondhi, Fractal symmetric phases of matter, SciPost Physics 6, 007 (2019).
  * Domany and Kinzel [1984] E. Domany and W. Kinzel, Equivalence of cellular automata to ising models and directed percolation, Phys. Rev. Lett. 53, 311 (1984).
  * Kinzel [1985] W. Kinzel, Phase transitions of cellular automata, Zeitschrift für Physik B Condensed Matter 58, 229 (1985).
  * Margolus [1984] N. Margolus, Physics-like models of computation, Physica D: Nonlinear Phenomena 10, 81 (1984).
  * Aharonov [2000] D. Aharonov, Quantum to classical phase transition in noisy quantum computers, Phys. Rev. A 62, 062311 (2000).
  * Li  _et al._ [2019] Y. Li, X. Chen, and M. P. Fisher, Measurement-driven entanglement transition in hybrid quantum circuits, Physical Review B 100, 134306 (2019).
  * Chan  _et al._ [2019] A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Unitary-projective entanglement dynamics, Phys. Rev. B 99, 224307 (2019).
  * Choi  _et al._ [2020] S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Quantum error correction in scrambling dynamics and measurement-induced phase transition, Phys. Rev. Lett. 125, 030505 (2020).
  * Gullans and Huse [2020a] M. J. Gullans and D. A. Huse, Dynamical purification phase transition induced by quantum measurements, Phys. Rev. X 10, 041020 (2020a).
  * Jian  _et al._ [2020] C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. W. Ludwig, Measurement-induced criticality in random quantum circuits, Phys. Rev. B 101, 104302 (2020).
  * Li and Fisher [2021] Y. Li and M. P. A. Fisher, Statistical mechanics of quantum error correcting codes, Phys. Rev. B 103, 104306 (2021).
  * Skinner  _et al._ [2019] B. Skinner, J. Ruhman, and A. Nahum, Measurement-induced phase transitions in the dynamics of entanglement, Physical Review X 9, 031009 (2019).
  * Sang  _et al._ [2021] S. Sang, Y. Li, T. Zhou, X. Chen, T. H. Hsieh, and M. P. Fisher, Entanglement negativity at measurement-induced criticality, PRX Quantum 2, 030313 (2021).
  * Liu  _et al._ [2022] H. Liu, T. Zhou, and X. Chen, Measurement-induced entanglement transition in a two-dimensional shallow circuit, Phys. Rev. B 106, 144311 (2022).
  * Bao  _et al._ [2022] Y. Bao, M. Block, and E. Altman, Finite time teleportation phase transition in random quantum circuits (2022), arXiv:2110.06963 [quant-ph] .
  * Li  _et al._ [2018] Y. Li, X. Chen, and M. P. A. Fisher, Quantum zeno effect and the many-body entanglement transition, Physical Review B 98, 10.1103/physrevb.98.205136 (2018).
  * Vasseur  _et al._ [2019] R. Vasseur, A. C. Potter, Y.-Z. You, and A. W. W. Ludwig, Entanglement transitions from holographic random tensor networks, Phys. Rev. B 100, 134203 (2019).
  * Yang  _et al._ [2022] Z.-C. Yang, Y. Li, M. P. A. Fisher, and X. Chen, Entanglement phase transitions in random stabilizer tensor networks, Phys. Rev. B 105, 104306 (2022).
  * Negari  _et al._ [2023] A.-R. Negari, S. Sahu, and T. H. Hsieh, Measurement-induced phase transitions in the toric code, arXiv preprint arXiv:2307.02292 (2023).
  * Liu  _et al._ [2023] H. Liu, V. Ravindranath, and X. Chen, Quantum entanglement phase transitions and computational complexity: Insights from ising models, arXiv preprint arXiv:2310.01699 (2023).
  * Pless [1998] V. Pless,  _Introduction to the theory of error-correcting codes_ , Vol. 48 (John Wiley & Sons, 1998).
  * Mezard and Montanari [2009] M. Mezard and A. Montanari,  _Information, physics, and computation_ (Oxford University Press, 2009).
  * Schuster  _et al._ [2023] T. Schuster, N. Tantivasadakarn, A. Vishwanath, and N. Y. Yao, A holographic view of topological stabilizer codes, arXiv preprint arXiv:2312.04617 (2023).
  * Yoshida [2013] B. Yoshida, Exotic topological order in fractal spin liquids, Phys. Rev. B 88, 125122 (2013).
  * Wolfram [1983] S. Wolfram, Statistical mechanics of cellular automata, Rev. Mod. Phys. 55, 601 (1983).
  * Kemper  _et al._ [2001] A. Kemper, A. Schadschneider, and J. Zittartz, Transfer-matrix density-matrix renormalization group for stochastic models: the domany-kinzel cellular automaton, Journal of Physics A: Mathematical and General 34, L279 (2001).
  * Atman  _et al._ [2002] A. P. F. Atman, R. Dickman, and J. G. Moreira, Scaling exponents of rough surfaces generated by the domany-kinzel cellular automaton, Phys. Rev. E 66, 016113 (2002).
  * Gottesman [1998] D. Gottesman, The heisenberg representation of quantum computers (1998), arXiv:quant-ph/9807006 [quant-ph] .
  * Okuda  _et al._ [2024] T. Okuda, A. P. Mana, and H. Sukeno, Anomaly inflow, dualities, and quantum simulation of abelian lattice gauge theories induced by measurements (2024), arXiv:2402.08720 [cond-mat.str-el] .
  * Negari  _et al._ [2024] A.-R. Negari, S. Sahu, and T. H. Hsieh, Measurement-induced phase transitions in the toric code, Phys. Rev. B 109, 125148 (2024).
  * Raussendorf  _et al._ [2003] R. Raussendorf, D. E. Browne, and H. J. Briegel, Measurement-based quantum computation on cluster states, Phys. Rev. A 68, 022312 (2003).
  * Raussendorf and Briegel [2001] R. Raussendorf and H. J. Briegel, A one-way quantum computer, Phys. Rev. Lett. 86, 5188 (2001).
  * Fattal  _et al._ [2004] D. Fattal, T. S. Cubitt, Y. Yamamoto, S. Bravyi, and I. L. Chuang, Entanglement in the stabilizer formalism (2004), arXiv:quant-ph/0406168 [quant-ph] .
  * Gullans and Huse [2020b] M. J. Gullans and D. A. Huse, Dynamical purification phase transition induced by quantum measurements, Physical Review X 10, 041020 (2020b).
  * Iaconis  _et al._ [2020] J. Iaconis, A. Lucas, and X. Chen, Measurement-induced phase transitions in quantum automaton circuits, Physical Review B 102, 10.1103/physrevb.102.224311 (2020).
  * Han and Chen [2023] Y. Han and X. Chen, Entanglement structure in the volume-law phase of hybrid quantum automaton circuits, Physical Review B 107, 10.1103/physrevb.107.014306 (2023).
  * Lyons  _et al._ [2023] A. Lyons, S. Choi, and E. Altman, Universal crossover in quantum circuits governed by a proximate classical error correction transition, Physical Review B 107, 10.1103/physrevb.107.174304 (2023).
  * Willsher  _et al._ [2022] J. Willsher, S.-W. Liu, R. Moessner, and J. Knolle, Measurement-induced phase transition in a chaotic classical many-body system, Physical Review B 106, 10.1103/physrevb.106.024305 (2022).
  * Bao  _et al._ [2020] Y. Bao, S. Choi, and E. Altman, Theory of the phase transition in random unitary circuits with measurements, Physical Review B 101, 10.1103/physrevb.101.104301 (2020).
  * Kirkpatrick and Selman [1994] S. Kirkpatrick and B. Selman, Critical behavior in the satisfiability of random boolean expressions, Science 264, 1297 (1994), https://www.science.org/doi/pdf/10.1126/science.264.5163.1297 .
  * Mézard  _et al._ [2003] M. Mézard, F. Ricci-Tersenghi, and R. Zecchina, Two solutions to diluted p-spin models and xorsat problems, Journal of Statistical Physics 111, 505 (2003).


# Extremal ternary self-dual codes of length 36 and symmetric 2-$(36,15,6)$
designs with an automorphism of order 2

Sanja Rukavina  Vladimir D. Tonchev  Department of Mathematical Sciences,
Michigan Technological University, Houghton, MI 49931, USA
E-mail<EMAIL_ADDRESS>ORCID: 0000-0003-1806-3571

###### Abstract

In this note we report the classification of all symmetric 2-$(36,15,6)$
designs that admit an automorphism of order 2 and their incidence matrices
generate an extremal ternary self-dual code. It is shown that up to
isomorphism, there exists only one such design, having a full automorphism
group of order 24, and the ternary code spanned by its incidence matrix is
equivalent to the Pless symmetry code.

* Corresponding author.

Keywords: Pless symmetry code, Hadamard matrix, symmetric 2-design,
automorphism group

Mathematical subject classification (2020): 05B05, 05B20, 94B05

##  1 Introduction

We assume familiarity with the basic facts and notions from error-correcting
codes, combinatorial designs and Hadamard matrices [2], [3], [9], [10], [17].
All codes considered in this paper are ternary.

The minimum distance, or equivalently, the minimum weight $d$ of a ternary
self-dual code of length $n$ divisible by 12 satisfies the upper bound $d\leq
n/4+3$ (cf., e.g. [9, 9.3]), and a self-dual code with minimum distance
$d=n/4+3$ is called extremal. The extremal self-dual codes support
combinatorial 5-designs by the Assmus-Mattson theorem [1], [9, 8.4].

Ternary extremal self-dual codes are known for the following lengths $n\equiv
0\pmod{12}$: $n=12$: the extended Golay code, unique up to equivalence;
$n=24$: there are exactly two inequivalent codes, the extended quadratic-
residue code [1] and the Pless symmetry code $C(11)$ [13], [14]; $n=36$: only
one code is known, namely the Pless symmetry code $C(17)$ [13], [14]; $n=48$:
two codes are known, the extended quadratic-residue code and the Pless
symmetry code $C(23)$; $n=60$: three codes are known: the extended quadratic-
residue code, the Pless symmetry code $C(29)$, and a code found by Nebe and
Villar [12].

Huffman [8] proved that any extremal ternary self-dual code of length 36 that
admits an automorphism of prime order $p>3$ is monomially equivalent to the
Pless symmetry code. More recently, Eisenbarth and Nebe [7] extended Huffman’s
result by proving that the Pless symmetry code is the unique (up to monomial
equivalence) ternary extremal self-dual code of length 36 that admits an
automorphism of order 3. In addition, it was proved in [7, Theorem 5.1] that
if $C$ is an extremal ternary self-dual code of length 36 then either $C$ is
equivalent to the Pless symmetry code or the full automorphism group of $C$ is
a subgroup of the cyclic group of order 8.

It is known [14] that the Pless symmetry code $C(q)$ of length $n=2q+2$, where
$q\equiv-1\pmod{3}$ is an odd prime power, contains a set of $n$ codewords of
weight $n$, which after replacing every entry equal to 2 by $-1$ form the rows
of a Hadamard matrix equivalent to the Paley-Hadamard matrix of type II. In
particular, the Pless symmetry code $C(17)$ contains the rows of a Hadamard
matrix $P$ of Paley type II, having a full automorphism group of order $4\cdot
17(17^{2}-1)=19584$, and the rows of $P$ span the code $C(17)$.

It was shown in [16] that the code $C(17)$ contains a second equivalence class
of Hadamard matrices of order 36 having as rows codewords of $C(17)$. Any
matrix $H$ from the second equivalence class has a full automorphism group of
order 72 and the rows of $H$ span the code $C(17)$. In addition, $H$ is
monomially equivalent to a regular Hadamard matrix $H^{\prime}$, such that
every row of $H^{\prime}$ has 15 entries equal to 1 and 21 entries equal to
$-1$ . The symmetric 2-$(36,15,6)$ design $D$ with $(0,1)$-incidence obtained
by replacing every entry $-1$ of $H^{\prime}$ with 0 has a trivial full
automorphism group, and the row span of its incidence matrix over $GF(3)$ is
equivalent to the Pless symmetry code $C(17)$.

In Section 2, we consider codes that are obtained from $C(17)$ by negating
code coordinates so that the resulting new code contains the all-one vector,
and examine Hadamard matrices formed by codewords of weight 36. This
examination reveals the surprising fact that the Pless symmetry code $C(17)$
is equivalent to a code spanned by the rows of a regular Hadamard matrix which
is monomially equivalent to the Paley-Hadamard matrix of type II, the
associate symmetric 2-$(36,15,6)$ design has a full automorphism group of
order 24, and its incidence matrix spans a code equivalent to $C(17)$.

Motivated by this phenomenon and the results from [7], in Section 3 we
classify all symmetric 2-$(36,15,6)$ designs that admit an automorphism of
order 2 and their incidence matrices span an extremal ternary self-dual code
of length 36. The results of this classification show that, up to isomorphism,
there is only one such design, being isomorphic to the design associated with
the regular presentation of the Paley-Hadamard matrix of type II found in
Section 2.

##  2 Hadamard matrices derived from the symmetry code of length 36

In this section we consider Hadamard matrices of order 36 whose rows are
obtained from codewords of full weight 36 in the ternary Pless symmetry code
$C(17)$ after replacing all codeword entries that are equal to 2 with $-1$.
The code $C(17)$ contains exactly 888 codewords of weight 36. In what follows,
we denote by $U$ the set of all 888 codewords in $C(17)$ of full weight. The
3-rank (that is, the rank over the finite field of order 3, $GF(3)$) of the
matrix having as rows the codewords from $U$ is 18, hence $U$ spans $C(17)$.
In the context of Hadamard matrices, we view the codewords from $U$ as vectors
with components $\pm 1$.

###### Lemma 2.1

A set $S$ of 36 codewords from $U$ is the row set of a Hadamard matrix of
order 36 if and only if the Hamming distance between every two codewords from
$S$ is 18.

Proof. The inner product of two vectors $u=(u_{1},\ldots,u_{36})$,
$v=(v_{1},\ldots,v_{36})$ of length 36 with entries $\pm 1$ is zero over the
field of real numbers if and only if there are exactly 18 indices $i$, $1\leq
i\leq 36$, such that $u_{i}v_{i}=-1$, hence $u_{i}\neq v_{i}$, and the
remaining 18 indices $j$, $1\leq j\leq 36$, satisfy $u_{j}v_{j}=1$, that is,
$u_{j}=v_{j}$. $\Box$

By Lemma 2.1, every Hadamard matrix of order 36 whose rows are codewords from
$U$ corresponds to a clique of size 36 in a graph $\Gamma$ having as vertices
the codewords from $U$, where two codewords are adjacent in $\Gamma$ if they
are at Hamming distance 18 from each other. Enumerating all Hadamard matrices
up to equivalence by a direct search for cliques of size 36 is computationally
infeasible because for every Hadamard matrix $H$, the graph $\Gamma$ contains
$2^{36}$ distinct 36-cliques that correspond to $2^{36}$ Hadamard matrices
that are equivalent to $H$ and are obtained by negating rows of $H$. Note that
negating a row is equivalent to a scalar multiplication of the corresponding
codeword by 2 (mod 3), thus any negation of a row preserves $U$.

A simple way to reduce the search is by considering normalized Hadamard
matrices.

###### Lemma 2.2

If $H$ is a Hadamard matrix whose rows are codewords from $U$ then the row set
of any normalized Hadamard matrix obtained from $H$ belongs to a code which is
monomially equivalent to $C(17)$.

Proof. The matrix $H$ is normalized with respect to column $j$ by negating (or
multiplying by 2 (mod 3)) all rows of $H$ having entry $-1$ in the $j$th
column. Any such negation preserves $U$. Similarly, $H$ is normalized with
respect to row $i$ by negating all columns of $H$ having entry $-1$ in the
$i$th row. These negations of the columns of $U$ either preserve $U$ or change
$U$ to a new set that spans a linear code which is monomially equivalent to
$C(17)$. $\Box$

The set $U$ of the Pless symmetry code $C(17)$, as defined in [13], [14], does
not contain the all-one codeword $\bar{1}=(1,\ldots,1)$, while it contains a
codeword $v$ with one entry equal to $-1$, located in the code coordinate
labeled by $\infty$, and 35 entries equal to 1. Negating (that is, multiplying
by 2 (mod 3)) the code coordinate of $C(17)$ labeled by $\infty$ transforms
the symmetry code $C(17)$ into a monomially equivalent code $L(17)$ which does
contain the all-one vector $\bar{1}$ [16]. If $x\in L(17)$ is a codeword of
full weight 36, we denote by $w_{i}(x)$ the number of entries of $x$ that are
equal to $i$ ($i=1,2$).

The complete weight distribution of the set $W$ consisting of all 888
codewords of $L(17)$ of weight 36 is listed in Table 1. The set $W$ is
available at

     https://pages.mtu.edu/~tonchev/W.txt

As shown in [16], the set $M$ of the 70 codewords $x\in W$ with weight
structure $(w_{1}(x),w_{2}(x))=(18,18)$ form the $(1,2)$-incidence matrix of a
Hadamard 3-$(36,18,8)$ design with full automorphism group of order
$272=2^{4}17$, while adding to $M$ the all-one vector $\bar{1}$ and the all-
two vector $\bar{2}=2\bar{1}$ gives a set of 72 vectors that comprises of the
rows of a normalized Hadamard matrix and its negative, being monomially
equivalent to the Paley-Hadamard matrix of type II, having full monomial
automorphism group of order $19584=2^{7}3^{2}17$.

In addition, it was shown in [16] that the set of 408 codewords $x\in W$ with
weight structure $(w_{1}(x),w_{2}(x))=(15,21)$ contains exactly 272 subsets of
size 36 that are regular Hadamard matrices of order 36. All these 272 regular
Hadamard matrices are pairwise equivalent, and have an automorphism group of
order 72, while the associated symmetric 2-$(36,15,6)$ design has trivial full
automorphism group, and its $(0,1)$-incidence matrix has 3-rank 18, hence
spans the code $L(17)$.

All inequivalent Hadamard matrices of order 36 whose rows are obtained from
codewords in $L(17)$ of weight 36 can be enumerated by using Lemma 2.2 as
follows. We consider $W$ as an $888\times 36$ matrix, and for every integer
$i$, $1\leq i\leq 888$, we define a matrix $W_{i}$, being the matrix obtained
by negating all columns of $W$ which have entry $-1$ in row $i$. Thus, the
$i$th row of $W_{i}$ is the all-one vector $\bar{1}$, and this row must be a
row of every normalized Hadamard matrix consisting of rows of $W_{i}$. We
refer to $W_{i}$ as a matrix obtained by switching of $W$ with respect to row
$i$. To reduce the search further, we consider only normalized Hadamard
matrices with first column being the all-one column.

Next we define a graph $\Gamma_{i}$ having as vertices all rows of $W_{i}$
with first entry 1 and exactly 18 entries equal to 1, where two vertices are
adjacent in $\Gamma_{i}$ if and only if the Hamming distance between the
corresponding rows of $W_{i}$ is 18.

###### Lemma 2.3

(a) Any set of 35 rows of $W_{i}$ that corresponds to a clique of size 35 in
$\Gamma_{i}$, together with the all-one row $\bar{1}$, is the set of rows of a
normalized Hadamard matrix.
(b) The maximum clique size in $\Gamma_{i}$ is 35.

Proof. Part (a) follows from Lemma 2.1. If $H$ is a normalized Hadamard matrix
of order $n=2d$, then deleting the all-one column from $H$ gives a
$2d\times(2d-1)$ matrix whose rows form a code $C$ of length $2d-1$ and
minimum distance $d$ over the alphabet $\\{1,-1\\}$ that meets the Plotkin
upper bound [9, 2.2], [17, 2.1]:

$|C|\leq 2d.$

This proves part (b). $\Box$

#$x$ | $(w_{1}(x),w_{2}(x))$
---|---
1 | (0,36)
408 | (15,21)
70 | (18,18)
408 | (21,15)
1 | (36,0)
Table 1: The complete weight distribution of $W$

One can compute representatives of the equivalence classes of Hadamard
matrices by an examination of the 35-cliques in the graphs $\Gamma_{i}$,
$1\leq i\leq 888$. Since every codeword $x\in W$ with weight structure
$(w_{1}(x),w_{2}(x))=(21,15)$ is the negative of the codeword $2x\in W$ with
weight structure $(15,21)$, it is sufficient to examine the graphs
$\Gamma_{i}$ such that the weight structure of the $i$th row of $W$ is
$(15,21)$ or $(18,18)$. The incidence structure with $(0,1)$-incidence matrix
obtained from $W$ by replacing all entries that are equal to 2 with zero, has
full automorphism group $G$ of order 272. The group $G$ partitions the set of
408 rows of $W$ with weight structure $(15,21)$ into four orbits of lengths
68, 68, 136, 136, and the set of 70 rows with weight structure $(18,18)$ into
four orbits of lengths 2, 17, 17, 34. Therefore, it is sufficient to examine
eight switchings of $W,$ one for each orbit, all together.

An examination of the matrices $W_{i}$ obtained by switching of $W$ with
respect to a row $i$ having weight structure $(18,18)$ shows that any such
matrix has the same complete weight distribution as $W$ (see Table 1), and the
resulting unique (up to a permutations of rows) normalized Hadamard matrix is
equivalent to the Paley-Hadamard matrix of type II.

An examination of the matrices $W_{i}$ obtained by switching of $W$ with
respect to any of the 408 rows having weight structure $(15,21)$ shows that
any such matrix has a complete weight distribution given in Table 2, where
$W_{408}$ is obtained by the switching of W with respect to row no. 408.

#$x$ | $(w_{1}(x),w_{2}(x))$
---|---
1 | (0,36)
93 | (12,24)
36 | (15,21)
628 | (18,18)
36 | (21,15)
93 | (24,12)
1 | (36,0)
Table 2: The complete weight distribution of $W_{408}$

Any graph $\Gamma_{i}$ associated with a matrix $W_{i}$ obtained by switching
of $W$ with respect to a row $i$ with weight structure $(15,21)$ has 314
vertices and contains exactly 24 cliques of size 35. The resulting 24
normalized Hadamard matrices are all equivalent to the Hadamard matrix with
full automorphism group of order 72 that was discovered in [16] as a regular
Hadamard matrix derived from the code $L(17)$, whose related symmetric
2-$(36,15,6)$ design has trivial automorphism group.

A surprising result of the switching of $W$ with respect to a row with weight
structure $(15,21)$ that was suggested by the weight distribution in Table 2
and was proved by inspection, is the following.

###### Theorem 2.4

(a) The 36 codewords of $W_{408}$ with weight structure $(15,21)$ form a
regular Hadamard matrix $H$ which is monomially equivalent to the Paley-
Hadamard matrix of type II.

(b) The symmetric 2-$(36,15,6)$ design $D$ associated with $H$ has a full
automorphism group of order 24.

(c) The incidence matrix of $D$ has 3-rank 18, and its linear span over
$GF(3)$ is a code equivalent to the Pless symmetry code $C(17)$.

###### Note 1

Every row of the regular Hadamard matrix $H$ from Theorem 2.4 contains 15
entries equal to 1 and 21 entries equal to $-1$. A $(0,1)$-incidence matrix
$A$ of the associated symmetric 2-$(36,15,6)$ design $D$ is obtained by adding
the all-one codeword $\bar{1}$ to every row of $H$, followed by a
multiplication of all rows by $2\pmod{3}$. Hence, the ternary code spanned by
the rows of $H$ contains also the rows of $A$.

A $(0,1)$-incidence matrix of $D$ is given in the Appendix.

##  3 Symmetric 2-$(36,15,6)$ designs with an involution and their ternary
codes

It was proved in [7, Theorem 5.1] that if $C$ is an extremal ternary self-dual
code of length 36 then either $C$ is equivalent to the Pless symmetry code or
the full automorphism group of $C$ is a subgroup of the cyclic group of order
8. This result and the fact that the Pless symmetry code is equivalent to
codes spanned by the incidence matrices of two nonisomorphic 2-$(36,15,6)$
designs, motivated our study of symmetric 2-$(36,15,6)$ designs with an
automorphism of order 2 (or, an involution), and the ternary codes spanned by
their incidence matrices.

In what follows we summarize the classification of symmetric 2-$(36,15,6)$
designs that are invariant under an involution and the ternary code of length
36 spanned by their incidence matrix is a self-dual code with minimum weight
12.

It is known that every finite group of order $2^{n}$ contains a subgroup of
order $2^{i}$ for every $i$ in the range $1\leq i\leq n$ (cf. [15, Theorem
6.5, page 116]). By this property, finding all nonisomorphic $2$-$(36,15,6)$
designs which are invariant under an involution and whose incidence matrix
spans a ternary extremal self-dual code, will also give the enumeration of all
nonisomorphic $2$-$(36,15,6)$ designs invariant under a nontrivial subgroup of
the cyclic group of order $8$ with this property.

For the construction of $2$-$(36,15,6)$ designs we use the method for
constructing orbit matrices with presumed action of an automorphism group,
which are then indexed to construct designs (see, e.g., [5], [6]). After
constructing $2$-$(36,15,6)$ designs, we check the $3$-rank of their incidence
matrices, and if it is equal to $18$ we determine the minimum weight of the
corresponding ternary code. In our work we use computers. In addition to our
own computer programs, we use computer programs by V. Ćepulić for the
construction of orbit matrices and the computer algebra system MAGMA [4] when
working with codes.

The first step in the construction of $2$-$(36,15,6)$ designs that admit an
involution is to determine the possible orbit lengths distributions. For that
we use the following statements.

###### Theorem 3.1

[11, Corollary 3.7] Suppose that a nonidentity automorphism $\sigma$ of a
symmetric $2$-$(v,k,\lambda)$ design fixes $f$ points. Then

$f\leq v-2(k-\lambda)\qquad{\rm and}\qquad
f\leq(\frac{\lambda}{k-\sqrt{k-\lambda}})v.$

###### Theorem 3.2

[11, Proposition 4.23] Suppose that ${\mathcal{D}}$ is a nontrivial symmetric
$2$-$(v,k,\lambda)$ design, with an involution $\sigma$fixing $f$ points and
blocks. If $f\neq 0$, then

$f\geq\left\\{\begin{tabular}[]{l l}$1+\frac{k}{\lambda}$,&if $k$ and
$\lambda$ are both even,\\\ $1+\frac{k-1}{\lambda}$,&otherwise.\\\
\end{tabular}\right.$

It follows that an involution acting on $2$-$(36,15,6)$ design could have $f$
fixed points, where $f\in\\{4,6,8,10,12,14,16,18\\}$.

Our analysis shows that orbit matrices do not exist for $f\in\\{6,14,18\\}$.
The results for the remaining cases are summarized in Table 3 and the
constructed orbit matrices are available at

     https://www.math.uniri.hr/~sanjar/structures/

The third row of Table 3 contains information on the maximal minimum weight
among the self-dual ternary codes spanned by incidence matrices of the
corresponding $2$-$(36,15,6)$ designs, where symbol $\times$ indicates that
all constructed designs have incidence matrices with 3-rank smaller than 18.

Number of fixed points | 4 | 8 | 10 | 12 | 16
---|---|---|---|---|---
Number of orbit matrices | 12991 | 670 | 56 | 311 | 83
max $d$ in self-dual codes | 12 | $\times$ | $\times$ | 9 | $\times$
Table 3: Symmetric 2-$(36,15,6)$ designs with an involution

All $2$-$(36,15,6)$ designs whose incidence matrix span an extremal ternary
code are mutually isomorphic. They are also isomorphic to the design discussed
in Section 2, whose incidence matrix is given in the Appendix.

The results of our classification of symmetric 2-$(36,15,6)$ designs with the
desired properties can be summarized as follows.

###### Theorem 3.3

(a) Up to isomorphism, there exists exactly one symmetric 2-$(36,15,6)$ design
$D$ that admits an automorphism of order 2 and its incidence matrix spans an
extremal ternary self-dual code of length 36.

(b) The full automorphism group $G$ of $D$ is of order 24, and $G$ is
isomorphic to the symmetric group $S_{4}$.

(c) The regular Hadamard matrix associated with $D$ is equivalent to the
Paley-Hadamard matrix of type II.

(d) The ternary code spanned by $D$ is equivalent to the Pless symmetry code.

Statements and Declarations

Competing Interests: The authors declare no conflict of interest.
Contributions: This is a joint collaboration with both authors contributing
substantially throughout.
Data Availability: The links to data generated and analysed are included in
this article.
Funding: The first author is supported by Croatian Science Foundation under
the project 6732.

## References

  * [1] E. F. Assmus, Jr., H. F. Mattson, Jr., New 5-designs, J. Combin. Theory, Ser. A 6 (1969), 122-151.
  * [2] T. Beth, D. Jungnickel, H. Lenz, Design Theory, 2nd Edition, Cambridge University Press, 1999.
  * [3] A. Betten, M. Braun, H. Fripertinger, A. Kerber, A. Kohnert, A. Wassermann, Error-Correcting Linear Codes: Classification by Isometry and Applications, Springer, Berlin (2006).
  * [4] W. Bosma, J. Cannon, Handbook of Magma Functions, Department of Mathematics, University of Sydney, 1994. http://magma.maths.usyd.edu.au/magma
  * [5] D. Crnković, S. Rukavina, Construction of block designs admitting an abelian automorphism group, Metrika, 62, no. 2-3 (2005), 175–183.
  * [6] V. Ćepulić, On symmetric block designs (45,12,3) with automorphisms of order 5, Ars Combin. 37 (1994), 33–48.
  * [7] S. Eisenbarth, G. Nebe, Self-dual codes over chain rings, Math. Comput. Sci. 14 (2020), 443 - 456.
  * [8] W. C. Huffman, On extremal self-dual ternary codes of lengths 28 to 40, IEEE Trans. Info. Theory 38 No. 4 (1992), 1395-1400.
  * [9] W. C. Huffman, V. Pless, Fundamentals of error-correcting codes, Cambridge University Press, Cambridge 2003.
  * [10] Y. J. Ionin, M. S. Shrikhande, Combinatorics of Symmetric Designs, Cambridge University Press, 2006.
  * [11] E. Lander, Symmetric Designs: An Algebraic Approach, Cambridge University Press, 1983.
  * [12] G. Nebe, D. Villar, An analogue of the Pless symmetry codes, in: Seventh International Workshop on Optimal Codes and Related Topics, Albena, Bulgaria, pp. 158 - 163 (2013).
  * [13] V. Pless, On a new family of symmetry codes and related new five-designs, Bull. Amer. Math. Soc. 75, No. 6 (1969), 1339-1342.
  * [14] V. Pless, Symmetry codes over $GF(3)$ and new five-designs, J. Combin. Theory, Ser. A 12 (1972), 119-142.
  * [15] H. E. Rose, A Course on Finite Groups, Springer-Verlag, London, 2009.
  * [16] V. D. Tonchev, On Pless symmetry codes, ternary QR codes, and related Hadamard matrices and designs, Des. Codes Cryptogr., to appear, doi.org/10.1007/s10623-021-00941-0.
  * [17] V. D. Tonchev, Combinatorial Configurations, John Wiley & Sons, Inc., New York 1988.

##  4 Appendix

     1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
     1 1 0 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
     1 0 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0
     0 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
     1 1 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 0 1 1 1 0 1 0 0 0
     1 1 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 1 0 1 1 1 0 1 0 1 0 0
     1 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 0 1 1
     1 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 1 1 1 0 0 1 0 0 1 1 1
     1 0 1 0 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 0 1 0 1 0 0 0 1 0 0 1 1 1 0 1 0
     1 0 1 0 0 1 0 0 0 1 0 0 1 1 0 0 1 0 1 0 1 0 1 0 0 0 1 0 0 0 1 1 0 1 0 1
     1 0 1 0 0 0 1 0 0 0 0 1 0 0 1 1 0 1 1 0 0 1 1 0 1 0 0 0 1 0 1 0 1 1 0 0
     1 0 1 0 0 0 0 1 0 0 1 0 0 0 1 1 1 0 0 1 1 0 0 1 0 1 0 0 0 1 0 1 1 1 0 0
     1 0 0 1 1 0 0 0 1 0 1 0 0 1 0 1 0 0 1 0 1 1 0 0 0 1 0 0 1 0 0 0 0 1 1 1
     1 0 0 1 0 1 0 0 0 1 0 1 1 0 1 0 0 0 0 1 1 1 0 0 1 0 0 0 0 1 0 0 1 0 1 1
     1 0 0 1 0 0 1 0 1 0 0 1 0 1 1 0 1 0 0 0 0 0 1 1 0 0 0 1 1 1 0 1 0 0 0 1
     1 0 0 1 0 0 0 1 0 1 1 0 1 0 0 1 0 1 0 0 0 0 1 1 0 0 1 0 1 1 1 0 0 0 1 0
     0 1 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 1 0 1 1 0 1 0 0 0 0 1 1 1 0 0 0 1 0 1
     0 1 1 0 0 1 0 0 0 0 1 1 0 1 0 0 1 0 1 0 0 1 0 1 0 0 1 0 1 1 0 0 1 0 1 0
     0 1 1 0 0 0 1 0 1 1 0 0 0 0 0 1 0 1 1 0 1 0 0 1 1 0 0 0 0 1 0 1 0 0 1 1
     0 1 1 0 0 0 0 1 1 1 0 0 0 0 1 0 1 0 0 1 0 1 1 0 0 1 0 0 1 0 1 0 0 0 1 1
     0 1 0 1 1 0 0 0 0 1 0 1 0 1 0 1 0 0 1 1 0 0 1 0 0 1 0 0 0 1 1 1 1 0 0 0
     0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 0 0 0 1 1 0 0 0 1 1 0 0 0 1 0 1 1 0 1 0 0
     0 1 0 1 0 0 1 0 0 1 1 0 0 1 1 0 1 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 1 1 0
     0 1 0 1 0 0 0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 1 1 0 1
     0 0 1 1 1 1 1 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 0 0 1
     0 0 1 1 1 1 0 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 1 1 0 1 1 0 0 1 0
     0 0 1 1 1 0 1 1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 1 0 1 0 1 0 0 0 1 1 0 1
     0 0 1 1 0 1 1 1 1 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 1 0 0 1 1 1 0
     0 0 0 0 1 1 0 0 1 0 0 0 0 1 1 1 1 1 0 1 0 0 1 0 1 0 1 0 0 1 0 0 0 1 1 0
     0 0 0 0 1 1 0 0 0 1 0 0 1 0 1 1 1 1 1 0 0 0 0 1 0 1 0 1 1 0 0 0 1 0 0 1
     0 0 0 0 1 0 0 1 1 1 0 1 0 0 1 0 0 0 1 0 1 1 0 1 0 0 1 1 0 1 1 0 0 1 0 0
     0 0 0 0 0 1 1 0 1 1 1 0 0 0 0 1 0 0 0 1 1 1 1 0 0 0 1 1 1 0 0 1 1 0 0 0
     0 0 0 0 1 0 1 0 0 1 1 1 1 0 0 0 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 0 1 1 0
     0 0 0 0 0 1 0 1 1 0 1 1 0 1 0 0 0 1 0 0 1 0 1 1 1 1 0 0 0 0 1 0 1 0 0 1
     0 0 0 0 0 0 1 1 0 0 1 0 1 1 0 1 1 0 1 1 0 1 0 0 1 0 0 1 0 1 1 0 0 0 0 1
     0 0 0 0 0 0 1 1 0 0 0 1 1 1 1 0 0 1 1 1 1 0 0 0 0 1 1 0 1 0 0 1 0 0 1 0

A symmetric 2-$(36,15,6)$ design associated with the Paley-Hadamard matrix of
type II