∇

# A Multicritical Point with Infinite Fractal Symmetries

Nayan Myerson-Jain, Kaixiang Su, and Cenke Xu  Department of Physics,
University of California, Santa Barbara, CA 93106, USA

###### Abstract

Recently a “Pascal’s triangle model” constructed with ${\rm U}(1)$ rotor
degrees of freedom was introduced, and it was shown that (i.) this model
possesses an infinite series of fractal symmetries; and (ii.) it is the parent
model of a series of $Z_{p}$ fractal models each with its own distinct fractal
symmetry. In this work we discuss a multi-critical point of the Pascal’s
triangle model that is analogous to the Rokhsar-Kivelson (RK) point of the
better known quantum dimer model. We demonstrate that the expectation value of
the characteristic operator of each fractal symmetry at this multi-critical
point decays as a power-law of space, and this multi-critical point is shared
by the family of descendent $Z_{p}$ fractal models. Afterwards, we generalize
our discussion to a $(3+1)d$ model termed the “Pascal’s tetrahedron model”
that has both planar and fractal subsystem symmetries. We also establish a
connection between the Pascal’s tetrahedron model and the ${\rm U}(1)$ Haah’s
code.

##  I Introduction

In recent years we have seen a very fruitful exploration of generalizing the
concept of “symmetry” in physics. It was demonstrated that introducing the
concept of higher form symmetries is very useful for studying gauge fields
Nussinov and Ortiz (2009); Aharony et al. (2013); Gukov and Kapustin (2013);
Kapustin and Thorngren (2013a, b); Kapustin and Seiberg (2014); Gaiotto et al.
(2015); Hsin et al. (2019); Seiberg (2020); McGreevy (2022); in particular the
analysis of the ’t Hooft anomaly of higher form symmetries can lead to very
clear and nontrivial conclusion for the low energy physics of gauge fields
with a topological term Gaiotto et al. (2017). Another set of generalized
symmetries have their conserved charges defined on a subset of the system
Nandkishore and Hermele (2019); Pretko et al. (2020), and these “subsystem
symmetries” can be further categorized into type-I, when the symmetric charges
are defined on regular subsystems such as lines and planes, and type-II, when
the symmetric charges reside on a fractal subsystem with often non-integer
Hausdorff dimensions. So far systems with type-I subsystem symmetries are
better understood, as not only can one study the lattice models with type-I
subsystem symmetries, for gapless states (or phases) with type-I subsystem
symmetry, it is also possible to describe the states with low energy effective
theories which are constructed only with the low energy subspace of the
Hilbert space rather than all the microscopic degrees of freedom of the
system. One classic example of such was constructed by Paramekanti, Balents,
Fisher in Ref. Paramekanti et al., 2002 (dubbed the XY-plaquette model). The
authors constructed a $(2+1)d$ quantum rotor model with one dimensional
conservation of rotor numbers, and the authors studied in great detail a
gapless state of this model. A low energy effective theory constructed with
low energy modes of the system was developed for this gapless state, and a lot
of physics can be studied systematically with the low energy effective theory.
It was demonstrated that standard methods such as renormalization group (RG)
can still be applied to the gapless state of the XY-plaquette model despite
its highly unusual symmetries. The XY-plaquette model along with its low
energy effective theory was later generalized to various models in higher
dimensions, and it has attracted much enthusiasm recently Xu and Fisher
(2007); Ma and Pretko (2018); You et al. (2020); Yuan et al. (2020); Gromov et
al. (2020); Seiberg and Shao (2021); Rudelius et al. (2021); Gorantla et al.
(2021); Burnell et al. (2021).

By the same standard, systems with type-II subsystem symmetries (or fractal
symmetries), are much less understood. Much success was accomplished through
studying exactly soluble lattice models such as the Haah’s code Haah (2011),
but a low energy effective theory which does not involve all the microscopic
degrees of freedom is more challenging compared with systems with type-I
subsystem symmetries. Of course, a low energy effective theory may not be
necessary if one only studies gapped states of type-II models, as under
renormalization group (RG) flow gapped states are expected to “flow to” an
exactly soluble lattice model as the fixed point of the entire gapped phase,
although a proper procedure of RG has not been developed for systems with
fractal symmetries. But recent numerics suggests that a continuous quantum
phase transition may exist for models with fractal symmetries Zhou et al.
(2021) (we note that earlier numerics suggested the opposite Vasiloiu et al.
(2020)), this calls for an understanding of gapless states with fractal
symmetries at the same level as regular systems. A low energy effective theory
will go a long way towards understanding the universal physics of such gapless
systems without relying too much on the lattice models.

Figure 1: The action of the Pascal’s triangle symmetry on a patch (Eq. 4)
where the triangular region outlined in blue is $\mathcal{R}(p)$. In the semi-
classical version of the model, three fractonlike excitations localized at the
triangular plaquettes at the corners of the Pascal’s triangle are created (the
top plaquette shown in red). The lattice sites are indexed by $(i_{1},i_{2})$
with lattice vectors $\mathbf{a}_{1}=\hat{x}$ and
$\mathbf{a}_{2}=\hat{x}/2+\sqrt{3}\hat{y}/2$.

In this work we construct a gapless state with an infinite series of fractal
symmetries. The model we will start with was called the “Pascal’s triangle
model”:

$\displaystyle
H=-K\sum_{\bigtriangledown}\cos(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})+\frac{U}{2}\sum_{i}\hat{n}_{i}^{2},$
(1)

where the boson phase $\hat{\theta}_{i}$ and boson number $\hat{n}_{i}$ are a
pair of conjugate variables which obey the usual algebra
$[\hat{\theta}_{i},\hat{n}_{j}]=\mathrm{i}\delta_{ij}$. This Hamiltonian
exhibits a conventional ${\rm U}(1)\times{\rm U}(1)$ global symmetry involving
rotating $\hat{\theta}$ on any two of the three sublattices of the triangular
lattice by an opposite angle. This model was introduced for the purpose of
understanding the quantum phase transition of the Sierpinski triangle model
Zhou et al. (2021), but it was pointed out in Ref. Myerson-Jain et al., 2022
that this model actually supports a much larger group of fractal symmetries.
In fact, the model possesses an infinite collection of distinct $Z_{p}$
fractal symmetries (one for each prime integer $p$) which together are called
the “Pascal’s triangle symmetry” that involve a staggered rotation of the
boson phase with support on a Pascal’s triangle mod $p$. The action of this
symmetry on a local patch of the lattice is generated by the so-called “patch
symmetry operator” of the Pascal’s triangle symmetry

$\displaystyle U_{p}$ $\displaystyle=$ $\displaystyle
e^{2\pi\mathrm{i}Q_{P}/p}$ (2) $\displaystyle Q_{p}$ $\displaystyle=$
$\displaystyle\sum_{i\in\mathcal{R}(p)}(-1)^{-i_{2}}{-i_{2}\choose
i_{1}}\hat{n}_{i}.$ (4)

The sum in the patch charge $Q_{p}$ is over a triangular region
$\mathcal{R}(p)$ of length $L=p^{k}-1$ (centered at the origin,
$\mathcal{R}(p)\equiv\\{(i_{1},i_{2}):-L\leq i_{2}\leq 0,0\leq i_{1}\leq-
i_{2}\\})$. The action of the patch symmetry operator is shown in Fig. 1. It
was demonstrated in Ref. Myerson-Jain et al., 2022 that for system sizes of
$L=p^{k}-1$, the $Z_{p}$ fractal symmetry becomes exact. The Pascal’s triangle
model also reduces to a family of descendent $Z_{p}$-clock model Hamiltonians
that exhibit only a single one of its fractal symmetries. For example if one
reduces the ${\rm U}(1)$ rotor degree of freedom to an Ising spin, the
Pascal’s triangle model reduces to the Sierpinski triangle model Newman and
Moore (1999); Yoshida (2013).

##  II Multi-critical point of the Pascal’s triangle model

###  II.1 Dual of the Pascal’s triangle model

The dual variables were introduced to understand the physics of the Pascal’s
triangle model Myerson-Jain et al. (2022) Eq. 1:

$\displaystyle\sum_{j\in\bigtriangledown}\hat{\theta}_{j}=\hat{\phi}_{\bar{j}},\
\ \ \ \hat{n}_{j}=-\sum_{\bar{j}\in\bar{\triangle}\ \mathrm{around}\
j}\hat{\Psi}_{\bar{j}}.$ (5)

where $j,\bar{j}$ label the sites of the original and dual lattices
respectively. $\hat{\Psi}$ is a discrete variable and $\hat{\phi}$ is its
conjugate compact variable that satisfies the bosonic algebra
$[\hat{\Psi}_{\bar{i}},\hat{\phi}_{\bar{j}}]=\mathrm{i}\delta_{ij}$. The
Hamiltonian Eq. 1 is now dual to a height model:

$H_{h}=-K\sum_{\bar{j}}\cos(\hat{\phi}_{\bar{j}})+\frac{U}{2}\sum_{\bar{\triangle}}(\hat{\Psi}_{1}+\hat{\Psi}_{2}+\hat{\Psi}_{3})^{2}.$
(6)

We analyze this dual height model following the standard procedure. We first
expand the height model to the quadratic order and add a vertex operator that
re-enforces the compactness of $\hat{\phi}$ (or discreteness of $\hat{\Psi}$),
then $\hat{\Psi}$ can be treated as a continuous variable:

$\displaystyle H_{h}^{g}\sim
H_{h}-\alpha\sum_{\bar{j}}\cos(2\pi\hat{\Psi}_{\bar{j}}).$ (7)

A Euclidean action describing this model can be attained by expanding the
Gaussian part of $H_{h}^{g}$ about the minima of the spectrum which occur at
the $\pm\mathbf{K}=\pm(\frac{4\pi}{3},0)$ points in the Brillouin zone (BZ).
After expanding the height operator

$\displaystyle\Psi(\mathbf{r})\sim
e^{i\mathbf{K}\cdot\mathbf{r}}\psi(\mathbf{r})+e^{-i\mathbf{K}\cdot\mathbf{r}}\psi^{\ast}(\mathbf{r})$
(8)

and taking the continuum limit, the low energy physics of the height model is
supposedly described by the Lagrangian

$\mathcal{L}=\frac{1}{2}(\partial_{\tau}\vec{\psi})^{2}+\frac{\rho_{2}}{2}(\nabla\vec{\psi})^{2}-\alpha\sum_{a}\cos(\vec{e}_{a}\cdot\vec{\psi}).$
(9)

To avoid potential confusion, we clarify that throughout this paper, the
phrase continuum limit simply means the limit where the system size is taken
to infinity while keeping a constant lattice spacing. Eq. 9 becomes an
effective description of the lattice model when the field $\psi$ is a smooth
function on the lattice. The two component vector field $\vec{\psi}$ is
$\vec{\psi}=(\Re{\psi},\Im{\psi})$ and $a=A,B,C$ label the three sublattices
of the dual triangular lattice with $e_{A}=2\pi(1,0)$,
$e_{B}=2\pi(-1/2,\sqrt{3}/2)$, $e_{C}=2\pi(-1/2,-\sqrt{3}/2)$. Analogous to
the analysis of instantons in compact $\text{QED}_{3}$, due to the linear
dispersion modes of the height field, the coefficient of the vertex operators
$\alpha$ is strongly relevant. As a consequence, we conclude that the gapless
modes of the Gaussian theory of the height model $H^{g}_{h}$ are unstable; or
in other words the generic Pascal’s triangle model should in general be in a
gapped quantum disordered phase, just like the compact QED3 should be in a
confined phase without fine-tuning. The connection between the Pascal’s
triangle model and QED3 will be made more explicit later.

###  II.2 Multi-critical point

Although the vertex operators (the $\alpha$-terms) in Eq. 9 destroy the
gapless modes of the Gaussian theory of the dual height model, they can be
suppressed to restore the gaplessness of the spectrum by fine-tuning the model
to a multi-critical point. This multi-critical point has a simple description
in terms of the dual height field:

$\mathcal{L}=\frac{1}{2}(\partial_{\tau}\vec{\psi})^{2}+\frac{\rho_{4}}{2}(\nabla^{2}\vec{\psi})^{2}-\alpha\sum_{a}\cos(\vec{e}_{a}\cdot\vec{\psi}).$
(10)

This Lagrangian is two copies of the $z=2$ quantum Lifshitz theory, which is
pertinent to the study of the RK point of the quantum dimer model in $(2+1)d$
Rokhsar and Kivelson (1988); FRADKIN and KIVELSON (1990); Fradkin et al.
(2004); Ardonne et al. (2004); Vishwanath et al. (2004). This theory is the
$(2+1)d$ analogue of the massless boson CFT in $(1+1)d$ and is known to be
stable against vertex operators for some range of $\rho_{4}$. In particular,
the $\beta$-function for the vertex operators to the leading order is

$\beta(\alpha)=\bigg{(}4-\frac{\pi}{2\sqrt{\rho_{4}}}\bigg{)}\alpha.$ (11)

Later we will further demonstrate that this multi-critical point is shared by
all the descendent $Z_{p}$ fractal models of the Pascal’s triangle model. We
can conclude that the vertex operators are irrelevant when
$\rho_{4}<\pi^{2}/64$ Fradkin et al. (2004).

Returning to lattice model Eq. 1, the multi-critical point can be reached by
turning on nearest neighbor density-density interactions

$\displaystyle\frac{J}{2}\sum_{\langle i,j\rangle}\hat{n}_{i}\hat{n}_{j}\sim$
(12) (13)
$\displaystyle\frac{J}{2}\sum_{\langle\bar{\triangle},\bar{\triangle}^{\prime}\rangle}(\hat{\Psi}_{1}+\hat{\Psi}_{2}+\hat{\Psi}_{3})(\hat{\Psi}_{1^{\prime}}+\hat{\Psi}_{2^{\prime}}+\hat{\Psi}_{3^{\prime}})$
(14)

which changes the band structure of the dual height model. A band structure
calculation shows that by fine-tuning $J$ to certain point, the dispersion of
the height field becomes quadratic.

However, we must be careful since there are possible relevant operators and
other marginal operators that may appear in RG that must be tuned to zero.
These operators must obey the point-group symmetries of the lattice
Hamiltonian. The scaling dimensions of space-time coordinates, and operators
in the Lifshitz theory are $[\partial_{x}]=[\partial_{y}]=1$,
$[\partial_{\tau}]=2$ and $[\psi_{1}]=[\psi_{2}]=0$. An example of a relevant
operator that must be tuned to zero is

$\displaystyle\mathcal{L}_{R}=\alpha^{\prime}\partial_{y}\psi\bigg{(}\frac{\sqrt{3}}{2}\partial_{x}\psi+\frac{1}{2}\partial_{y}\psi\bigg{)}\bigg{(}\frac{\sqrt{3}}{2}\partial_{x}\psi-\frac{1}{2}\partial_{y}\psi\bigg{)}+\text{c.c.}$
(15)

This operator is similar to the one discussed in the context of quantum dimer
model on the honeycomb lattice Fradkin et al. (2004). Several other operators
such as $\alpha^{\prime\prime}(\nabla\psi)^{2}(\nabla\psi^{\ast})^{2}$ are
marginal at the Lifshitz theory by power-counting, but they may become
marginally relevant after considering quantum fluctuations Vishwanath et al.
(2004). These perturbations had better also be tuned to zero to ensure the
system stays at the multi-critical point.

###  II.3 Characteristic features at the multi-critical point

To probe the physics at the multi-critical point, we would like to compute the
expectation value of some characteristic operators. In particular we need an
operator whose expectation value can distinguish different phases of a quantum
many-body system. For a system with a symmetry, the so called “patch operator”
defined as symmetry transformation operation localized on a patch can serve as
a diagnosis of the phase Ji and Wen (2020). This operator was also referred to
as the “disorder operator” for the example of Ising model in previous
literature Kadanoff and Ceva (1971), and later generalized to systems with
type-I subsystem symmetries Wu et al. (2021a) as the order diagnosis operator.
The application of the patch operator can be exemplified in the $1d$ quantum
Ising model: the patch operator has short range, power-law, and long range
expectation value in the Ising ordered phase, critical point, and the Ising
disordered phase respectively (the exact opposite from the correlation
function of the Ising order parameter). For more general cases, the
expectation value of the patch operator should have a qualitative enhancement
in a quantum disordered phase compared with the critical point. In the last
few years the patch operators have been evaluated for various systems, it was
shown that useful information can be extracted from the expectation value of
the patch operator Zhao et al. (2021); Wu et al. (2021b); Wang et al. (2021a,
b).

Here we compute the expectation value of the patch operator of the Pascal’s
triangle symmetry Eq. 4. As explained, opposed to the more standard order
parameter, a long-range expectation value signifies a disordered phase. Power-
law decay, however, still identifies a critical point. In terms of the dual
variables of the height model, the Pascal’s triangle symmetry charge on a
patch $\mathcal{R}(p)$ of length $r=p^{k}$ with the left corner at the origin
becomes

$\displaystyle
Q_{p}(r)=\hat{\Psi}_{0,p^{k}}+(-1)^{p^{k}-1}\sum_{\bar{i}_{1}=0}^{p^{k}}{p^{k}\choose\bar{i}_{1}}\hat{\Psi}_{\bar{i}_{1},0}$
(16)

where the support of the charge on a patch completely cancels in the bulk of
the patch in the dual variables, and the remaining support is on the top
corner and bottom row of the triangular region. However, since
${p^{k}\choose\bar{i}_{1}}$ is always an integer multiple of $p$ for
$\bar{i}_{1}\neq 0,p^{k}$, the patch symmetry operator simplifies even further
to involve only three $\hat{\Psi}$ operators at the corners of the patch:

$\displaystyle
U_{p}(r)=\exp\left(\frac{2\pi\mathrm{i}}{p}(\hat{\Psi}_{0,p^{k}}+(-1)^{p^{k}-1}(\hat{\Psi}_{0,0}+\hat{\Psi}_{p^{k},0}))\right).$
(17)

Let us use lattice vectors of the triangular lattice $\mathbf{a}_{1}=\hat{x}$
and $\mathbf{a}_{2}=\hat{x}/2+\sqrt{3}\hat{y}/2$ where $\hat{x}$ and $\hat{y}$
are the usual Cartesian unit vectors. $\hat{\Psi}_{0,r=p^{k}}$ is at the top
corner of the triangle in the dual representation, while $\hat{\Psi}_{0,0}$
and $\hat{\Psi}_{r=p^{k},0}$ are located at the bottom left and right corners
of the Pascal’s triangle respectively. For simplicity, we first consider the
case where $p>3$ when $(-1)^{r-1}=+1$. We will return to the case of $p=2$ and
$p=3$ later. After expanding around the $\pm\mathbf{K}=\pm(4\pi/3,0)$ points
we have that

$\displaystyle\hat{\Psi}_{0,0}$
$\displaystyle=\hat{\Psi}(\mathbf{r}_{2}=0)\sim 2\psi_{1}(\mathbf{r}_{2});$
$\displaystyle\hat{\Psi}_{0,r}$
$\displaystyle=\hat{\Psi}(\mathbf{r}_{1}=p^{k}\mathbf{a}_{2})$
$\displaystyle\sim 2\cos(\frac{2\pi
p^{k}}{3})\psi_{1}(\mathbf{r}_{1})-2\sin(\frac{2\pi
p^{k}}{3})\psi_{2}(\mathbf{r}_{1});$ $\displaystyle\hat{\Psi}_{r,0}$
$\displaystyle=\hat{\Psi}(\mathbf{r}_{3}=p^{k}\mathbf{a}_{1})$
$\displaystyle\sim 2\cos(\frac{4\pi
p^{k}}{3})\psi_{1}(\mathbf{r}_{3})-2\sin(\frac{4\pi
p^{k}}{3})\psi_{2}(\mathbf{r}_{3}).$

For $p\neq 3$, $\cos(2\pi p^{k}/3)=\cos(4\pi p^{k}/3)=-1/2$ and $\sin(4\pi
p^{k}/3)=-\sin(2\pi p^{k}/3)=\mp\sqrt{3}/2$. So we have after expanding
$\hat{\Psi}_{0,0}+\hat{\Psi}_{r,0}+\hat{\Psi}_{0,r}\sim[2\psi_{1}(\mathbf{r}_{2})-\psi_{1}(\mathbf{r}_{1})-\psi_{1}(\mathbf{r}_{3})]\mp\sqrt{3}[\psi_{2}(\mathbf{r}_{1})-\psi_{2}(\mathbf{r}_{3})]$.
Since $\psi_{1}$ and $\psi_{2}$ are uncoupled in the Lagrangian, the
expectation of $U_{p}(r)$ factorizes into the product of two simpler
correlators

$\displaystyle\langle
e^{\frac{2\pi\mathrm{i}}{p}[(2\psi_{1}(x_{2})-\psi_{1}(x_{1})-\psi_{1}(x_{3}))\mp\sqrt{3}(\psi_{2}(x_{1})-\psi_{2}(x_{3}))]}\rangle$
(18) $\displaystyle\sim$ $\displaystyle
e^{\frac{4\pi^{2}}{p^{2}}[2\langle\psi_{1}(x_{2})\psi_{1}(x_{1})\rangle+2\langle\psi_{1}(x_{2})\psi_{1}(x_{3})\rangle-\langle\psi_{1}(x_{1})\psi_{1}(x_{3})\rangle]}$
(20) $\displaystyle\times$ $\displaystyle
e^{\frac{12\pi^{2}}{p^{2}}\langle\psi_{2}(x_{1})\psi_{2}(x_{3})\rangle}.$ (22)

Here the Euclidean propagator
$\Delta(x-x^{\prime})=\langle\psi_{i}(x)\psi_{i}(x^{\prime})\rangle$ is
computed according to a single copy of the Lifshitz theory. It is known that
the propagator is logarithmic,

$\displaystyle\Delta(x-x^{\prime})=\int^{\Lambda}\frac{\differential[3]{q}}{(2\pi)^{3}}\text{
}\frac{e^{\mathrm{i}q\cdot(x-x^{\prime})}}{\omega^{2}+\rho_{4}q^{4}}=-\frac{\ln(\Lambda\absolutevalue{x-x^{\prime}})}{4\pi\sqrt{\rho_{4}}}$
(23)

$\Lambda$ is the momentum cutoff. Since the propagator only depends on the
distance $\absolutevalue{x-x^{\prime}}$, and our fields lie at the corners of
an equilateral of length $r$, $\langle U_{p}(r)\rangle$ becomes

$\displaystyle\langle
U_{p}(r)\rangle\sim\exp\left(\frac{24\pi^{2}}{p^{2}}\Delta(r)\right)\sim
r^{-\frac{6\pi}{p^{2}\sqrt{\rho_{4}}}}.$ (24)

Hence, the patch symmetry operator decays as a power-law with exponent
$\delta_{p}=\frac{6\pi}{p^{2}\sqrt{\rho_{4}}}$.

For $p=2$, the calculation is exactly the same but we must be careful since
$(-1)^{r-1}=-1$ now. However, since the original height fields $\hat{\Psi}$
are quantized, as the operator
$U_{p}(r)=e^{\mathrm{i}\pi(\hat{\Psi}_{0,r}-(\hat{\Psi}_{0,0}+\hat{\Psi}_{r,0}))}=e^{\mathrm{i}\pi(\hat{\Psi}_{0,r}+\hat{\Psi}_{0,0}+\hat{\Psi}_{r,0})}$,
the latter operator having a sensible algebriac decay computable in the
Lifshitz theory with exponent $\delta_{2}=\frac{3\pi}{2\sqrt{\rho_{4}}}$.

Now for $p=3$, $\hat{\Psi}_{0,0},\hat{\Psi}_{r,0},\hat{\Psi}_{0,r}$ all lie on
the same sublattice. As a consequence, after expansion we have

$\displaystyle\hat{\Psi}_{0,0}$
$\displaystyle=\hat{\Psi}(\mathbf{r}_{2}=0)\sim 2\psi_{1}(\mathbf{r}_{2});$
$\displaystyle\hat{\Psi}_{0,r}$
$\displaystyle=\hat{\Psi}(\mathbf{r}_{1}=3^{k}\mathbf{a}_{2})\sim
2\psi_{1}(\mathbf{r}_{1});$ $\displaystyle\hat{\Psi}_{r,0}$
$\displaystyle=\hat{\Psi}(\mathbf{r}_{3}=3^{k}\mathbf{a}_{1})\sim
2\psi_{1}(\mathbf{r}_{3}).$

Again making use of the fact that the initial fields were quantized we can
compute the correlation of
$U_{3}(r)=e^{\frac{2\pi\mathrm{i}}{3}(\hat{\Psi}_{0,r}+\hat{\Psi}_{0,0}+\hat{\Psi}_{r,0})}=e^{\frac{2\pi\mathrm{i}}{3}(\hat{\Psi}_{0,r}+\hat{\Psi}_{0,0}-2\hat{\Psi}_{r,0})}$.
The latter form of the operator is an observable which has sensible decay once
computed in the Lifshitz theory. This is because the sum of coefficients in
front of the fields vanishes in this form, a necessary requirement to find a
non-trivial result for the correlation function. After expansion, we must
compute in the Lifshitz theories

$\displaystyle\langle
e^{\frac{4\pi\mathrm{i}}{3}[\psi_{1}(x_{1})+\psi_{1}(x_{2})-2\psi_{1}(x_{3})]}\rangle$
(25) $\displaystyle\sim$ $\displaystyle
e^{\frac{16\pi^{2}}{9}[2\langle\psi_{1}(x_{1})\psi_{1}(x_{3})\rangle+2\langle\psi_{1}(x_{2})\psi_{1}(x_{3})\rangle-\langle\psi_{1}(x_{2})\psi_{1}(x_{1})\rangle]},$
(27) $\displaystyle=$ $\displaystyle e^{\frac{48\pi^{2}}{9}\Delta(r)}\sim
r^{-\frac{4\pi}{3\sqrt{\rho_{4}}}}.$ (29)

This is a slightly stronger decay than we found for $p\neq 3$, though the
decay remains power-law. Hence, we conclude that the decay of the patch
symmetry operator at the multi-critical point is $\langle U_{p}(r)\rangle\sim
r^{-\frac{6\pi}{p^{2}\sqrt{\rho_{4}}}}$ with this exponent being twice as
large for $p=3$.

We can compare the behavior of the patch operator at the multi-critical point
versus in the disordered phase, which as we argued before is the generic phase
of the model Eq. 1. Inside the disordered phase, the vertex operator $\alpha$
is relevant, and $\psi$ is pinned at the minima of the vertex operators, with
suppressed fluctuation. In this case, $\Delta(x-x^{\prime})$ decays
exponentially with distance, hence $\langle U_{p}(r)\rangle$ is expected to
saturate to a finite value with long distance. This is consistent with the
intuition one gather from other models: in the disordered phase, the disorder
operator usually has enhanced correlation compared with the critical point.

###  II.4 Multi-criticality of descendent $Z_{p}$ models

We have shown that upon fine-tuning the Pascal’s triangle model, it can reach
a multi-critical point. We can consider if there exists a similar multi-
critical point for the descendent fractal models with one of the fractal
symmetries. For example, a descendent fractal model can be constructed by
adding an operator to the theory that energetically favors the eigenvalues of
the boson phase to be $\frac{2\pi}{p}Z_{p}$ valued. The appropriately
perturbed Pascal’s triangle model is

$\displaystyle H_{p}=H-g\sum_{i}\cos(p\hat{\theta}_{i}).$ (30)

Where $H$ is given by Eq. 1. If the $Z_{p}$ anisotropy is irrelevant at the
multi-critical point discussed above, then we can conclude that all $Z_{p}$
fractal models attained by (weakly) breaking the rotor degrees of freedom down
share the same multi-critical point as the Pascal’s triangle model upon some
fine-tuning.

To analyze the relevance of the $Z_{p}$ anisotropy, we need to evaluate the
perturbation of the $g$ operator in Eq. 30 at the multi-critical point
discussed above. If an operator is relevant, it would make infrared divergent
contribution to the partition function at certain order of perturbation. For
example, in the XY-plaquette model, If one turns on the same $Z_{p}$
anisotropy $\sum_{i}g\cos(p\hat{\theta}_{i})$, the fourth order perturbation
of $g$ will have nonzero expectation value at the fixed point with $g=0$, as
long as the four operators are located at the four corners of a rectangle
Paramekanti et al. (2002), as required by the type-I subsystem symmetry.
Though this expectation value decays with the size of the rectangle, it was
shown that it can lead to infrared divergent contribution to the partition
function, especially when the $Z_{p}$ operators are paired up as dipoles.
Hence the $Z_{p}$ anisotropy can generate relevant dipole operators. But in
our case, one can verify that, in order to preserve the fractal symmetries at
the multi-critical point, at any finite order of perturbation of $g$, the
expectation value
$\langle\cos(p\hat{\theta}_{i_{1}})\cdots\cos(p\hat{\theta}_{i_{n}})\rangle$
always vanishes at large distance due to violation of the Pascal’s triangle
symmetry. Another way to look at the difference between our case and the XY-
plaquette model is that, the nonzero four point correlation function of the
XY-plaquette model can be viewed as the product of the
$\cos(\nabla_{x}\nabla_{y}\hat{\theta})$ term in the Hamiltonian over a
rectangular area; but in our case a product of
$\cos(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})$ over certain area
$A$ does not reduce to a simple finite composite operator when $A$ is large.

We stress here that, our argument for $g\cos(p\hat{\theta})$ being irrelevant
is based on perturbative expansion of $g$ at the fine-tuned fixed point where
$\hat{\theta}$ takes continuous value. This argument is only supposed to hold
for weak $Z_{p}$ anisotropy.

###  II.5 Relation of the Pascal’s triangle model to a gauge theory

Apart from the duality between the Pascal’s triangle model and the height
model, there exists another separate representation of the Pascal’s triangle
model in terms of a gauge theory. As we discussed in previous subsections, the
dual height model of the Pascal’s triangle model has two components of scalar
fields at low energy, then following the standard “boson-vortex” duality
procedure Peskin (1978); Dasgupta and Halperin (1981); Fisher and Lee (1989),
the dual of the height model should be two ${\rm U}(1)$ gauge fields. In this
section we explore the relation between the Pascal’s triangle model and a
${\rm U}(1)\times{\rm U}(1)$ gauge theory.

We define the canonical conjugate fields $\hat{E}_{i}$ and $\hat{A}_{i}$ on
the sites of a honeycomb lattice as

$\displaystyle\hat{E}_{i\in\mathcal{A}}=-(\hat{\Psi}_{1}+\hat{\Psi}_{2}+\hat{\Psi}_{3})_{\bar{\triangle}},$
(31) (32)
$\displaystyle\hat{E}_{i\in\mathcal{B}}=-(\hat{\Psi}_{1}+\hat{\Psi}_{2}+\hat{\Psi}_{3})_{\bar{\bigtriangledown}}$
(33) (34)
$\displaystyle\hat{\phi}_{\bar{j}}=(\hat{A}_{1}+\hat{A}_{2}+\hat{A}_{3}+\hat{A}_{4}+\hat{A}_{5}+\hat{A}_{6})_{\hexagon}$
(35)

where $\mathcal{A}$ and $\mathcal{B}$ are the two sublattices of the honeycomb
lattice. This honeycomb lattice has the original triangular lattice as one of
its sublattices, say $\mathcal{A}$. The electric fields $\hat{E}_{i}$ are
forced to obey the Gauss’ law constraint

$\displaystyle G_{\hexagon}=\sum_{i\in\hexagon}\eta_{i}\hat{E}_{i}\equiv
0\text{ with }\eta_{i}=\begin{cases}+1&i\in\mathcal{A}\\\
-1&i\in\mathcal{B}\end{cases}.$ (36)

The Hamiltonian of the resultant model, which we will refer to as the “six-
boson model”, is

$\displaystyle
H_{\text{sb}}=-K\sum_{\hexagon}\cos(\sum_{i\in\hexagon}\hat{A}_{i})+\frac{U}{2}\sum_{i\in\mathcal{A}}\hat{E}_{i}^{2}+\sum_{i,i^{\prime}\in\mathcal{A}}\frac{J}{2}\hat{E}_{i}\hat{E}_{i^{\prime}}.$
(37)

Figure 2: This figure shows an example of a gauge transformation that fixes a
single given $\hat{A}_{i\in\mathcal{B}}\equiv 0$, and leaves the rest of the
gauge fields on the $\mathcal{B}$ sublattice intact. Successive gauge
transformations of this form will fix all $\hat{A}_{i\in\mathcal{B}}=0$. Such
an accumulative gauge transformation that achieves this is given by Eq. 39 for
every $\alpha_{\hexagon}$ where $\tilde{\mathcal{R}}(\infty)$ is centered at
the top vertex of $\hexagon$.

In addition, this six-boson model is directly related to the Pascal’s triangle
model with or without fine-tuning by gauge fixing. The gauge symmetry
generated by the Gauss operator $G_{\hexagon}$ transforms the gauge fields
according to

$\displaystyle\hat{A}_{i}\rightarrow\hat{A}_{i}+\eta_{i}\sum_{\hexagon\ni
i}\alpha_{\hexagon}.$ (38)

There is actually a valid gauge where the gauge fields on one of the two
sublattices is completely gauged away. A gauge transformation that sets any
single $\hat{A}_{i\in\mathcal{B}}$ to zero and does not change any other
$\hat{A}_{i\in\mathcal{B}}$ is shown in Fig. 2. Using this gauge
transformation, an accumulative gauge transformation can be constructed that
fixes the gauge field to be zero at every site in the $\mathcal{B}$
sublattice. This is achieved by making the following choice for every
$\alpha_{\hexagon}$ appropriately:

$\displaystyle\alpha_{\hexagon}=\sum_{i\in\tilde{\mathcal{R}}(\infty)}(-1)^{i_{2}}{i_{2}\choose-
i_{1}}\hat{A}_{i\in\mathcal{B}},$ (39)

Specifically for $\alpha_{\hexagon}$ where $\hexagon$ is the honeycomb whose
top corner is at the origin, the appropriate region being summed over is
$\tilde{\mathcal{R}}(\infty)=\\{(i_{1},i_{2})\in\mathcal{B}:i_{2}\geq\infty,0\leq-
i_{1}\leq i_{2}\\}$.

Despite the form of the gauge transformation in Eq. 38, the gauge group of the
six-boson model is actually ${\rm U}(1)\times{\rm U}(1)$. This should be
expected from duality as each of the two scalar fields of the height model is
dual to a ${\rm U}(1)$ gauge field in $(2+1)d$. As a consequence, there are
two sets of Gauss’ law constraints that generate two separate gauge
transformations: one set where each $\hexagon$ lies on two of the sublattices
of the dual triangular lattice and a second set where each $\hexagon$ lies on
a different choice of two sublattices. The gauge group can also be directly
verified in the continuum limit of the Gaussian approximation of Eq. 37. After
the gauge fields are expanded around the $\pm\mathbf{K}$ points in the BZ, the
six-boson model reduces to two copies of QED3 up to a field redefinition.

Other correlation functions of the Pascal’s triangle model at the multi-
critical point may be computed using the six-boson model. In particular,
Wilson loops of the six-boson model end up reducing to correlation functions
of products of the three-body operator
$e^{i(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})}$ in the
$\hat{A}_{i\in\mathcal{B}}\equiv 0$ gauge. These Wilson loops are given by
staggered sums of the six-body interactions $\sum_{i\in\hexagon}\hat{A}_{i}$
in some region so that the support in the bulk cancels and what remains is a
loop of staggered gauge fields on the boundary. More precisely,

$\displaystyle\langle W[\eta]\rangle_{\text{sb}}$ $\displaystyle=$
$\displaystyle\bigg{\langle}\exp(i\sum_{\hexagon}\eta_{\hexagon}\bigg{(}\sum_{i\in\hexagon}\hat{A}_{i}\bigg{)})\bigg{\rangle}_{\text{sb}}$
(40) $\displaystyle\overset{\hat{A}_{i\in\mathcal{B}}\equiv
0}{\longrightarrow}$
$\displaystyle\bigg{\langle}\prod_{\bigtriangledown}\exp(i\eta_{\bigtriangledown}(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3}))\bigg{\rangle}_{\text{Pascal}}$
(42)

after a relabeling of $\hat{A}_{i\in\mathcal{A}}$ as $\hat{\theta}_{i}$.

In addition, the $z=2$ Lifshitz theory is known to have a duality to a gauge
theory. In particular, the field theory of the fine-tuned height model at the
multi-critical point is dual to

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\frac{1}{2\rho_{4}}\bigg{(}A_{x}^{(1)}\frac{\partial_{\tau}^{2}}{\partial_{y}^{2}}A_{x}^{(1)}+A_{y}^{(1)}\frac{\partial_{\tau}^{2}}{\partial_{x}^{2}}A_{y}^{(1)}\bigg{)}+\frac{1}{2}(\nabla\times\vec{A}^{(1)})^{2}$
(43) $\displaystyle+$
$\displaystyle(\vec{A}^{(1)}\leftrightarrow\vec{A}^{(2)})$ (45)

for two ${\rm U}(1)$ gauge fields $\vec{A}^{(1)}$ and $\vec{A}^{(2)}$. This
implies that the low energy theory of a fine-tuned version of the six-boson
model, which reduces to the appropriate fine-tuned Pascal’s triangle model, is
equal to Eq. 45 up to a redefinition of the gauge fields.

Figure 3:  On the triangular lattice, a staggered product of
$e^{i(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})}$ and
$e^{-i(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})}$ on the
$\tilde{A}$ and $\tilde{B}$ sublattices of the dual triangular lattice inside
patch $\mathcal{P}$ has support only on $\partial\mathcal{P}$. The region
inside the dashed gray line is $\mathcal{P}$, with blue dots representing
sites in $\tilde{A}\cap\mathcal{P}$ and red dots representing sites in
$\tilde{B}\cap\mathcal{P}$. This product reduces to some staggered
$\hat{\theta}$ variables leftover on $\partial\mathcal{P}$ as shown by the
blue and red crosses which is
$e^{i\int_{\partial\mathcal{P}}\gamma_{i}\hat{\theta}_{i}}$ in Eq. 46.

It is known that the Wilson loops decay as a perimeter law in this theory, Eq.
45. As this is the case, correlation functions of products of the three-body
term in the Pascal’s triangle model will also decay as a perimeter law. On the
triangular lattice, we can consider a local patch let us denote as
$\mathcal{P}$. This patch contains a set of downward plaquettes which in
themselves form the sites of the dual triangular lattice. If we now take a
product of the three-body operators
$e^{i(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})}$ on two of the
three dual sublattices of the plaquettes inside $\mathcal{P}$ staggered on one
of the sublattices, we recover an operator that is a product of $e^{\pm
i\hat{\theta}}$ with support only on $\partial\mathcal{P}$. That is,

$\displaystyle\prod_{\bigtriangledown\in\tilde{A}\cap\mathcal{P}}e^{i(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})}\prod_{\bigtriangledown\in\tilde{B}\cap\mathcal{P}}e^{-i(\hat{\theta}_{1}+\hat{\theta}_{2}+\hat{\theta}_{3})}=e^{i\int_{\partial\mathcal{P}}\gamma_{i}\hat{\theta}_{i}}$
(46)

where $\gamma_{i}=+1,-1,0$ depending on whether $\hat{\theta}_{i}$ on the
boundary lies on a downward triangle on the $\tilde{A},\tilde{B},\tilde{C}$
sublattices of the dual triangular lattice. However, this operator is really
just a Wilson loop of the six-boson model in the gauge
$\hat{A}_{i\in\mathcal{B}}\equiv 0$ of the form Eq. 42 and hence decays as a
perimeter law with the number of boson phases on $\partial\mathcal{P}$. This
is explicitly shown in Fig. 3. But in the quantum disordered phase, which is
the generic phase of either the Pascal’s triangle model or the six-boson
model, the Wilson loop will decay with an area law. In terms of the gauge
theory language the quantum disordered phase is the confined phase. The area
law can be shown by starting with the ground state with large $U$, and perform
perturbation of $K$ operator in the Hamiltonian.

##  III Pascal’s tetrahedron model

###  III.1 The Hamiltonian and its symmetries

There is a natural generalization of the Pascal’s triangle model to $(3+1)d$
which we refer to as the Pascal’s tetrahedron model. On the cubic lattice, the
Hamiltonian of the Pascal’s tetrahedron model reads

$\displaystyle
H=-K\sum_{i}\cos(\sum_{\mu=0,\hat{x},\hat{y},\hat{z}}\hat{\theta}_{i+\mu})+\frac{U}{2}\sum_{i}\hat{n}_{i}^{2}.$
(47)

The four-body interaction
$\cos(\hat{\theta}_{i}+\hat{\theta}_{i+\hat{x}}+\hat{\theta}_{i+\hat{y}}+\hat{\theta}_{i+\hat{z}})$
involves the sum of boson phases over a tetrahedron on the cubic lattice; it
is the natural generalization of the three-body interaction of Eq. 1. The
“Pascal’s tetrahedron symmetries” are also made up of an infinite series of
$Z_{p}$ fractal symmetries. The action of the fractal symmetries on a finite
tetrahedron with top corner at the origin of length $L=p^{k}-1$ is generated
by the following patch symmetry operator

$\displaystyle U_{p}$ $\displaystyle=$ $\displaystyle
e^{2\pi\mathrm{i}Q_{p}/p}$ (48) $\displaystyle Q_{p}$ $\displaystyle=$
$\displaystyle\sum_{i\in\mathcal{R}(p)}(-1)^{-i}\frac{(-i_{x}-i_{y}-i_{z})!}{(-i_{x})!(-i_{y})!(-i_{z})!}\hat{n}_{i}.$
(50)

The action of this operator rotates the boson phase according to
$\hat{\theta}_{i}\rightarrow\hat{\theta}_{i}+(2\pi/p)(-1)^{i}(-i_{x}-i_{y}-i_{z})!/(-i_{x})!(-i_{y})!(-i_{z})!$
on the tetrahedral region
$\mathcal{R}(p)\equiv\\{(i_{x},i_{y},i_{z}):i_{x},i_{y},i_{z}\leq
0,i_{x}+i_{y}+i_{z}\in[-L,0]\\}$. Successive transformations of this form over
extensively large regions will give rise to exact fractal symmetries, e.g. an
exact fractal symmetry with $p=2$ can be generated by three of these
transformations for lattice size $L=2^{k}-1$.

Despite Eq. 47 being the straightforward generalization of the Pascal’s
triangle model to $(3+1)d$, the Pascal’s tetrahedron model manifests another
type-I planar subsystem symmetry of which there is no analogue in $(2+1)d$. We
may decompose the cubic lattice into four sublattices each spanned by vectors
$\mathbf{b}_{1}=\hat{x}+\hat{y}-\hat{z},\mathbf{b}_{2}=\hat{x}+\hat{z}-\hat{y},\mathbf{b}_{3}=\hat{y}+\hat{z}-\hat{x}$.
The first sublattice is given by
$\text{span}\\{\mathbf{b}_{1},\mathbf{b}_{2},\mathbf{b}_{3}\\}$ and the other
three are obtained from subsequent shifts by $\hat{x},\hat{y}$ and $\hat{z}$
respectively. Now, consider a plane $\mathfrak{L}$ on the cubic lattice with
normal $\hat{n}_{1}=\frac{1}{\sqrt{2}}(\hat{y}+\hat{z})$. We construct this
plane by $\mathfrak{L}=\mathfrak{L}_{1}\cup\mathfrak{L}_{2}$ where
$\mathfrak{L}_{1}=\text{span}\\{\mathbf{b}_{1}+\mathbf{b}_{2},\mathbf{b}_{2}\\}$
and $\mathfrak{L}_{2}$ is attained from $\mathfrak{L}_{1}$ by a shift of
$\hat{x}$. The Pascal’s tetrahedron model is invariant under the following
transformation on $\mathfrak{L}$:

$\displaystyle{\rm
U}(1)^{\text{(sub)}}:\hat{\theta}_{i\in\mathfrak{L}_{1}}\rightarrow\hat{\theta}_{i\in\mathfrak{L}_{1}}+\alpha,\
\
\theta_{i\in\mathfrak{L}_{2}}\rightarrow\hat{\theta}_{i\in\mathfrak{L}_{2}}-\alpha$
(51)

Translations of the support of this transformation in the $\hat{y}$ and
$\hat{z}$ directions will be exact symmetries as well. Furthermore, there
exist two other sets of mutually orthogonal planes with normals
$\hat{n}_{2}=\frac{1}{\sqrt{2}}(\hat{x}+\hat{y}),\hat{n}_{3}=\frac{1}{\sqrt{2}}(\hat{x}+\hat{z})$
of which a staggered rotation such as Eq. 51 leaves the Pascal’s tetrahedron
model invariant. These collection of planar transformations form the entire
planar ${\rm U}(1)^{\text{(sub)}}$ symmetry of Eq. 47.

In the classical limit of the Pascal’s tetrahedron model ($U=0$), there exist
both domain walls and fractonlike excitations. Firstly there exist fractonlike
excitations created at the tetrahedra located at the corners of a Pascal’s
tetrahedron mod $p$ by acting with the fractal transformation in a local
patch. Like the $(2+1)d$ case, these excitations cost $K[1-\cos(2\pi/p)]$
energy above the $\theta_{i}=0$ configuration. In addition, there are one-
dimensional loop excitations analogous to domain walls of the $(2+1)d$ XY
model created by acting with the ${\rm U}(1)^{(\text{sub})}$ symmetry on a
patch of a plane such as $\mathfrak{L}$. Starting with the uniform
configuration $\theta_{i}\equiv 0$, these loop excitations cost an energy of
$\sim N_{\text{loop}}K[1-\cos\alpha]$ above such a ground-state configuration.

There is a dual of the Pascal’s tetrahedron model similar to that of Eq. 6. If
we define variables on the dual cubic lattice whose lattice sites are located
at the center of the cubes of the original lattice

$\displaystyle\hat{\theta}_{i}+\hat{\theta}_{i+\hat{x}}+\hat{\theta}_{i+\hat{y}}+\hat{\theta}_{i+\hat{z}}$
$\displaystyle=$ $\displaystyle\hat{\phi}_{\bar{i}},$ (52)
$\displaystyle\hat{\Psi}_{\bar{i}}+\hat{\Psi}_{\bar{i}-\hat{x}}+\hat{\Psi}_{\bar{i}-\hat{y}}+\hat{\Psi}_{\bar{i}-\hat{z}}$
$\displaystyle=$ $\displaystyle-\hat{n}_{i}.$ (54)

Such a mapping preserves the bosonic algebra, and yields a dual height model
of the Pascal’s tetrahedron model

$\displaystyle
H_{h}=-K\sum_{\bar{i}}\cos(\hat{\phi}_{\bar{i}})+\frac{U}{2}\sum_{\bar{i}}\bigg{(}\sum_{\mu=0,\hat{x},\hat{y},\hat{z}}\hat{\Psi}_{\bar{i}-\mu}\bigg{)}^{2}.$
(55)

This duality involving the $Z_{2}$, or Sierpinski-tetrahedron, version of this
model was pointed out in Ref. Garrahan, 2002; Devakul et al., 2019.

This dual height model does not explicitly possess the Pascal’s tetrahedron
symmetries of the original model as the $\hat{\Psi}$ fields are discrete,
though it does retain a dual version of the ${\rm U}(1)^{\text{(sub)}}$
symmetry. This dual subsystem symmetry is realized in the spectrum via nodal
lines of zero energy modes. It is straightforward to calculate the spectrum of
the height model as
$\epsilon(\mathbf{q})=U[2+\cos(q_{x})+\cos(q_{y})+\cos(q_{z})+\cos(q_{x}-q_{y})+\cos(q_{x}-q_{z})+\cos(q_{y}-q_{z})]$.
Despite the modest form of the Hamiltonian and the spectrum, we know that the
height model has a dual ${\rm U}(1)^{(\text{sub})}$ symmetry and so its
spectrum must display nodal lines where $\epsilon(\mathbf{q})=0$. In fact
there are six nodal line segments in the BZ, e.g.
$\\{\mathbf{q}=(t-\pi,t,\pi):t\in[0,\pi]\\}$. If expanded along this line, the
dispersion of the height field becomes

$\displaystyle\epsilon(\mathbf{q})$ $\displaystyle\sim$
$\displaystyle-2Uq_{x}q_{y}+2U\cos(t)(q_{x}q_{z}-q_{y}q_{z})$ (56)
$\displaystyle+$ $\displaystyle U(q_{x}^{2}+q_{y}^{2}+q_{z}^{2}).$ (58)

As a consequence, an analysis of a potential multicritical point of the
Pascal’s tetrahedron model must involve a similar height field expansion as in
the analysis in $(2+1)d$, but around each of the nodal lines. We leave the
detailed discussion of the RK point of tetrahedron model to future studies.

Figure 4: By making a staggered product of
$\exp(\mathrm{i}(\hat{A}^{(1)}_{i}+\dots+\hat{A}^{(2)}_{i+\hat{y}+\hat{z}}))$
term in Eq. 73 on the cubes labelled as blue and red which belong to different
dual sublattices, both gauge fields $\hat{A}^{(1)}$ and $\hat{A}^{(2)}$ cancel
out in the bulk; the remaining fields are located on the boundary. The
remaining $\hat{A}^{(1)}$ and $\hat{A}^{(2)}$ fields are labelled as crosses
in the top and bottom figures respectively.

###  III.2 Relation to the U(1) Haah’s code

Like the Pascal’s triangle model, the Pascal’s tetrahedron model is related to
a gauge theory by gauge fixing. However, while in the $(2+1)d$ this theory was
a conventional ${\rm U}(1)\times{\rm U}(1)$ gauge theory, we will show that
the appropriate gauge theory in $(3+1)d$ is the pure (i.e. uncoupled from
matter, the Gauss’ operator is non-dynamical) version of the ${\rm U}(1)$
Haah’s code. The dual height model in Eq. 55 has a second duality. We can
define two new electric and gauge fields according to

$\displaystyle\hat{E}_{i}^{(1)}$ $\displaystyle=$
$\displaystyle-(\hat{\Psi}_{\bar{i}}+\hat{\Psi}_{\bar{i}-\hat{x}}+\hat{\Psi}_{\bar{i}-\hat{y}}+\hat{\Psi}_{\bar{i}-\hat{z}}),$
(59) $\displaystyle\hat{E}_{i}^{(2)}$ $\displaystyle=$
$\displaystyle-(\hat{\Psi}_{\bar{i}}+\hat{\Psi}_{\bar{i}-\hat{x}-\hat{y}}+\hat{\Psi}_{\bar{i}-\hat{x}-\hat{z}}+\hat{\Psi}_{\bar{i}-\hat{y}-\hat{z}}),$
(61) $\displaystyle\hat{\phi}_{\bar{i}}$ $\displaystyle=$
$\displaystyle\hat{A}^{(1)}_{i}+\hat{A}^{(1)}_{i+\hat{x}}+\hat{A}^{(1)}_{i+\hat{y}}+\hat{A}^{(1)}_{i+\hat{z}}$
(63) $\displaystyle+$
$\displaystyle\hat{A}^{(2)}_{i}+\hat{A}^{(2)}_{i+\hat{x}+\hat{y}}+\hat{A}^{(2)}_{i+\hat{x}+\hat{z}}+\hat{A}^{(2)}_{i+\hat{y}+\hat{z}}.$
(65)

The $Z_{2}$ version of this mapping is equivalent to the “F-S” duality
outlined in Ref. Vijay et al., 2016. There exists a linear combination of
interactions of the form
$(\hat{\Psi}_{\bar{i}}+\hat{\Psi}_{\bar{i}-\hat{x}}+\hat{\Psi}_{\bar{i}-\hat{y}}+\hat{\Psi}_{\bar{i}-\hat{z}})$
and
$(\hat{\Psi}_{\bar{i}}+\hat{\Psi}_{\bar{i}-\hat{x}-\hat{y}}+\hat{\Psi}_{\bar{i}-\hat{x}-\hat{z}}+\hat{\Psi}_{\bar{i}-\hat{y}-\hat{z}})$
that vanishes, and as a result the electric fields of the dual must obey the
corresponding Gauss’ law constraint this time being

$\displaystyle G_{i}$ $\displaystyle=$
$\displaystyle\hat{E}^{(1)}_{i+\hat{x}+\hat{y}+\hat{z}}+\hat{E}^{(1)}_{i+\hat{x}}+\hat{E}^{(1)}_{i+\hat{y}}+\hat{E}^{(1)}_{i+\hat{z}}$
(66) $\displaystyle-$
$\displaystyle\hat{E}^{(2)}_{i+\hat{x}+\hat{y}+\hat{z}}-\hat{E}^{(2)}_{i+\hat{x}+\hat{y}}+\hat{E}^{(2)}_{i+\hat{x}+\hat{z}}-\hat{E}^{(2)}_{i+\hat{y}+\hat{z}}\equiv
0.$ (68)

The Hamiltonian of the resulting gauge theory takes the form of the rotor
limit of pure Haah’s code

$\displaystyle H$ $\displaystyle=$
$\displaystyle-K\sum_{i}\text{cos}(\hat{A}^{(1)}_{i}+\hat{A}^{(1)}_{i+\hat{x}}+\hat{A}^{(1)}_{i+\hat{y}}+\hat{A}^{(1)}_{i+\hat{z}}$
(69) $\displaystyle+$
$\displaystyle\hat{A}^{(2)}_{i}+\hat{A}^{(2)}_{i+\hat{x}+\hat{y}}+\hat{A}^{(2)}_{i+\hat{x}+\hat{z}}+\hat{A}^{(2)}_{i+\hat{y}+\hat{z}})$
(71) $\displaystyle+$
$\displaystyle\frac{U}{2}\sum_{i}(\hat{E}_{i}^{(1)})^{2}.$ (73)

Reducing the degrees of freedom to Ising spins via
$(\sigma^{z})^{(a)}_{j}=e^{i\hat{A}^{(a)}_{j}}$ and
$(\sigma^{x})^{(a)}_{j}=e^{\pi i\hat{E}_{j}^{(a)}}$, leaves us with a
perturbed version of the pure $Z_{2}$ Haah’s code. Fractional gauge fluxes of
this model created at the corners of products of $e^{i\hat{E}_{i}^{(1)}}$ or
$e^{i\hat{E}^{(2)}_{i}}$ in the shape of Pascal’s tetrahedra mod $p$ will be
completely immobile fractons. As each of these fractons costs an energy
$\Delta E=K[1-\cos(2\pi/p)]$ in the classical limit ($U=0$).

The Gauss operator generates a gauge transformation acting on the gauge fields
of

$\displaystyle\hat{A}^{(1)}_{i}$ $\displaystyle\rightarrow$
$\displaystyle\hat{A}^{(1)}_{i}+\alpha_{i-\hat{x}}+\alpha_{i-\hat{y}}+\alpha_{i-\hat{z}}+\alpha_{i-\hat{x}-\hat{y}-\hat{z}},$
(74) $\displaystyle\hat{A}^{(2)}_{i}$ $\displaystyle\rightarrow$
$\displaystyle\hat{A}^{(2)}_{i}-\alpha_{i-\hat{x}-\hat{y}}-\alpha_{i-\hat{x}-\hat{z}}-\alpha_{i-\hat{y}-\hat{z}}-\alpha_{i-\hat{x}-\hat{y}-\hat{z}}.$
(76)

The pure ${\rm U}(1)$ Haah’s code can be reduced to the Pascal’s tetrahedron
model upon gauge fixing and tuning of the electric couplings. We may choose

$\displaystyle\alpha_{j-\hat{x}-\hat{y}-\hat{z}}=\sum_{i\in\tilde{\mathcal{R}}(\infty)}(-1)^{i}\frac{(i_{x}+i_{y}+i_{z})!}{i_{x}!i_{y}!i_{z}!}\hat{A}_{i}^{(2)},$
(77)

with tetrahedral region $\tilde{\mathcal{R}}(\infty)$ being that of
$\mathcal{R}(\infty)$ with the opposite orientation centered at site $j$. Such
a transformation sets $\hat{A}^{(2)}_{j}\equiv 0$ at every site as with this
choice,
$\alpha_{j}+\alpha_{j+\hat{x}}+\alpha_{j+\hat{y}}+\alpha_{j+\hat{z}}=\hat{A}^{(2)}_{j}$
yielding exactly the Pascal’s tetrahedron model. This should be expected given
the dualities we have just described; our pure ${\rm U}(1)$ Haah’s code is
dual to the height model Eq. 55 which in turn is dual to Eq. 47, and hence
must be dynamically equivalent to the Pascal’s tetrahedron model.

Since the height field in the dual representation of the tetrahedron model is
a scalar field, then one would expect that the dual of a scalar field is a
2-form gauge field in three dimensional space. And a 2-form gauge field should
have “Wilson membranes” rather than Wilson lines. Indeed, one can show that
gauge invariant operators constructed as staggered products of the
$K$-operator,
$\exp(\mathrm{i}(\hat{A}_{i}^{(1)}+\dots+\hat{A}_{i+\hat{y}+\hat{z}}^{(2)}))$
in Eq. 73, in a three dimensional bulk of the system reduces to surface
products of $\hat{A}^{(1)}$ and $\hat{A}^{(2)}$ (Fig. 4). Specifically, one
can consider a box shaped patch of the cubic lattice. As each cube inside the
patch belongs to one of the four sublattices of the dual cubic lattice, a
product of $K$-operators staggered on two of the dual sublattices reduces to
an operator which only has remaining support on the boundary of the patch.
This operator appropriately reduces to a product of boson phases in the
boundary of the patch in the Pascal’s tetrahedron model upon gauge fixing.

##  IV Discussion

In this work we constructed a gapless multi-critical point from a lattice
model with an infinite series of fractal symmetries, which together are called
the Pascal’s triangle symmetry. This multi-critical point can be reached by
tuning a finite number of parameters of the generic Pascal’s triangle model,
and it is analogous to the RK point of the quantum dimer model in the sense
that it is gapless with a $z=2$ dispersion. At this multi-critical point,
physics can be analyzed conveniently in the dual description of the system;
and we demonstrate that the characteristic “patch operators” of all the
fractal symmetries have an expectation value that decay with a power law of
space, which is usually a sign of criticality. Potential generalizations to
three dimensional systems are also discussed, but we leave an in-depth
discussion of a potential RK point of the Pascal’s tetrahedron model to a
future work. As we mentioned earlier, gapless states with type-I subsystem
symmetries have been studied quite extensively since Ref. Paramekanti et al.,
2002, we hope our work can draw attention to the study of gapless systems with
type-II fractal symmetries.

We would like to acknowledge Sagar Vijay and Shang Liu for insightful
discussions. This work is supported by NSF Grant No. DMR-1920434, and the
Simons Investigator program. Part of the discussions in this work originated
from an unpublished work with Ashvin Vishwanath many years back.

## References

  * Nussinov and Ortiz (2009) Z. Nussinov and G. Ortiz, Annals of Physics 324, 977 (2009), ISSN 0003-4916, URL http://dx.doi.org/10.1016/j.aop.2008.11.002.
  * Aharony et al. (2013) O. Aharony, N. Seiberg, and Y. Tachikawa, Journal of High Energy Physics 2013 (2013), ISSN 1029-8479, URL http://dx.doi.org/10.1007/JHEP08(2013)115.
  * Gukov and Kapustin (2013) S. Gukov and A. Kapustin (2013), eprint 1307.4793.
  * Kapustin and Thorngren (2013a) A. Kapustin and R. Thorngren (2013a), eprint 1308.2926.
  * Kapustin and Thorngren (2013b) A. Kapustin and R. Thorngren (2013b), eprint 1309.4721.
  * Kapustin and Seiberg (2014) A. Kapustin and N. Seiberg, Journal of High Energy Physics 2014 (2014), ISSN 1029-8479, URL http://dx.doi.org/10.1007/JHEP04(2014)001.
  * Gaiotto et al. (2015) D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, Journal of High Energy Physics 2015 (2015), ISSN 1029-8479, URL http://dx.doi.org/10.1007/JHEP02(2015)172.
  * Hsin et al. (2019) P.-S. Hsin, H. T. Lam, and N. Seiberg, SciPost Physics 6 (2019), ISSN 2542-4653, URL http://dx.doi.org/10.21468/SciPostPhys.6.3.039.
  * Seiberg (2020) N. Seiberg, SciPost Physics 8 (2020), ISSN 2542-4653, URL http://dx.doi.org/10.21468/SciPostPhys.8.4.050.
  * McGreevy (2022) J. McGreevy, _Generalized symmetries in condensed matter_ (2022), URL https://arxiv.org/abs/2204.03045.
  * Gaiotto et al. (2017) D. Gaiotto, A. Kapustin, Z. Komargodski, and N. Seiberg, Journal of High Energy Physics 2017 (2017), URL https://doi.org/10.1007%2Fjhep05%282017%29091.
  * Nandkishore and Hermele (2019) R. M. Nandkishore and M. Hermele, Annual Review of Condensed Matter Physics 10, 295 (2019), eprint https://doi.org/10.1146/annurev-conmatphys-031218-013604, URL https://doi.org/10.1146/annurev-conmatphys-031218-013604.
  * Pretko et al. (2020) M. Pretko, X. Chen, and Y. You, International Journal of Modern Physics A 35, 2030003 (2020), eprint https://doi.org/10.1142/S0217751X20300033, URL https://doi.org/10.1142/S0217751X20300033.
  * Paramekanti et al. (2002) A. Paramekanti, L. Balents, and M. P. A. Fisher, Phys. Rev. B 66, 054526 (2002), URL https://link.aps.org/doi/10.1103/PhysRevB.66.054526.
  * Xu and Fisher (2007) C. Xu and M. P. A. Fisher, Phys. Rev. B 75, 104428 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.75.104428.
  * Ma and Pretko (2018) H. Ma and M. Pretko, Phys. Rev. B 98, 125105 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.125105.
  * You et al. (2020) Y. You, Z. Bi, and M. Pretko, Phys. Rev. Research 2, 013162 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.013162.
  * Yuan et al. (2020) J.-K. Yuan, S. A. Chen, and P. Ye, Phys. Rev. Research 2, 023267 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.023267.
  * Gromov et al. (2020) A. Gromov, A. Lucas, and R. M. Nandkishore, Phys. Rev. Research 2, 033124 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033124.
  * Seiberg and Shao (2021) N. Seiberg and S.-H. Shao, SciPost Physics 10 (2021), URL https://doi.org/10.21468%2Fscipostphys.10.2.027.
  * Rudelius et al. (2021) T. Rudelius, N. Seiberg, and S.-H. Shao, Phys. Rev. B 103, 195113 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.195113.
  * Gorantla et al. (2021) P. Gorantla, H. T. Lam, N. Seiberg, and S.-H. Shao, Phys. Rev. B 104, 235116 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.104.235116.
  * Burnell et al. (2021) F. J. Burnell, T. Devakul, P. Gorantla, H. T. Lam, and S.-H. Shao, _Anomaly inflow for subsystem symmetries_ (2021), URL https://arxiv.org/abs/2110.09529.
  * Haah (2011) J. Haah, Phys. Rev. A 83, 042330 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.83.042330.
  * Zhou et al. (2021) Z. Zhou, X.-F. Zhang, F. Pollmann, and Y. You (2021), eprint 2105.05851.
  * Vasiloiu et al. (2020) L. M. Vasiloiu, T. H. E. Oakes, F. Carollo, and J. P. Garrahan, Phys. Rev. E 101, 042115 (2020), URL https://link.aps.org/doi/10.1103/PhysRevE.101.042115.
  * Myerson-Jain et al. (2022) N. E. Myerson-Jain, S. Liu, W. Ji, C. Xu, and S. Vijay, Phys. Rev. Lett. 128, 115301 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.115301.
  * Newman and Moore (1999) M. E. J. Newman and C. Moore, Phys. Rev. E 60, 5068 (1999), URL https://link.aps.org/doi/10.1103/PhysRevE.60.5068.
  * Yoshida (2013) B. Yoshida, Phys. Rev. B 88, 125122 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.125122.
  * Rokhsar and Kivelson (1988) D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. 61, 2376 (1988), URL https://link.aps.org/doi/10.1103/PhysRevLett.61.2376.
  * FRADKIN and KIVELSON (1990) E. FRADKIN and S. KIVELSON, Modern Physics Letters B 04, 225 (1990), eprint https://doi.org/10.1142/S0217984990000295, URL https://doi.org/10.1142/S0217984990000295.
  * Fradkin et al. (2004) E. Fradkin, D. A. Huse, R. Moessner, V. Oganesyan, and S. L. Sondhi, Phys. Rev. B 69, 224415 (2004), URL https://link.aps.org/doi/10.1103/PhysRevB.69.224415.
  * Ardonne et al. (2004) E. Ardonne, P. Fendley, and E. Fradkin, Annals of Physics 310, 493 (2004), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/S0003491604000247.
  * Vishwanath et al. (2004) A. Vishwanath, L. Balents, and T. Senthil, Phys. Rev. B 69, 224416 (2004), URL https://link.aps.org/doi/10.1103/PhysRevB.69.224416.
  * Ji and Wen (2020) W. Ji and X.-G. Wen, Phys. Rev. Research 2, 033417 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.033417.
  * Kadanoff and Ceva (1971) L. P. Kadanoff and H. Ceva, Phys. Rev. B 3, 3918 (1971), URL https://link.aps.org/doi/10.1103/PhysRevB.3.3918.
  * Wu et al. (2021a) X.-C. Wu, W. Ji, and C. Xu, Journal of Statistical Mechanics: Theory and Experiment 2021, 073101 (2021a), URL https://doi.org/10.1088/1742-5468/ac08fe.
  * Zhao et al. (2021) J. Zhao, Z. Yan, M. Cheng, and Z. Y. Meng, Phys. Rev. Research 3, 033024 (2021), URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.033024.
  * Wu et al. (2021b) X.-C. Wu, C.-M. Jian, and C. Xu, SciPost Physics 11 (2021b), URL https://doi.org/10.21468%2Fscipostphys.11.2.033.
  * Wang et al. (2021a) Y.-C. Wang, M. Cheng, and Z. Y. Meng, Phys. Rev. B 104, L081109 (2021a), URL https://link.aps.org/doi/10.1103/PhysRevB.104.L081109.
  * Wang et al. (2021b) Y.-C. Wang, N. Ma, M. Cheng, and Z. Y. Meng, _Scaling of disorder operator at deconfined quantum criticality_ (2021b), URL https://arxiv.org/abs/2106.01380.
  * Peskin (1978) M. E. Peskin, Annals of Physics 113, 122 (1978), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/000349167890252X.
  * Dasgupta and Halperin (1981) C. Dasgupta and B. I. Halperin, Phys. Rev. Lett. 47, 1556 (1981), URL https://link.aps.org/doi/10.1103/PhysRevLett.47.1556.
  * Fisher and Lee (1989) M. P. A. Fisher and D. H. Lee, Phys. Rev. B 39, 2756 (1989), URL https://link.aps.org/doi/10.1103/PhysRevB.39.2756.
  * Garrahan (2002) J. P. Garrahan, Journal of Physics: Condensed Matter 14, 1571 (2002), URL https://doi.org/10.1088/0953-8984/14/7/314.
  * Devakul et al. (2019) T. Devakul, Y. You, F. J. Burnell, and S. L. Sondhi, SciPost Phys. 6, 7 (2019), URL https://scipost.org/10.21468/SciPostPhys.6.1.007.
  * Vijay et al. (2016) S. Vijay, J. Haah, and L. Fu, Phys. Rev. B 94, 235157 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.235157.

##  Appendix A Review of the XY-plaquette Model

When we discuss gapless states with subsystem symmetries, it is always useful
to first think about the XY-plaquette model Paramekanti et al. (2002). The XY-
plaquette model is a quantum rotor model with a type-I subsystem symmetry. The
Hamiltonian of the model on the square lattice reads

$\displaystyle
H_{\text{plaquette}}=-K\sum_{\square}\cos(\nabla_{x}\nabla_{y}\hat{\theta})+\frac{U}{2}\sum_{i}\hat{n}_{i}^{2}.$
(78)

The operator $\nabla_{x}\nabla_{y}\hat{\theta}_{i}$ is defined as
$\sum_{i\in\square}(-1)^{i}\hat{\theta}_{i}$ on each unit square plaquette.
This model admits a ${\rm U}(1)^{(\text{sub})}$ symmetry involving rotating
the boson phase along lines of constant $i_{x}$ or $i_{y}$:

$\displaystyle{\rm U}(1)^{\text{(sub)}}$ $\displaystyle:$
$\displaystyle\theta_{i_{x},i_{y}=\text{const}}\rightarrow\theta_{i_{x},i_{y}=\text{const}}+\alpha(i_{y}),$
(79) or $\displaystyle:$
$\displaystyle\theta_{i_{x}=\text{const},i_{y}}\rightarrow\theta_{i_{x}=\text{const},i_{y}}+\alpha(i_{x}).$
(81)

As per usual, the boson number and phase obey the canonical commutation
relation $[\hat{\theta}_{i},\hat{n}_{j}]=\mathrm{i}\delta_{ij}$. This model is
known to exhibit multiple dual representations. We will discuss two of these
here as they are analogous to the dualities involving the Pascal’s triangle
model discussed in the main text. Firstly, there is a dual to the height model
that can be attained by defining variables that live on the dual square
lattice

$\displaystyle\sum_{j\in\square}(-1)^{j}\hat{\theta}_{j}=\hat{\phi}_{\bar{j}},\
\ \ \ \hat{n}_{j}=-\sum_{\bar{j}\in\bar{\square}\text{ around}\
j}(-1)^{\bar{j}}\hat{\Psi}_{\bar{j}}.$ (82)

The $\bar{j}$ label the sites of the dual square lattice and live at the
center of the plaquettes of the original square lattice. The dual Hamiltonian
is a height version of the XY-plaquette model

$\displaystyle
H_{h}=-K\sum_{\bar{j}}\cos(\hat{\phi}_{\bar{j}})+\frac{U}{2}\sum_{\bar{\square}}\bigg{(}\nabla_{x}\nabla_{y}\Psi\bigg{)}^{2}$
(83)

This height dual exhibits a dual ${\rm U}(1)^{(\text{sub})}$ symmetry –
shifting the height fields $\hat{\Psi}_{\bar{j}}$ by an integer along a line
of constant $\bar{j}_{x}$ or $\bar{j}_{y}$ leaves the Hamiltonian invariant.

The XY-plaquette model has a gapless phase that allows us to take a Gaussian
expansion of the cosine term in Eq. 78. The (continuum limit) effective theory
of the Gaussian theory reads

$\displaystyle\mathcal{L}_{\text{plaquette}}=\frac{1}{2}(\partial_{\tau}\theta)^{2}+\frac{\rho}{2}(\partial_{x}\partial_{y}\theta)^{2}.$
(84)

We remind the readers that in the current paper the continuum limit means
infinite system size, constant lattice spacing. Due to the special subsystem
symmetry of the XY-plaquette model (and its various generalizations), one
cannot take lattice spacing $a$ to zero, which would lead to unphysical
conclusions. This effective theory is manifestly non-Lorentzian as it retains
part of the type-I subsystem symmetry: the Lagrangian is invariant under any
shift

$\displaystyle\theta(x,y,\tau)\rightarrow\theta(x,y,\tau)+f_{1}(x)+f_{2}(y),$
(85)

where $f_{1}$ and $f_{2}$ are completely arbitrary functions. This Gaussian
theory is gapless, as the dispersion of the modes is
$\omega(\mathbf{q})\sim\sqrt{\rho}\absolutevalue{q_{x}q_{y}}$. The subsystem
symmetry of the effective theory does not necessarily include the entire
subsystem symmetry of the original lattice model, as $f_{1}$ and $f_{2}$ are
presumably smooth functions of space. In fact, some of the low energy modes of
the original lattice near the boundary of the Brillouin zone are ignored in
the effective theory.

There is another dual of the height model Eq. 83 that is relevant for our
discussion. Defining new variables on the links of the original square lattice

$\displaystyle\hat{E}_{j,j+\hat{x}}$ $\displaystyle=$
$\displaystyle\nabla_{y}\hat{\Psi}_{\bar{j}},\ \
\hat{E}_{j,j+\hat{y}}=-\nabla_{x}\hat{\Psi}_{\bar{j}},$ (86)
$\displaystyle\hat{\phi}_{\bar{j}}$ $\displaystyle=$
$\displaystyle(\nabla\times\hat{A}_{j,j+\mu})_{\text{around }\bar{j}}.$ (88)

The electric fields $\hat{E}$ are forced to obey the usual Gauss’ law
constraint of $\text{QED}_{3}$,

$\displaystyle G_{j}=(\nabla_{x}\hat{E}_{x}+\nabla_{y}\hat{E}_{y})_{j}\equiv
0$ (89)

The dual Hamiltonian, consequently is a fine-tuned version of the compact
$\text{QED}_{3}$

$\displaystyle
H=-K\sum_{\square}\cos(\nabla\times\hat{A})_{\square}+\frac{U}{8}\left(\nabla_{x}\hat{E}_{x}-\nabla_{y}\hat{E}_{y}\right)^{2}.$
(90)

To bridge the connection between the two representations of the XY-plaquette
model, i.e. Eq. 78 and Eq. 90, we need to fix the gauge of $\hat{A}$. We
choose a gauge such that $\nabla_{y}\hat{A}_{x}+\nabla_{x}\hat{A}_{y}=0$,
which can always be met by choosing a right gauge transformation
$\hat{A}\rightarrow\hat{A}+\nabla\varphi$. For example, for a configuration of
$\hat{A}$, one can choose function $\varphi(x,y)$ to be
$-\frac{1}{2}\int_{0}^{x}\int_{0}^{y}dx^{\prime}dy^{\prime}(\nabla_{x}\hat{A}_{y}+\nabla_{y}\hat{A}_{x})$.
Then the relation between Eq. 78 and Eq. 90 is that
$\nabla_{x}\hat{\theta}=2\hat{A}_{x}$, $\nabla_{y}\hat{\theta}=-2\hat{A}_{y}$,
$2\hat{n}=\nabla_{x}\hat{E}_{x}-\nabla_{y}\hat{E}_{y}$. Note that since
$\nabla_{x}\hat{E}_{x}=-\nabla_{y}\hat{E}_{y}$ due to the Gauss law
constraint, $\hat{n}$ takes its value in all integers. The subsystem
transformation $\hat{\theta}(x)\rightarrow\hat{\theta}(x)+f_{1}(x)+f_{2}(y)$
becomes part of the residual gauge transformation:
$\hat{A}\rightarrow\hat{A}+\nabla\varphi$, where
$\varphi=(f_{1}(x)-f_{2}(y))/2$. In particular, the nonzero four point
correlation function of $\cos(\mathrm{i}\hat{\theta})$ where $\hat{\theta}$
locate on four corners of a rectangle becomes a Wilson loop of the gauge field
$\hat{A}$.


# A new goodness of fit test for uniform distribution with censored
observations

Sudheesh K Kattumannila and Sreedevi E Pb,†
aIndian Statistical Institute, Chennai, India,
bSNGS College, Pattambi, India

###### Abstract.

Using fixed point characterization, we develop a new goodness of fit test for
uniform distribution. We also discuss how the right censored observations can
be incorporated in the proposed test procedure. We study the asymptotic
properties of the proposed test statistics. A Monte Carlo simulation is
carried out to evaluate the finite sample performance of the tests. We
illustrate the test procedures using real data sets.
Key words: Right censoring; Stein’s identity; U-statistics.

Corresponding author email<EMAIL_ADDRESS>
##  1\. Introduction

Uniform distribution is a widely used statistical model in various fields of
applied as well as theoretical statistics. Due to the range property of
standard uniform distribution, many Monte Carlo simulation algorithms use a
sample from uniform distribution to generate random samples from other
distributions. In view of the probability integral transformation, the simple
goodness of fit problem of testing that a sample is from a particular
continuous distribution is equivalent to testing that the transformed sample
is from a uniform distribution on the interval $(0,1)$. This motivates us to
develop a goodness of fit test for uniform distribution.

Tests for uniform distribution has great significance in all discipline due to
the property that the distribution with maximum entropy is standard uniform.
Various tests were developed in literature to test the hypothesis of
uniformity. We refer interested readers to Kimball (1947), Sherman (1950),
Quesenberry and Miller Jr. (1977), Hegazy and Green (1975), Young (1982),
Cheng and Spiring (1987) and Frozini et al. (1987). Under the uniform
distribution assumption, Asgharian et al. (2002) estimated survival of
dementia patients and de Una-Alvarez (2004) estimated the distribution of
unemployment spells of women. In many situations, we deal with a uniform
distribution with truncation and/or censoring. A test for goodness of fit for
a uniform truncation model was developed by Mandel and Betensky (2007).
Recently, Cho et al. (2021) developed a test for uniformity of exchangeable
random variables on the circle. However, goodness of fit for uniform
distribution under right censored case is not developed yet. Using Stein’s
type characterisation we develop goodness of fit tests for uniform
distribution with complete as well as censored data.

Stein’s identity and its application has been studied in literature. See
Sudheesh (2009) and Sudheesh and Isha (2016) for a generalized Stein’s
identity and its application. Making use of Stein’s type identity, Ebner and
Liebenberg (2020) obtain the following fixed point characterization for beta
distribution.

######  Theorem 1.

Let $X$ be a random variable (r.v.) taking values in $[0,1]$. Then $X\sim
B(\alpha,\beta)$ for $\alpha,\beta>0$ if and only if

$(\alpha+\beta)E(XI(X>t))=\alpha
E(I(X>t))+\frac{t^{\alpha}(1-t)^{\beta}}{B(\alpha,\beta)},\,\,0\leq t\leq 1,$

where $I(A)$ denotes the indicator function of a set $A$.

######  Corollary 1.

Let $X$ be a r.v. taking values in $[0,1]$. Then $X\sim U(0,1)$ if and only if

$\vspace{-0.1in}2E(XI(X>t))=E(I(X>t))+{t(1-t)},\,\,0\leq t\leq 1.$

Using the fixed point characterization given in Corollary 1, we develop a
goodness of fit test for uniform distribution. The rest of the paper is
organised as follows. In Section 2, we develop a new non-parametric test for
uniform distribution. We study the asymptotic properties of the test
statistic. In Section 3, we discuss how to incorporate right censored
observations in the proposed methodology. The results of the Monte Carlo
simulation study are reported in Section 4. In Section 5, the test procedures
are illustrated with application on real data sets. Finally, in Section 6, we
give the concluding remarks.

##  2\. Test statistic

In this section, we develop a goodness of fit test for uniform distribution
for complete data. Based on a random sample $X_{1},...,X_{n}$ from $U(0,1)$,
we are interested in testing the null hypothesis

$H_{0}:F\text{ has uniform distribution}.$

against

$H_{1}:F\text{ does not follow uniform distribution}.$

For this purpose, we define a departure measure which discriminate between
null and alternative hypothesis. Consider $\Delta(F)$ given by

$\displaystyle\Delta(F)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\left(2E(XI(X>t))-E(I(X>t))-{t(1-t)}\right)dF(t).$

In view of Corollary 1, $\Delta(F)$ is zero under $H_{0}$ and non zero under
$H_{1}$. Hence $\Delta(F)$ can be considered as a measure of departure from
the null hypothesis $H_{0}$ towards the alternative hypothesis $H_{1}$.

Next we express $\Delta(F)$ in a simple form as

$\displaystyle\Delta(F)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\left(2E(XI(X>t))-E(I(X>t))-t(1-t)\right)dF(t)$ (1)
$\displaystyle=$ $\displaystyle
2\int_{0}^{1}\int_{0}^{1}yI(y>t)dF(y)dF(t)-E(X)-E(X(1-X))$ $\displaystyle=$
$\displaystyle 2\int_{0}^{1}\int_{t}^{1}ydF(y)dF(t)-2E(X)+E(X^{2})$
$\displaystyle=$ $\displaystyle
2\int_{0}^{1}y\int_{0}^{y}dF(t)dF(y)-2E(X)+E(X^{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{1}2yF{(y)}dF(y)-2E(X)-E(X^{2})$ $\displaystyle=$
$\displaystyle E\big{(}max(X_{1},X_{2})-2X+X^{2}\big{)}.$

Consider a symmetric kernel
$h_{1}(X_{1},X_{2})=\frac{1}{2}(2\max(X_{1},X_{2})-2X_{1}-2X_{2}+X_{1}^{2}+X_{2}^{2})$,
then $E(h_{1}(X_{1},X_{2}))=\Delta(F)$. Hence we propose a U-statistic based
test statistic given by

$\widehat{\Delta}=\frac{2}{n(n-1)}\sum_{i=1}^{n}\sum_{j=1,j<i}^{n}h_{1}(X_{i},X_{j}),$

for testing uniform distribution. The $\widehat{\Delta}$ is an unbiased and
consistent estimator of $\Delta(F)$ (Lehmann, 1951). Let
$X_{(i)},\,\,i=1,\ldots,n$ be the order statistics based on $n$ independently
and identically observations from $U(0,1)$. Then we can represent
$\widehat{\Delta}$ in a simple form as

$\widehat{\Delta}=\frac{1}{n(n-1)}\sum_{i=1}^{n}\left(2(i-n)+(n-1)X_{(i)}\right)X_{(i)}.$

We reject the null hypothesis $H_{0}$ against the alternative for large value
of $|\widehat{\Delta}|$. Next we find an asymptotic critical region of the
test using normal approximation. The following result is useful in this
direction.

######  Theorem 2.

As $n\rightarrow\infty$, $\sqrt{n}(\widehat{\Delta}-\Delta)$ converges in
distribution to normal with mean zero and variance $\sigma^{2}$, where
$\sigma^{2}$ is given by

$\sigma^{2}=Var\Big{(}2XF(X)+2\int_{X}^{1}ydF(y)-2X+X^{2}\Big{)}.$ (2)

Proof: Asymptotic normality of $\widehat{\Delta}$ follows from the central
limit theorem for U-statistics. The asymptotic variance is $4\sigma_{1}^{2}$
where $\sigma_{1}^{2}$ is given by (Lee, 1990)

$\sigma_{1}^{2}=Var\left[E\left(h(X_{1},X_{2}|X_{1})\right)\right].$ (3)

Consider

$\displaystyle E[\max(X_{1},X_{2})|X_{1}=x]$ $\displaystyle=$ $\displaystyle
E[xI(X_{2}\leq x)+X_{2}I(X_{2}>x)]$ $\displaystyle=$ $\displaystyle
xP(X_{2}\leq x)+\int_{0}^{1}yI(x<y)dF(y)$ $\displaystyle=$ $\displaystyle
xF(x)+\int_{x}^{1}ydF(y).$

Hence, from (3) we obtain the variance expression given in (2) and the proof
of the theorem follows.
Under the null hypothesis $H_{0}$, we know that $\Delta{(F)}=0$. Hence we have
the following corollary.

######  Corollary 2.

Under $H_{0}$, as $n\rightarrow\infty$, $\sqrt{n}\widehat{\Delta}$ converges
in distribution to normal with mean zero and variance $\frac{1}{45}$.

Proof: Under $H_{0}$, we have

$\displaystyle\sigma^{2}$ $\displaystyle=$ $\displaystyle
Var\Big{(}2XF(X)+2\int_{X}^{1}ydF(y)-2X+X^{2}\Big{)}$ $\displaystyle=$
$\displaystyle Var\Big{(}2X^{2}+2\int_{X}^{1}ydy-2X+X^{2}\Big{)}$
$\displaystyle=$ $\displaystyle 4Var\Big{(}X(X-1)\Big{)}=\frac{1}{45}.$

A distribution free test for testing uniform distribution can be constructed
using Corollary 1. We reject the null hypothesis $H_{0}$ against the
alternative hypothesis $H_{1}$ at a significance level $\alpha$, if

${\sqrt{45n}|\widehat{\Delta}|}>Z_{\alpha/2},$

where $Z_{\alpha}$ is the upper $\alpha$-percentile point of the standard
normal distribution. The performance of the test is evaluated through Monte
Carlo simulation study and the result of the same is reported in Section 4.

##  3\. Test for right censored case

Next we discuss how the censored observations can be incorporated in the
proposed testing procedure. Consider the right-censored data $(Y,\delta)$,
with $Y=min(X,C)$ and $\delta=I(X\leq C)$, where $C$ is the censoring time. We
assume censoring times and lifetimes are independent. Now we need to address
the testing problem discussed in Section 2 based on $n$ independent and
identical observations $\\{(Y_{i},\delta_{i}),1\leq i\leq n\\}$. We develop
the test using the same departure measure $\Delta(F)$ given in (1). Consider a
symmetric kernel defined by
$h(Y_{1},Y_{2})=\frac{1}{2}(2\max(Y_{1},Y_{2})-2Y_{1}-2Y_{2}+Y_{1}^{2}+Y_{2}^{2})$.
An estimator of (1) in the right censored case is given by (Datta et al.,
2010)

$\widehat{\Delta}_{c}=\frac{2}{n(n-1)}\sum_{i=1}^{n}\sum_{j<i;j=1}^{n}\frac{h(Y_{1},Y_{2})\delta_{i}\delta_{j}}{\widehat{K}_{c}(Y_{i})\widehat{K}_{c}(Y_{j})},$

provided $\widehat{K}_{c}(Y_{i})>0$ and $\widehat{K}_{c}(Y_{j})>0$, with
probability 1 where $\widehat{K}_{c}$ is the Kaplan-Meier estimator of
$K_{c}$, the survival function of $C$. We reject $H_{0}$ in favour of $H_{1}$
for large values of $|\widehat{\Delta}_{c}|$.

Next we obtain the limiting distribution of $\widehat{\Delta}_{c}$. Let
$N_{i}^{c}(t)=I(Y_{i}\leq t,\delta_{i}=0)$ be the counting process corresponds
to the censoring variable $C_{i}$, $R_{i}(t)=I(Y_{i}\geq t)$. Also let
$\lambda_{c}$ be the hazard rate of $C$. The martingale associated with this
counting process $N_{i}^{c}(t)$ is given by

$M_{i}^{c}(t)=N_{i}^{c}(t)-\int_{0}^{t}R_{i}(u)\lambda_{c}(u)du.$

Denote $h_{1}(x)=E(h(X_{1},X_{2}|X_{1}=x))$, $G(x,y)=P(X_{1}\leq x,Y_{1}\leq
y,\delta=1),x\in\mathcal{X}$, $\bar{H}(t)=P(Y_{1}>t)$ and

$w(t)=\frac{1}{\bar{H}(t)}\int_{\mathcal{X}\times[0,\infty)}{\frac{h_{1}(x)}{K_{c}(y-)}I(y>t)dG(x,y)}.$

######  Theorem 3.

If $E(Y_{1}^{2})<\infty$ and the integrals
$\int_{\mathcal{X}\times[0,\infty)}{\frac{h^{2}_{1}(x)}{K_{c}^{2}(y)}dG(x,y)}$
and $\int_{0}^{\infty}w^{2}(t)\lambda_{c}(t)dt$ are finite, then as
$n\rightarrow\infty$, $\sqrt{n}(\widehat{\Delta}_{c}-\Delta)$ converges in
distribution to Gaussian with mean zero and variance $4\sigma_{c}^{2}$, where
$\sigma_{c}^{2}$ is given by

$\sigma_{c}^{2}=Var\left(\frac{\Big{(}2XF(X)+2\int_{X}^{1}ydF(y)-2X+X^{2}\Big{)}\delta_{1}}{4K_{c}(Y_{1}-)}+\int
w(t)dM_{1}^{c}(t)\right).$

######  Corollary 3.

Assume that the condition stated in Theorem 3 holds. Under $H_{0}$, as
$n\rightarrow\infty$, the distribution of $\sqrt{n}\widehat{\Delta}$ is
Gaussian with mean zero and variance $4\sigma_{c0}^{2}$, where
$\sigma_{c0}^{2}$ is given by

$\sigma_{c0}^{2}=Var\Big{(}\frac{X(X-1)\delta_{1}}{K_{c}(Y_{1}-)}+\int
w(t)dM_{1}^{c}(t)\Big{)}.$

Next we find an estimator of $\sigma_{c0}^{2}$ using the reweighed techniques.
An estimator of $\sigma_{c0}^{2}$ is given by

$\widehat{\sigma}_{c0}^{2}=\frac{4}{(n-1)}\sum_{i=1}^{n}(V_{i}-\bar{V})^{2},$

where

$V_{i}=\frac{\widehat{h}_{1}(X_{i})\delta_{i}}{\widehat{K}_{c}(Y_{i})}+\widehat{w}(X_{i})(1-\delta_{i})-\sum_{j=1}^{n}\frac{\widehat{w}(X_{i})I(X_{i}>X_{j})(1-\delta_{i})}{\sum_{i=1}^{n}I(X_{i}>X_{j})},$
$\bar{V}=\frac{1}{n}\sum_{i=1}^{n}V_{i},\quad\widehat{h}_{1}(X)=\frac{1}{n}\sum_{i=1}^{n}\frac{h(X,Y_{i})\delta_{i}}{\widehat{K}_{c}(Y_{i}-)},\quad
R(t)=\frac{1}{n}\sum_{i=1}^{n}I(Y_{i}>t)$

and

$\widehat{w}(t)=\frac{1}{R(t)}\sum\limits_{i=1}^{n}\frac{\widehat{h}_{1}(X_{i})\delta_{i}}{\widehat{K}_{c}(Y_{i})}I(X_{i}>t).$

Under right censored situation, we reject the null hypothesis $H_{0}$ against
the alternative hypothesis $H_{1}$ at a significance level $\alpha$, if

$\frac{\sqrt{n}|\widehat{\Delta}_{c}|}{\widehat{\sigma}_{0c}}>Z_{\alpha/2}.$

The results of the Monte Carlo Simulation which asses the finite sample
performance of the test is also reported in Section 4.

##  4\. Empirical evidence

To evaluate the finite sample performance of our test procedures, we conduct a
Monte Carlo simulation study using R software. The simulation is repeated ten
thousand times. For complete data, to show the competitiveness of our test, we
compare the empirical power our test with the existing test procedures.

###  4.1. Uncensored case

We find the empirical type I error and power of the proposed tests for
different choices of alternatives. We compare our new test with the classical
procedures based on the empirical distribution function Kolmogorov-Smirnov
test and the well used tests for uniform distribution proposed by Sherman
(1950), Quesenberry and Miller Jr. (1977) and Frozini et al. (1987). Next we
briefly discuss about these tests. The Kolmogorov-Smirnov test statistic is
given by $KS=\text{max}\\{D^{+},D^{-}\\}$ where

$D^{+}=\max_{j=1,2,...,n}\left(\frac{j}{n}-\widehat{F}(X_{(j)})\right)~{}~{}~{}~{}\text{and}~{}~{}~{}~{}D^{-}=\max_{j=1,2,...,n}\left(\widehat{F}(X_{(j)})-\frac{j-1}{n}\right)$

with $\widehat{F}(.)$ as the empirical distribution function of the observed
data. Frozini’s test statistic has the form

$F=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\left|X_{(j)}-\frac{j-0.5}{n}\right|.$

Sherman test statistic is given by

$S=\frac{1}{2}\sum_{j=1}^{n+1}\left|(X_{(j)}-X_{(j-1)})-\frac{1}{n+1}\right|.$

with $X_{(0)}=0$ and $X_{(n+1)}=1$. The Q-Statistic proposed by Quesenberry
and Miller Jr. (1977) is given by

$Q=\sum_{j=1}^{n+1}\left(X_{(j)}-X_{(j-1)}\right)^{2}+\sum_{j=1}^{n}(X_{(j+1)}-X_{(j)})(X_{(j)}-X_{(j-1)}).$

Table 1. Empirical type I error of the tests

t1 | $\widehat{\Delta}$ | $KS$ | $F$ | $S$ | $Q$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
25 | 0.0108 | 0.0546 | 0.0113 | 0.0486 | 0.0108 | 0.0513 | 0.0109 | 0.0516 | 0.0085 | 0.0475
50 | 0.0106 | 0.0534 | 0.0108 | 0.0491 | 0.0094 | 0.0464 | 0.0088 | 0.0485 | 0.0113 | 0.0483
75 | 0.0104 | 0.0492 | 0.0096 | 0.0509 | 0.0104 | 0.0474 | 0.0095 | 0.0498 | 0.0108 | 0.0511
100 | 0.0102 | 0.0502 | 0.0103 | 0.0508 | 0.0102 | 0.0486 | 0.0095 | 0.0502 | 0.0096 | 0.0504

Table 2. Empirical power the tests: $U(0,1.2)$

| $\widehat{\Delta}$ | $KS$ | $F$ | $S$ | $Q$ |
---|---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level |
25 | 0.5335 | 0.6950 | 0.2374 | 0.3221 | 0.2643 | 0.4635 | 0.6012 | 0.7182 | 0.6721 | 0.7682 |
50 | 0.8146 | 0.9068 | 0.3282 | 0.5927 | 0.4956 | 0.7171 | 0.8932 | 0.9321 | 0.8934 | 0.9302 |
75 | 0.9285 | 0.9682 | 0.5599 | 0.7997 | 0.6682 | 0.8562 | 0.9431 | 0.9732 | 0.9473 | 0.9651 |
100 | 0.9882 | 0.9921 | 0.7481 | 0.9212 | 0.8036 | 0.9387 | 0.9865 | 0.9921 | 0.9832 | 0.9972 |

Table 3. Empirical power the tests: Exp$(1)$

| $\widehat{\Delta}$ | $KS$ | $F$ | $S$ | $Q$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
25 | 0.9998 | 0.9999 | 0.7744 | 0.8925 | 0.9939 | 0.9981 | 0.9982 | 1.0000 | 0.9999 | 1.0000
50 | 1.0000 | 1.0000 | 0.9851 | 0.9976 | 0.9998 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
75 | 1.0000 | 1.0000 | 0.9999 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
100 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000

Table 4. Empirical power the tests: Gamma$(1,2)$

| $\widehat{\Delta}$ | $KS$ | $F$ | $S$ | $Q$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
25 | 0.9999 | 1.0000 | 0.7841 | 0.8981 | 0.9994 | 0.9988 | 0.9999 | 1.0000 | 0.9999 | 1.0000
50 | 1.0000 | 1.0000 | 0.9863 | 0.9974 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
75 | 1.0000 | 1.0000 | 0.9997 | 0.9999 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
100 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000

Table 5. Empirical power the tests: Weibull$(1,2)$

| $\widehat{\Delta}$ | $KS$ | $F$ | $S$ | $Q$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
25 | 0.9999 | 1.0000 | 0.7789 | 0.8982 | 0.9943 | 0.9983 | 0.9999 | 1.0000 | 0.9999 | 1.0000
50 | 1.0000 | 1.0000 | 0.9846 | 0.9980 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
75 | 1.0000 | 1.0000 | 0.9999 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
100 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000

Table 6. Empirical power the tests: Pareto$(1,1)$

| $\widehat{\Delta}$ | $KS$ | $F$ | $S$ | $Q$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
25 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
50 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
75 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
100 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000

For finding the empirical type I error, we generated samples of sizes
$n=25,50,75$ and $100$ from $U(0,1)$. Empirical type I error of the proposed
test and above given tests for uniform distribution is given in Table 1. We
consider different choices of alternatives for finding empirical power of the
tests and the results of the simulation study are given in Tables 2-6.

The comparative study shows the competitiveness of the newly proposed test to
classical procedures. Type I error of the test attains the nominal level and
the test shows very good power. For the choice of alternative distribution
$U(0,1.2)$, the test has good power which shows that our test captures even a
slight deviation from the null distribution. Also, even for small sample size
$n=25$, we obtain high power which reaches to unity by $n=50$ for other
choices of distributions.

###  4.2. Censored case

The performance of proposed test procedures incorporating right censored
observations is also validated through a Mote Carlo simulation study. We
consider a mild censoring situation where 20% of lifetimes are censored and
heavy censored situation where 40% of the lifetimes are censored. In censored
case, we choose sample sizes $n=50,75,100$ and $200$ and simulation is
repeated ten thousand times.

For finding empirical type I error we simulate observations from $U(0,1)$.
Various choices of alternative are consider for finding empirical power. In
both these cases, we generate censored observations from $U(0,c)$ where $c$ is
chosen such way that $P(X>C)=0.2(0.4)$. Results of the simulation study with
right censored observations are given in Tables 7 and 8. Type I error reaches
the chosen significance level in both mild and heavy censored situations. As
in case for complete data, the test statistic points out even a small
deviation ($U(0,1.2)$) from the null hypothesis. Power of the test reaches one
for other choices of alternative distributions even for sample size $n=50$.

Table 7. Empirical type I error and power the test with 20% censoring

| $U(0,1)$ | $U(0,1.2)$ | $Exp(1)$ | $Weibull(1,2)$ | $Pareto(1,1)$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
50 | 0.0088 | 0.0474 | 0.6993 | 0.8303 | 0.9999 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
75 | 0.0096 | 0.0506 | 0.8635 | 0.9387 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
100 | 0.0106 | 0.0489 | 0.9397 | 0.9761 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000
200 | 0.0106 | 0.0499 | 0.9984 | 0.9995 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000

Table 8. Empirical type I error and power the test for 40% censoring

| $U(0,1)$ | $U(0,1.2)$ | $Exp(1)$ | $Weibull(1,2)$ | $Pareto(1,1)$
---|---|---|---|---|---
$n$ | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level | 1% level | 5% level
50 | 0.0079 | 0.0469 | 0.5174 | 0.6892 | 0.9771 | 0.9856 | 1.0000 | 1.0000 | 1.0000 | 1.0000
75 | 0.0091 | 0.0488 | 0.6537 | 0.8112 | 0.9972 | 0.9982 | 1.0000 | 1.0000 | 1.0000 | 1.0000
100 | 0.0109 | 0.0508 | 0.8025 | 0.9209 | 0.9996 | 0.9998 | 1.0000 | 1.0000 | 1.0000 | 1.0000
200 | 0.0103 | 0.0502 | 0.9839 | 0.9963 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000

##  5\. Data analysis

###  5.1. Uncensored case

The newly proposed test procedures are applied to a real data set for
illustration. We consider the data set discussed in Illowsky and Dean (2018)
in Page 317, Table 5.1. The data set consist of smiling times of 55 babies
measured in seconds. The data originally follows a uniform distribution
$U(0,23)$. Using the proposed procedure, we obtain $\hat{\Delta}$ as 8157.285.
The test static value lies in critical region and we reject the hypothesis
that the data follows $U(0,1)$. Next, to ensure the validity of the proposed
procedure, we standardize the same data to $U(0,1)$. For this transformed data
the value of the test statistic is obtained as $\hat{\Delta}=-1.2523$, which
belongs to the acceptance region. Hence we accept the null hypothesis that the
data follows $U(0,1)$.

###  5.2. Censored case

We consider a real data set on life times of 10 pieces of equipment installed
in a system. We consider the data given in Table 1.2 of Page 4 in Lawless
(2011). Out of the 10 observed lifetimes 3 are randomly right censored. We
analysed the data using the procedures developed for censored data. The value
of test statistic is obtained as $\hat{\Delta}_{c}=9858.522$, which clearly
indicates that the data does not follow uniform distribution. We then
standardize the same data to $U(0,1)$. For the transformed data we obtain
$\hat{\Delta}_{c}=-0.3451$ and we accept the null hypothesis that the data
comes from $U(0,1)$.

##  6\. Concluding Remarks

Using the recently introduced Stein’s characterization, a simple non-
parametric test based on U-statistic theory is developed for testing uniform
distribution. The test is distribution free. We study the asymptotic
properties of the test statistic. Even though several tests are available in
literature to test the uniformity, none of these tests incorporated censored
samples. We discussed how to incorporate censored data in our methodology. An
extensive Monte Carlo simulation study is carried to validate the finite
sample performance of the tests procedures. Power comparisons show that the
performance of our test is competent with the existing tests in the complete
case. Also the test has well controlled error rate even for small sample
sizes. In censored case, even with high percentage of censored data (40%) our
test performs well in terms of empirical power and attains the size of the
test.

## References

  * [1] Asgharian, M., M’Lan, C. E. and Wolfson, D. B. (2002). Length-biased sampling with right censoring: An unconditional approach. Journal of the American Statistical Association, 97, 201–209.
  * [2] Cho, S., Choi, Y. G., Lim, J., Lee, W. J., Bai, H. J. and Kwon, S. (2021). Test for uniformity of exchangeable random Variables on the circle. Journal of the Korean Statistical Society, 50, 615–630.
  * [3] Cheng, S. W. and Spiring, F. A. (1987). A test to identify the uniform distribution, with applications to probability plotting and other distributions. IEEE transactions on reliability, 36, 98–105.
  * [4] Datta, S., Bandyopadhyay, D. and Satten, G. A. (2010). Inverse probability of censoring weighted U‐statistics for right‐censored data with an application to testing hypotheses. Scandinavian Journal of Statistics, 37, 680–700.
  * [5] de Una-Alvarez, J. (2004). Nonparametric estimation under length-biased sampling and type I censoring: A moment based approach. Annals of the Institute of Statistical Mathematics, 56, 667–681
  * [6] Ebner, B. and Liebenberg, S. C. (2020). On a new test of fit to the beta distribution. Stat, e341.
  * [7] Frozini, B. V., Revesz, P. and Sarkadi, K. (1987). On the distribution and power of goodness-of-fit statistic with parametric and nonparametric applications,“Goodness-of-fit”. Amsterdam-Oxford-New York: North-Holland Publ. Comp, 133.
  * [8] Hegazy, Y. A. S. and Green, J. R. (1975). Some new goodness‐of‐fit tests using order statistics. Journal of the Royal Statistical Society: Series C, 24, 299–308.
  * [9] Illowsky, B. and Dean, S. (2018). Introductory Statistics, OpenStax, Houston.
  * [10] Kimball, B. F. (1947). Some basic theorems for developing tests of fit for the case of the non-parametric probability distribution function, I. The Annals of Mathematical Statistics, 18, 540–548.
  * [11] Lawless, J. F. (2011).  Statistical Models and Methods for Lifetime Data. John Wiley & Sons, New York.
  * [12] Lehmann, E. L. (1951). Consistency and unbiasedness of certain nonparametric tests. The Annals of Mathematical Statistics, 22, 165–179.
  * [13] Lee, A. J. (1990). U-Statistics: Theory and Practice, Marcel Dekker Inc., New York.
  * [14] Mandel, M. and Betensky, R. A. (2007). Testing goodness of fit of a uniform truncation model. Biometrics, 63, 405–412.
  * [15] Quesenberry, C. P. and Miller Jr, F. L. (1977). Power studies of some tests for uniformity. Journal of Statistical Computation and Simulation, 5, 169–191.
  * [16] Sherman, B. (1950). A random variable related to the spacing of sample values. The Annals of Mathematical Statistics, 21, 339–361.
  * [17] Sudheesh, K. K. (2009). On Stein’s identity and its applications. Statistics and Probability Letters, 79, 1444–1449.
  * [18] Sudheesh, K. K. and Dewan, I. (2016). On generalized moment identity and its applications: A unified approach. Statistics, 50, 1149–1160.
  * [19] Young, D. L. (1982). The linear nearest neighbour statistic. Biometrika, 69, 477–480.