here, this means that the completeness of paths of accelerated observers, in
addition to the geodesic ones, must be guaranteed for a given space-time to
give an accurate description of the physical world. This becomes yet another
test for the regularity of a given space-time: to verify that every non-
geodesic observer with arbitrary motion and, in particular, with an arbitrary
(but bound) acceleration finds infinitely extended paths towards the future
and the past. For instance, there are examples (admittedly, toy-like) in which
a geodesically complete space-time contains time-like curves with (bound)
acceleration but finite total proper length Geroch:1968ut , and one would like
to make sure that such a thing does not occur in their beloved “regular”
space-times.

Under the presence of external forces, the right-hand side of the geodesic
equation (121) for a massive particle no longer vanishes:

$\frac{Du^{\mu}}{d\tau}=\frac{d^{2}x^{\mu}}{d\tau^{2}}+\Gamma_{\alpha\beta}^{\mu}\frac{dx^{\alpha}}{d\tau}\frac{dx^{\beta}}{d\tau}=a^{\mu}\neq
0$ (135)

where $\tau$ is the proper time of a time-like observer. The acceleration
undergone by our observer can be computed by considering an orthonormal
Frenet-Serret basis $\\{u_{(0)},u_{(1)},u_{(2)},u_{(3)}\\}$ Friedman:2015ufz
(parenthesis are introduced here to stress the coordinate-free character of
these vectors) to the tangent space satisfying the equation Olmo:2017fbc

$\frac{Du_{(a)}}{ds}=u_{(b)}{A^{(b)}}_{(a)}\hskip 5.69046pt;\hskip
5.69046pt{A^{(b)}}_{(a)}=\left(\begin{array}[]{cccc}0&k(s)&0&0\\\
k(s)&0&\tau_{1}(s)&0\\\ 0&-\tau_{1}(s)&0&\tau_{2}(s)\\\ 0&0&-\tau_{2}(s)&0\\\
\end{array}\right)\ ,$ (136)

where the curvature $k(u)$ is interpreted as the linear acceleration
experienced by the observer in the direction of $u_{(1)}$, while the first
$\tau_{1}(u)$ and the second $\tau_{2}(u)$ torsions of the curve correspond to
the rotational acceleration $\sqrt{\tau_{1}^{2}+\tau_{2}^{2}}$ along the axis
given by the vector $\omega=\tau_{2}\lambda_{(1)}+\tau_{1}\lambda_{(3)}$.
Obviously, when there is neither linear acceleration nor torsion, Eq.(135)
recovers the usual geodesic equation. In a spherically symmetric space-time
written in suitable coordinates as
$ds^{2}=-A(y)dt^{2}+A^{-1}(y)dy^{2}+r^{2}(y)d\Omega^{2}$, making use of the
conservation of energy and angular momentum, and restricting ourselves to
linear accelerations only for simplicity, Eq.(135) can be written, after a
straightforward algebraic exercise, as Olmo:2017fbc

$(u^{y})^{2}+A(y)\left(1+\frac{L^{2}}{r^{2}(y)}\right)=\left[E+\int^{y}_{y_{0}}\frac{k(s)dy^{\prime}}{\sqrt{1+\frac{L^{2}}{r^{2}(y^{\prime})}}}\right]^{2}\
.$ (137)

where $u^{y}=dy/du$ and $k(y)$ is the total acceleration of the curve in this
simplified model. This is the natural extension of the Eq.(122) to accelerated
observers, covering every space-time considered in this chapter. We can thus
check the behaviour of any such set of accelerated observers in any pick of
geodesically complete sample found above, and for convenience we shall
consider (once again) the quadratic/EiBI gravity by its agreeable properties.
By considering a test charged (Maxwell) particle which experiences an
acceleration (in modulus) of size $k(y)=\frac{Qq}{mr^{2}(y)}$, and switching
by convenience to the variable $x$ defined by the relation $dy=dx/\Omega_{+}$,
then Eq.(137) becomes

$\left(\frac{u^{x}}{\Omega_{+}}\right)^{2}+A(x)\left(1+\frac{L^{2}}{r^{2}(x)}\right)=\left(\tilde{E}+I^{BI}_{L}(r)\right)^{2}\
,$ (138)

where $\tilde{E}=E-I^{BI}_{L}(r_{i})$, with $r_{i}$ the initial location of
the particle and the function

$I^{BI}_{L}(r)=\frac{Qq}{m}\int^{r}\frac{dr}{r\Omega_{-}^{\frac{1}{2}}\sqrt{r^{2}+L^{2}}}=-\frac{Qq}{m}\frac{\text{EllipticF}\left[\sin^{-1}\left(\sqrt{\frac{L^{2}+r_{c}^{2}}{L^{2}+r^{2}}}\right),\frac{L^{2}-r_{c}^{2}}{L^{2}+r_{c}^{2}}\right]}{\sqrt{r_{c}^{2}+L^{2}}}\
.$ (139)

encodes the acceleration term in this scenario. Let us recall the constant
character of the radial function at the bounce $x=0$ ($r=r_{c}$), which in
turn makes the function $I^{BI}_{l}(r_{c})$ to become a constant there and,
therefore, the single effect induced by the acceleration of our point-like
charge is to shift its total energy, therefore having no effect in the
extensibility of geodesics. We thus conclude the completeness of all null and
time-like paths (in the latter both geodesic and accelerated) in these
geometries.

Let us note that in those cases in which the bounce lies at the future (past)
boundary of the space-time (quadratic $f(\mathcal{R})$ case coupled to
different fluids), one expects that the introduction of bounded accelerations
will not alter the conclusion on every time-like geodesic being bounced back
by the infinite potential barrier. And what if, by the sake of curiosity, we
consider an unbounded acceleration? In such a case, a similar computation as
above reveals that such observers would be able to get to the bounce location
and gain access to the other side of the wormhole. Such scenarios are, of
course, unlikely to be physically achievable.

###  4.6 Further tests?

Up to now we have considered tests with idealized point-like particles
(completeness of geodesic and of accelerated paths), and the effects of tidal
forces from unbound curvatures upon extended observers. Is there any other
further test one could bring to the table?. There is, of course, the important
and physically relevant case of the interaction of test waves with potentially
pathological space-time regions. In a regular geometry it seems natural to
demand that the propagation of any wave through it should be well defined at
all times, i.e., that one can establish a well-posed problem of transmission
and reflection coefficients. This includes, in particular, the case of
gravitational waves, where we find the notable (and perhaps non-intuitive)
fact that in Palatini theories of gravity, these waves propagate on the
effective gravitational geometry provided by $q_{\mu\nu}$ rather than the
space-time one, see BeltranJimenez:2017uwv  for details. For the sake of this
chapter, however, we shall simplify the problem down to scalar waves
propagating upon the latter geometry, a scenario in which relevant insights on
the regularity of these space-times can be obtained.

Let us thus consider a massive scalar field with equation of motion
$(\Box-m^{2})\phi=0$. Using the usual mode decomposition, $\phi_{\omega
ml}=e^{-i\omega t}Y_{lm}(\theta)f_{\omega l}(x)/r(x)$, one can find the
expression of the equation for its radial component in any space-time of
interest. Considering again the quadratic/EiBI gravity case, and formulating
the problem in terms of the coordinate $y$ given by $dy/dx=1/(A\Omega_{+})$,
one finds a Schrödinger-like equation of the form Olmo:2015dba

$-f_{yy}+V_{eff}(y)f=\omega^{2}f\hskip 2.84544pt;\hskip
2.84544ptV_{eff}=\frac{r_{yy}}{r}+A(r)\left(m^{2}+\frac{l(l+1)}{r^{2}}\right)$
(140)

One can thus proceed to define an incoming wave from past null infinity (for
horizonless solutions), or from the event/Cauchy horizon (for black holes) and
study its behaviour as the bounce location, $r=r_{c}$, is approached. This
problem is better captured by the computation of the transmission and
reflection factors, as well as with the cross section of the former. The final
result is that no pathologies over these factors, or any other illness of
these space-times, are found Olmo:2015dba . This implies that, in addition to
regular black holes outside their event horizons, which is the relevant case
for gravitational wave astronomy, within these space-times one can also
consider the propagation of waves over regular naked wormholes, thus naturally
connecting the theoretical issue of space-time singularities with the
astrophysically-relevant one of gravitational wave imprints from these
objects.

##  5 Closing thoughts

Well inside the twentieth-first century, space-time singularities are still a
cornerstone upon which many attempts to extend GR are built. The canonical
viewpoint within this approach is that space-time singularities are an
artifice of the classical GR description, and that an extension of it via
additional scalar objects, supposedly emerging as low-energy effective (Planck
scale-suppressed) contributions from a quantized version of it, should come to
save the day. In this approach one faces the new opportunities offered by a
larger flexibility in choosing the main geometrical ingredients imbued in the
enlarged theory and, in particular, those associated to the affine connection:
curvature, torsion, and non-metricity. In this chapter we have granted the
mutually desired divorce to the metric and the affine connection, and
considered a theory whose building blocks are made up of these two independent
entities. The field equations of the theories constructed this way are thus
found by independent variation of the action with respect to each of them,
which requires the introduction of workable methods to solve the connection
equations, a problem which is missing in the metric formulation. We succeeded
at this regard by introducing an auxiliary, connection-compatible metric, in
such a way that the field equations can be written as a set of Einstein-like
ones for it. This procedure can be systematized for a large class of models
built upon scalars out of contractions of the metric with the (symmetric part
of the) Ricci tensor, and consequently dubbed them as Ricci-based theories.
Revolving around these methods we can solve such field equations either by
direct attack, or by using the Einstein frame to construct an equivalent GR
action with non-linearly coupled fields, the latter procedure allowing to find
new solutions of any RBG theory under consideration using a seed solution from
the GR side, via purely algebraic transformations. The latter is dubbed as the
mapping method, which should be a basic element in the toolkit of the metric-
affine practitioner.

The spectrum of (spherically symmetric) solutions found via any of the two
procedures above, is ample both in the gravity and matter sectors. For the
former, $f(\mathcal{R})$ gravity, quadratic gravity, and EiBI gravity, are
theories extensively studied in the literature (including their higher-
dimensional and three-dimensional versions), while in the latter we considered
electromagnetic fields (both Maxwell and non-linear), different classes of
anisotropic fluids, and scalar fields. We actually found different classes of
configurations, from black holes with two and a single (degenerated or not)
horizons to naked (horizonless) objects of several types including de Sitter
cores. Besides the modified horizon structure, the new gravitational
corrections have an impact in the innermost structure, where depending on the
gravity/matter combination one may find a bounce in the radial function, which
is the most interesting case for the hopes of achieving singularity-removal.
Indeed, we identified two main mechanisms for the removal of singularities,
namely, either pushing the bouncing region to the future (or past) boundary of
the space-time beyond the reach of any geodesics, or by showing that when the
bounce is reachable by some sets of geodesics, these can be naturally extended
across of it to reach arbitrarily large values of their affine parameter. In
the latter case, one can wonder what is then the effect of curvature
divergences in those cases where they appear at the bounce, since the
extensibility of geodesics is utterly numb to the presence of them. In
addition, the de Sitter cores familiar of the GR case, which are free of
curvature divergences, come out naturally as a nice surprise within this
formalism. Subsequently, we brought to the table four additional criteria as
seekers of potential pathologies: the presence of unbound tidal forces upon
extended observers, the loss of causal contact upon any two components of a
congruence of geodesics, the extension of the geodesic completeness criterion
to observers with (bound) acceleration, and tests with (scalar) waves. We
extensively discussed these criteria for several of the solutions above,
identifying some in which all of them are satisfied, thus representing non-
singular geometries from this point of view. Nonetheless, while these are
necessary conditions for a space-time to be physically reasonable on its full
extent, this does not guarantee that other types of pathological behaviours
may arise beyond the limited imagination of the authors. In absence of robust
theoretical results telling us when a space-time is regular rather than when
it is singular, this is how far we have been able to go.

To conclude, different opinions can be found in the community regarding the
interpretation and the severity of the difficulties (or lack thereof) with the
existence of singularities within GR, but undoubtedly they are a great
reservoir of ideas to challenge (even if timorously) Einstein’s interpretation
of the geometrical ingredients and physical principles underlying the
gravitational interaction. In the end, only experiment will tell if regular
black holes contain any clue to it.

###### Acknowledgements.

DRG is funded by the Atracción de Talento Investigador programme of the
Comunidad de Madrid (Spain) No. 2018-T1/TIC-10431. This work is supported by
the Spanish Grants PID2019-108485GB-I00, and PID2020-116567GB-C21, funded by
MCIN/AEI/10.13039/ 501100011033 (“ERDF A way of making Europe” and “PGC
Generación de Conocimiento”), the project PROMETEO/2020/079 (Generalitat
Valenciana), the project H2020-MSCA-RISE-2017 Grant FunFiCO- 777740, the FCT
projects No. PTDC/FIS-PAR/31938/2017 and PTDC/FIS-OUT/29048/2017, and the
Edital 006/2018 PRONEX (FAPESQ-PB/CNPQ, Brazil, Grant 0015/2019). This article
is based upon work from COST Action CA18108, supported by COST (European
Cooperation in Science and Technology).

## References

  * (1) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Astrophys. J. Lett. 848, L13 (2017).
  * (2) V. I. Afonso, C. Bejarano, J. Beltran Jimenez, G. J. Olmo and E. Orazi, Class. Quant. Grav. 34, 235003 (2017).
  * (3) V. I. Afonso, G. J. Olmo, E. Orazi and D. Rubiera-Garcia, Eur. Phys. J. C 78, 866 (2018).
  * (4) V. I. Afonso, G. J. Olmo, E. Orazi and D. Rubiera-Garcia, Phys. Rev. D 99, 044040 (2019).
  * (5) V. I. Afonso, G. J. Olmo, E. Orazi and D. Rubiera-Garcia, JCAP 12, 044 (2019).
  * (6) V. I. Afonso, G. J. Olmo and D. Rubiera-Garcia, JCAP 08, 031 (2017).
  * (7) V. I. Afonso, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 97, 021503 (2018).
  * (8) R. Aldrovandi and J. G. Pereira, Teleparallel Gravity: An Introduction (Springer, 2013).
  * (9) S. Ansoldi, [arXiv:0802.0330 [gr-qc]].
  * (10) S. Bahamonde and J. Gigante Valcarcel, Eur. Phys. J. C 81, 495 (2021).
  * (11) C. Bambi, A. Cardenas-Avendano, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 93, 064016 (2016).
  * (12) C. Bambi, D. Rubiera-Garcia and Y. Wang, Phys. Rev. D 94, 064002 (2016).
  * (13) M. Banados and P. G. Ferreira, Phys. Rev. Lett. 105 011101 (2010) [erratum: Phys. Rev. Lett. 113 119901 (2014)].
  * (14) M. Banados, C. Teitelboim and J. Zanelli, Phys. Rev. Lett. 69, 1851 (1992).
  * (15) D. Bazeia, L. Losano, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 90, 044011 (2014).
  * (16) D. Bazeia, L. Losano, G. J. Olmo, D. Rubiera-Garcia and A. Sanchez-Puente, Phys. Rev. D 92, 044018 (2015).
  * (17) C. Bejarano, A. Delhom, A. Jiménez-Cano, G. J. Olmo and D. Rubiera-Garcia, Phys. Lett. B 802, 135275 (2020).
  * (18) C. Bejarano, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 95, 064043 (2017).
  * (19) J. Beltrán Jiménez, D. De Andrés and A. Delhom, Class. Quant. Grav. 37, 225013 (2020).
  * (20) J. Beltrán Jiménez and A. Delhom, Eur. Phys. J. C 79, 656 (2019).
  * (21) J. Beltrán Jiménez, L. Heisenberg and T. S. Koivisto, Universe 5, 173 (2019).
  * (22) J. Beltran Jimenez, L. Heisenberg and G. J. Olmo, JCAP 11, 004 (2014).
  * (23) J. Beltran Jimenez, L. Heisenberg, G. J. Olmo and D. Rubiera-Garcia, JCAP 10, 029 (2017).
  * (24) J. Beltran Jimenez, L. Heisenberg, G. J. Olmo and D. Rubiera-Garcia, Phys. Rept. 727, 1 (2018).
  * (25) N. D. Birrell and P. C.W. Davies, Quantum Fields in Curved Space (Cambridge University Press, Cambridge, England, 1982).
  * (26) M. Born and L. Infeld, Proc. R. Soc. London. A 144, 425 (1934).
  * (27) C. Burrage and J. Sakstein, Living Rev. Rel. 21, 1 (2018).
  * (28) R. Carballo-Rubio, F. Di Filippo, S. Liberati and M. Visser, Class. Quant. Grav. 37, 14 (2020).
  * (29) R. Carballo-Rubio, F. Di Filippo, S. Liberati and M. Visser, Phys. Rev. D 101, 084047 (2020).
  * (30) C.J.S. Clarke, A. Królak, Journal of Geometry and Physics, 2, 127 (1985).
  * (31) E. Curiel, Stanford Encyclopedia of Philosophy.
  * (32) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010).
  * (33) J. Diaz-Alonso and D. Rubiera-Garcia, Phys. Rev. D 81, 064021 (2010).
  * (34) I. Dymnikova, Int. J. Mod. Phys. D 12, 1015 (2003).
  * (35) I. Dymnikova, Class. Quant. Grav. 21, 4417 (2004).
  * (36) G. F. R. Ellis and B. G. Schmidt, Gen. Rel. Grav. 8, 951 (1977).
  * (37) E. S. Fradkin and A. A. Tseytlin, Phys. Lett. B 163, 123 (1985).
  * (38) Y. Friedman and T. Scarr, Gen. Rel. Grav. 47, 121 (2015).
  * (39) R. P. Geroch, Annals Phys. 48, 526 (1968).
  * (40) M. Guerrero and D. Rubiera-Garcia, Phys. Rev. D 102, 024005 (2020).
  * (41) M. Guerrero, G. Mora-Pérez, G. J. Olmo, E. Orazi and D. Rubiera-Garcia, JCAP 11, 025 (2021).
  * (42) J.B. Hartle, S.W. Hawking, Commun. Math. Phys. 26, 87 (1972).
  * (43) S. I. Kruglov, Mod. Phys. Lett. A 32, 1750092 (2017).
  * (44) S. Jana, G. K. Chakravarty and S. Mohanty, Phys. Rev. D 97, 084011 (2018).
  * (45) J.P.S. Lemos, V.T. Zanchin, Phys. Rev. D 77, 064003 (2008).
  * (46) J. P. S. Lemos and V. T. Zanchin, Phys. Rev. D 83, 124005 (2011).
  * (47) F. S. N. Lobo, G. J. Olmo and D. Rubiera-Garcia, JCAP 07, 011 (2013).
  * (48) S.D. Majumdar, Phys. Rev. D. 72, 930 (1947).
  * (49) J. W. Maluf, Annalen Phys. 525, 339 (2013).
  * (50) P. Martin-Moruno and M. Visser, Fundam. Theor. Phys. 189, 193 (2017).
  * (51) C. Menchon, G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 96, 104028 (2017).
  * (52) J. M. Nester and H. J. Yo, Chin. J. Phys. 37, 113 (1999).
  * (53) B. C. Nolan, Phys. Rev. D 60, 024014 (1999).
  * (54) Y. N. Obukhov and J. G. Pereira, Phys. Rev. D 67, 044016 (2003).
  * (55) G. J. Olmo, Phys. Rev. Lett. 95, 261102 (2005).
  * (56) G. J. Olmo, Phys. Rev. D 72, 083505 (2005).
  * (57) G. J. Olmo, Int. J. Mod. Phys. D 20, 413 (2011).
  * (58) G. J. Olmo, E. Orazi and D. Rubiera-Garcia, Eur. Phys. J. C 80, 1018 (2020).
  * (59) G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 84, 124059 (2011).
  * (60) G. J. Olmo and D. Rubiera-Garcia, Eur. Phys. J. C 72, 2098 (2012).
  * (61) G. J. Olmo and D. Rubiera-Garcia, Phys. Rev. D 86, 044014 (2012).
  * (62) G. J. Olmo and D. Rubiera-Garcia, Universe 1, 173 (2015).
  * (63) G. J. Olmo, D. Rubiera-Garcia and H. Sanchis-Alepuz, Eur. Phys. J. C 74, 2804 (2014).
  * (64) G. J. Olmo, D. Rubiera-Garcia and A. Sanchez-Puente, Eur. Phys. J. C 76, 143 (2016).
  * (65) G. J. Olmo, D. Rubiera-Garcia and A. Sanchez-Puente, Class. Quant. Grav. 33, 115007 (2016).
  * (66) G. J. Olmo, D. Rubiera-Garcia and A. Sanchez-Puente, Class. Quant. Grav. 35, 055010 (2018).
  * (67) E. Orazi, Int. J. Mod. Phys. D 29, 2041010 (2020).
  * (68) A. Ori, Phys. Rev. D 61, 064016 (2000).
  * (69) A. Papapetrou, Proc. R. Irish Acad. A51, 191 (1947).
  * (70) L. Parker and D. J. Toms, Quantum Field Theory in Curved Spacetime: Quantized Fields and Gravity (Cambridge University Press, Cambridge, England, 2009).
  * (71) R. Penrose, Phys. Rev. Lett. 14, 57 (1965).
  * (72) T. A. Roman, Phys. Rev. D 33, 3526 (1986).
  * (73) J. M. M. Senovilla and D. Garfinkle, Class. Quant. Grav. 32, 124008 (2015).
  * (74) F. J. Tipler, Phys. Rev. D 15, 942 (1977).
  * (75) M. Visser, Lorentzian Wormholes (Springer-Verlag, NY, 1996).
  * (76) D. N. Vollick, Phys. Rev. D 69, 064030 (2004).
  * (77) C. M. Will, Living Rev. Rel. 17, 4 (2014).
  * (78) M. Wyman, Phys. Rev. D 24, 839 (1981).


# Pascal’s Triangle Fractal Symmetries

Nayan E. Myerson-Jain  Department of Physics, University of California, Santa
Barbara, CA 93106, USA  Shang Liu  Kavli Institute for Theoretical Physics,
University of California, Santa Barbara, CA 93106, USA  Wenjie Ji  Department
of Physics, University of California, Santa Barbara, CA 93106, USA  Cenke Xu
Department of Physics, University of California, Santa Barbara, CA 93106, USA
Sagar Vijay  Department of Physics, University of California, Santa Barbara,
CA 93106, USA

###### Abstract

We introduce a model of interacting bosons exhibiting an infinite collection
of fractal symmetries – termed “Pascal’s triangle symmetries” – which provides
a natural ${\rm U}(1)$ generalization of a spin-(1/2) system with Sierpinski
triangle fractal symmetries introduced in Ref. Newman and Moore (1999). The
Pascal’s triangle symmetry gives rise to exact degeneracies, as well as a
manifold of low-energy states which are absent in the Sierpinski triangle
model. Breaking the ${\rm U}(1)$ symmetry of this model to $Z_{p}$, with prime
integer $p$, yields a lattice model with a unique fractal symmetry which is
generated by an operator supported on a fractal subsystem with Hausdorff
dimension $d_{H}=\ln(p(p+1)/2)/\ln p$. The Hausdorff dimension of the fractal
can be probed through correlation functions at finite temperature. The phase
diagram of these models at zero temperature in the presence of quantum
fluctuations, as well as the potential physical construction of the ${\rm
U}(1)$ model are discussed.

Introduction: In recent years, generalizations of the notion of symmetry have
significantly broadened our view on states of matter. Many previously
developed phases, such as $Z_{2}$ topological order Read and Sachdev (1991);
Wen (1991); Moessner and Sondhi (2001); Kitaev (2003), and the $(3+1)d$
algebraic spin liquid with photon excitations Wen (2003); Moessner and Sondhi
(2003); Hermele et al. (2004), were thought to be beyond the notion of
spontaneous symmetry breaking (SSB) of ordinary global symmetries. But it has
been realized in recent years that these phases still have a unified
description as the SSB of generalized higher-form symmetries 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); Córdova et al. (2019). The notion of subsystem symmetry has
further enriched our understanding along this line. Long-range-entangled
quantum phases with fractionalized excitations that exhibit restricted
mobility – termed fracton orders – exhibit emergent subsystem symmetries Vijay
et al. (2016, 2015); Haah (2011); Nandkishore and Hermele (2019); Pretko et
al. (2020). Subsystem symmetries can be categorized into type-I and II Vijay
et al. (2016), where a type-I subsystem symmetry has generators and conserved
charges defined on regular submanifolds of the system, such as lines and
planes, while type-II subsystem symmetries have conserved charges defined on a
fractal-shaped subsystem Yoshida (2013); Haah (2011), often with non-integer
spatial dimensions.

The simplest model with a fractal subsystem symmetry is the Sierpinski-
triangle model, which was first introduced for the purpose of studying glassy
dynamics in the absence of disorder Newman and Moore (1999):

$H_{\mathrm{ST}}=\sum_{\bigtriangledown}-K\sigma^{z}_{1}\sigma^{z}_{2}\sigma^{z}_{3}.$
(1)

Here $\sigma^{z}_{i}$ is an Ising spin defined on each site of a triangular
lattice, and the sum is only over the downward-facing triangular plaquettes of
the lattice. This simple model has the following desirable features. (1) The
model has an exotic “fractal symmetry”, which becomes most explicit when the
system is defined on a $L\times L$ lattice with $L=2^{k}-1$: the Hamiltonian
is invariant under flipping spins on a Sierpinski-triangle fractal subsystem;
(2) at finite temperature, the three point correlation function
$\langle\sigma^{z}_{0,0}\sigma^{z}_{r,0}\sigma^{z}_{0,r}\rangle$ is nonzero
only when $r=2^{k}$, and scales as $\sim\exp(-r^{d_{H}})$ with $d_{H}=\ln
3/\ln 2$, which is the Hausdorff dimension of the Sierpinski-triangle Yoshida
(2013); (3) with a transverse field $\sum_{i}-h\sigma^{x}_{i}$, the system
becomes a quantum Sierpinski-triangle model, and at $h=K$ there is a quantum
phase transition at zero temperature Vasiloiu et al. (2020); Zhou et al.
(2021), which separates the “fractal-ordered” phase that spontaneously breaks
the fractal symmetry ($K>h$), and a disordered phase with ($h>K$).

Figure 1: Pascal’s Triangle Symmetries: The ${\rm U}(1)$ parent model (2) has
a family of fractal symmetries, generated by a staggered rotation of the boson
phase $\theta_{i}\rightarrow\theta_{i}+\frac{2\pi}{p}m_{i}$ over a triangular
region of side length which is a power of any prime number $p$. For each $p$,
the fractal symmetry becomes exact for system size $L^{2}$ with $L=p^{k}-1$.
When acting on a classical ground-state of the parent model, this
transformation generates excitations at the corners of the triangular region.
The non-trivial action of this rotation can be visualized as Pascal’s triangle
modulo $p$, which is a fractal with Hausdorff dimension
$d_{H}(p)={\ln(p(p+1)/2)}/{\ln p}$.

In this work we introduce generalizations of both the classical and quantum
Sierpinski-triangle models. These models are obtained from a ${\rm U}(1)$
parent model with “Pascal’s triangle” (also called Yang Hui triangle in China)
symmetries, a family of symmetry transformations along a fractal region which
are exact in a system with periodic boundary conditions, and for particular
system sizes. Even when these symmetries are not exact, the presence of an
“approximate” Pascal’s triangle symmetry gives rise to low-energy states which
are absent in the Sierpinski triangle model Newman and Moore (1999).
Descendant models are obtained by reducing the ${\rm U}(1)$ degree of freedom
of the parent model to $Z_{p}$, with prime integer $p$. The $Z_{p}$ models
have their own fractal symmetry that are deduced from the Pascal’s triangle
symmetry, and their degenerate excitations have an emergent fractal structure
with Hausdorff dimension $d_{H}=\ln(p(p+1)/2)/\ln p$.

The ${\rm U}(1)$ parent model: The Hamiltonian of the ${\rm U}(1)$
generalization of the Sierpinski-triangle model reads

$\displaystyle H_{{\rm
U}(1)}=\sum_{\bigtriangledown}-t\cos\left(\theta_{1}+\theta_{2}+\theta_{3}\right).$
(2)

It is straightforward to see that Eq. 2 has two conventional ${\rm U}(1)$
global symmetries:

$\displaystyle{\rm U}(1)_{1}:\theta_{j\in A}\rightarrow\theta_{j\in
A}+\alpha,\ \ \theta_{j\in B}\rightarrow\theta_{j\in B}-\alpha,$ (3) (4)
$\displaystyle{\rm U}(1)_{2}:\theta_{j\in A}\rightarrow\theta_{j\in A}+\beta,\
\ \theta_{j\in C}\rightarrow\theta_{j\in C}-\beta.$ (5)

$A,B,C$ are the three sublattices of the triangular lattice. The ground states
of Eq. 2 spontaneously break the two ${\rm U}(1)$ symmetries. Starting with
one of the ground states, say $\theta=0$ uniformly on the entire lattice, a
class of ground states can be generated by rotating $\theta$ globally
according to Eq. 5. Any ground state obtained this way still has a uniform
order of $\theta$ on each of the three sublattices, hence the ground states
generated through Eq. 5 have a conventional “$\sqrt{3}\times\sqrt{3}$” order,
which is the order often observed on the triangular lattice antiferromagnet.

Besides the two ordinary ${\rm U}(1)$ global symmetries, this model Eq. 2
actually contains an infinite series of $Z_{p}$ distinct fractal symmetries,
one for each prime number $p$, $p\geq 2$. In contrast to the classical
Sierpinski-triangle model, where the $Z_{2}$ fractal subsystem symmetry is
taken as a sequence of spin flips in the shape of a Sierpinski-triangle, the
series of $Z_{p}$ fractal symmetries exhibited by the ${\rm U}(1)$ parent
model are in the shape of a Pascal’s triangle modulo $p$. For example when
$p=2$, this Pascal’s triangle symmetry reduces down to the familiar $Z_{2}$
fractal symmetry of the Sierpinski-triangle model; for $p=3$ the Pascal’s
triangle modulo 3 reduces to another fractal shape (Fig. 1).

The exact series of fractal transformations of Eq. 2 can be written down as a
staggered rotation of the $\theta_{i}$’s in the shape of a Pascal’s triangle
modulo $p$, which has a side-length of $p^{k}-1$, where $k$ is any integer
greater than zero. The precise form of the transformation is

$\displaystyle\theta_{i}\longrightarrow\theta_{i}+\frac{2\pi}{p}(-1)^{i_{x}+i_{y}}{i_{x}+i_{y}\choose
i_{y}}.$ (6)

at the points $(i_{x},i_{y})$ for which $0\leq i_{y}\leq i_{x}$ and
$i_{x}+i_{y}\in[0,p^{k}-1]$. As shown in the Supplemental Material (SM),
transformations of this form can be used to generate exact symmetries when the
system is placed on an $L\times L$ lattice with periodic boundary conditions
and with $L=p^{k}-1$.

Any $Z_{p}$ fractal transformation of the ${\rm U}(1)$ parent model which is
not an exact symmetry of the model generates fully immobile defects which are
analogous to fractons. From the uniform $\theta_{i}=0$ ground-state,
transforming the ${\rm U}(1)$ degrees of freedom according to Eq. 6 in the
shape of a local Pascal’s triangle of size $p^{k}-1$ creates three defects of
energy $t(1-\cos({2\pi}/{p}))$, one at each downward-facing triangular
plaquette located at the corners of the Pascal’s triangle, as shown in the SM,
and as indicated schematically in Fig. 1. If we treat these defects as point-
like quasiparticles localized on their downward facing plaquettes, individual
defects cannot be moved by any rotation of $\theta_{i}$’s without creating
more excitations and are hence completely immobile. We note that since $p$ can
be arbitrarily large, the excitations created by a $Z_{p}$ fractal
transformation can cost a vanishingly small energy in a thermodynamically
large system.

At finite temperature, the ${\rm U}(1)$ parent model is completely disordered
similarly to the Sierpinski-triangle model Newman and Moore (1999). This can
be most easily seen from a duality mapping of the ${\rm U}(1)$ degrees of
freedom on the vertices to new ${\rm U}(1)$ degrees of freedom on downward
facing plaquettes
$(\theta_{1}+\theta_{2}+\theta_{3})_{\bigtriangledown}\rightarrow\phi_{\bigtriangledown}$,
where $\phi_{\bigtriangledown}$ is defined on the dual site located at the
center of each downward facing triangular plaquette (Fig. 2), and $\phi$ is
still compact (periodically defined). The dual of the ${\rm U}(1)$ parent
model is

$\displaystyle H_{{\rm
U}(1)}^{d}=\sum_{\bigtriangledown}-t\cos(\phi_{\bigtriangledown}).$ (7)

Since each $\phi$ is decoupled from one another, the partition function
factorizes into a product of local partition functions for each individual
$\phi$ which does not support any phase transition.

The three-body interactions of the Sierpinski-triangle model, as well as the
${\rm U}(1)$ parent model look artificial. Ref. Myerson-Jain et al. (2021)
proposed to realize the Sierpinski-triangle model with the Rydberg atoms with
only two-body Van der Waals interactions. In the SM we present a more natural
construction of the ${\rm U}(1)$ parent model through a set-up with only two-
body interactions.

The $Z_{p}$ models: From the ${\rm U}(1)$ parent model Eq. 2, models with a
single fractal symmetry that are natural extensions of the Sierpinski-triangle
model can be constructed. This is done by breaking the ${\rm U}(1)$ degrees of
freedom down to $Z_{p}$ clock degrees of freedom
$\sigma_{i}=e^{i\theta_{i}},\theta_{i}\in\frac{2\pi}{p}Z_{p}$:

$\displaystyle H_{Z_{p}}$ $\displaystyle=$
$\displaystyle\sum_{\bigtriangledown}-\frac{t}{2}\sigma_{1}\sigma_{2}\sigma_{3}+\text{h.c.}$
(8)

The model Eq. 8 extends many properties of the Sierpinski-triangle model to a
more general series of $Z_{p}$ “Pascal’s triangle models” which reduces to Eq.
1 when $p=2$. Generally, since the $Z_{p}$ fractal symmetry of the Pascal’s
triangle models are descended from the ${\rm U}(1)$ parent model, the fractal
symmetry transformation in these models is realized by Eq. 6 with the
appropriate choice of $p$. These series of models also display the fracton-
like defects associated to fractal excitations in the shape of a Pascal’s
triangle modulo $p$ that cost energy $t(1-\cos(\frac{2\pi}{p}))$ each as well
as spontaneously breaking the $Z_{p}$ fractal symmetry, yielding a ground-
state degeneracy of $p^{L-1}$ when $L=p^{k}-1$ (see SM for derivation).

Spontaneous symmetry breaking of the $Z_{p}$ fractal symmetries in the
Pascal’s triangle models can be diagnosed by the behavior of a three-point
correlation function. Making use of the duality of these models, we can define
plaquette degrees of freedom
$\tau_{\bigtriangledown}=(\sigma_{1}\sigma_{2}\sigma_{3})_{\bigtriangledown}$.
The dual Hamiltonian is

$\displaystyle
H_{Z_{p}}^{d}=\sum_{\bigtriangledown}-\frac{t}{2}\tau_{\bigtriangledown}+\text{h.c.}$
(9)

In the thermodynamic limit, each $\sigma$-variable can be represented as an
infinite staggered product of dual $\tau$-variables in the shape of a Pascal’s
triangle modulo $p$. The three-point function
$\mathcal{C}_{3}(r)=\langle\sigma_{0,0}\sigma_{r,0}\sigma_{0,r}\rangle$, after
being rewritten in terms of the dual variables only has compact support when
$r=p^{k}$ and hence must vanish else wise. The three-point function factors
into a product of single-site expectation values
$\langle\tau\rangle,\dots,\langle\tau^{p-1}\rangle$. From the form of the
Hamiltonian, a general expression for the three-point function for arbitrary
$p$ prime can be derived (see SM for details)

$\displaystyle\mathcal{C}_{3}(r=p^{k})=\prod_{m=1}^{\frac{p-1}{2}}\langle\tau^{m}\rangle^{N_{m,p-m}(k)},$
(10)

$N_{m,p-m}(k)$ is the number of times $m$ and $p-m$ appear in a Pascal’s
triangle modulo $p$ with length $p^{k}-1$. Such an expression is complex, but
a complete set of recurrence relations is constructed in the SM for
$N_{m,p-m}(k)$ lending Eq. 10 to efficient numerical evaluation. For small $p$
this can be done analytically, e.g. for $p=3$, the three-point function is

$\displaystyle\mathcal{C}_{3}(r=3^{k})=\langle\tau\rangle^{r^{d_{H}}}$ (11)
$\displaystyle=$ $\displaystyle\bigg{(}\frac{e^{\beta t}-e^{-\beta
t/2}}{e^{\beta t}+2e^{-\beta t/2}}\bigg{)}^{r^{d_{H}}}\sim e^{-\alpha
r^{d_{H}}},$ (13)

where $d_{H}=\frac{\ln 6}{\ln 3}$ is the Hausdorff dimension of a Pascal’s
triangle modulo $3$. In general there are natural bounds on the decay of the
three-point function: $\mathcal{C}_{3}(r=p^{k})$ will decay hyper-
exponentially as $\sim e^{-\alpha r^{\gamma}}$ for sufficiently large $r$,
where the exponent of the decay $\gamma$ is bounded such that

$\displaystyle d_{H}\leq\gamma\leq\bigg{(}\frac{p-1}{2}\bigg{)}d_{H},\text{
}p>2.$ (14)

As demonstrated by Eqs. 10 and 14, the decay of the three-point function can
be complicated for general $p$ prime. However, a modified version of the
$Z_{p}$ Pascal’s triangle models in Eq. 8 can be proposed for which the three-
point function at finite temperature always decays as a fractal area-law. If
we consider an equal-weight summation of plaquette terms

$\displaystyle\mathcal{H}_{\text{p}}=\sum_{\bigtriangledown}\sum_{m=0}^{\frac{p-1}{2}}-\frac{t}{2}(\sigma_{1}\sigma_{2}\sigma_{3})^{m}+\text{h.c.}$
(15)

This model retains the $Z_{p}$ fractal symmetry of Eq. 8 as it only includes
products of the original plaquette terms. As such, the duality
$(\sigma_{1}\sigma_{2}\sigma_{3})_{\bigtriangledown}\rightarrow\tau_{\bigtriangledown}$
still exists and the dual of $\mathcal{H}_{p}$ is

$\displaystyle\mathcal{H}_{\text{p}}^{d}=\sum_{\bigtriangledown}\sum_{m=0}^{\frac{p-1}{2}}-\frac{t}{2}\tau_{\bigtriangledown}^{m}+\text{h.c.}$
(16)

The manner in which the three-point correlation for $\mathcal{H}_{\text{p}}$
is calculated remains the same as for what it was in Eq. 8 with the exception
that $\langle\tau^{m}\rangle$ no longer depends on power, $m$. As a result,
$\mathcal{C}_{3}(r=p^{k})$ decays as a fractal area-law no matter what value
$p$ takes:

$\displaystyle\mathcal{C}_{3}(r=p^{k})=\langle\tau\rangle^{r^{d_{H}}},$ (17)
$\displaystyle=$ $\displaystyle\bigg{(}\frac{e^{(\frac{p+1}{2})\beta
t}-e^{\frac{\beta t}{2}}}{e^{(\frac{p+1}{2})\beta t}+(p-1)e^{\frac{\beta
t}{2}}}\bigg{)}^{r^{d_{H}}}\sim e^{-\alpha r^{d_{H}}}.$ (19)

The $Z_{p}$ Pascal’s triangle model with prime integer $p$ can be further
extended to $Z_{N}$ models where $N$ is a composite positive integer. These
new composite $Z_{N}$ models have more than one fractal symmetry. In fact,
there is a distinct $Z_{p}$ fractal symmetry for each unique prime divisor of
$N$, e.g. the $Z_{N=6}$ model has a $Z_{2}$ Sierpinski-triangle fractal
symmetry and a $Z_{3}$ Pascal’s triangle modulo $3$ fractal symmetry. The
behavior of the three-point correlations of the $Z_{N}$ models is further
discussed in the SM.

Figure 2: The U(1) Model and its Dual Description: The models we consider in
this work involve sum of all the downward facing triangles (shaded in green);
The dual of the ${\rm U}(1)$ model (Eq. 7, Eq. 23) is defined on the dual
triangular lattice, whose sites are the center of each downward facing
triangles of the original lattice.

The Quantum Phase diagram: So far we have only discussed the classical version
of the models. To turn quantum fluctuations on in Eq. 2, one can modify the
model as

$\displaystyle H_{Q-{\rm
U}(1)}=\sum_{\bigtriangledown}-t\cos\left(\theta_{1}+\theta_{2}+\theta_{3}\right)+\sum_{j}\frac{U}{2}n_{j}^{2},$
(20)

where $n_{j}$ is the boson number operator defined on each site of the
triangular lattice, which are conjugate to the boson phase
$[n_{i},\theta_{j}]=\mathrm{i}\delta_{ij}$.

As shown previously, a class of ground states of the classical ${\rm U}(1)$
model Eq. 2 have the conventional $\sqrt{3}\times\sqrt{3}$ order, which
spontaneously breaks the Pascal’s triangle symmetry, and the two ${\rm U}(1)$
symmetries in Eq. 5. We now investigate whether this classical order is stable
against quantum fluctuation, i.e. whether it is stable against the $U$ term in
Eq. 20. We argue that all symmetries of the Hamiltonian (20) are restored by
quantum fluctuations.

In an ordered phase that spontaneously breaks the ${\rm U}(1)$ global
symmetry, one can ignore the fact that the phase angle $\theta$ is a compact
boson (i.e. $\theta\sim\theta+2\pi$), and hence expand the cosine function of
Eq. 20 to the lowest nontrivial order. This procedure leads to an approximate
Gaussian Hamiltonian:

$\displaystyle H^{g}_{Q-{\rm
U}(1)}=\sum_{i,j}t\theta_{i}\theta_{j}+\sum_{j}3t\theta_{j}^{2}+\frac{U}{2}n_{j}^{2}.$
(21)

The band structure of $\theta$ based on this Gaussian Hamiltonian has minima
at $\pm\mathbf{K}=\pm(4\pi/3,0)$, which is consistent with the
$\sqrt{3}\times\sqrt{3}$ order of the classical Hamiltonian. The spectrum of
the Gaussian Hamiltonian is gapless.

The Gaussian expansion of the Hamiltonian ignored the compactness of $\theta$.
In a quantum model constructed with $\theta$, $\theta$ being a compact boson
is equivalent to the constraint that its quantum conjugate variable $n$ take
discrete values. To check the stability of a semiclassical state of $\theta$
under quantum fluctuation, one needs to investigate whether the compactness of
$\theta$, or equivalently the discrete nature of $n$ would destabilize the
semiclassical state described by the Gaussian Hamiltonian Eq. 21. For example,
the $(2+1)d$ quantum dimer model on the square lattice can be mapped to a
compact ${\rm U}(1)$ gauge theory Rokhsar and Kivelson (1988); FRADKIN and
KIVELSON (1990), a Gaussian expansion would lead to gapless photons. But the
compactness of the gauge field is always relevant in a semiclassical photon
state unless the system is at a fine-tuned multicritical point (the so-called
RK point Rokhsar and Kivelson (1988)), hence the Gaussian state is generally
unstable against quantum fluctuations. This effect is also referred to as
confinement of a lattice gauge theory. The standard method of the analysis
relies on the dual formalism of Eq. 20 and Eq. 21. The dual model is defined
on the dual triangular lattice (Fig. 2), by introducing the following
variables:

$\displaystyle\sum_{j\in\bigtriangledown}\theta_{j}=\phi_{\bar{j}},\ \ \ \
n_{j}=\sum_{\bar{j}\in\triangle\ \mathrm{around}\ j}\Psi_{\bar{j}}.$ (22)

where $\bar{j}$ denotes the sites of the dual triangular lattice.
$\phi_{\bar{j}}$ and $\Psi_{\bar{j}}$ are canonically conjugate variables.
$\Psi_{\bar{j}}$ takes half-integer values, and $\phi_{\bar{j}}$ is compact.
The dual Hamiltonian reads

$\displaystyle H_{Q-{\rm
U}(1)}^{d}=\sum_{\bar{j}}-t\cos(\phi_{\bar{j}})+\sum_{\bar{\triangle}}\frac{U}{2}(\Psi_{\bar{1}}+\Psi_{\bar{2}}+\Psi_{\bar{3}})^{2}.$
(23)

Instead of directly dealing with the discrete variable $\Psi$, we may view
$\Psi_{\bar{j}}$ as taking continuous values, and $\phi_{\bar{j}}$ as its non-
compact conjugate variable. The discrete nature of $\Psi$ can be enforced
through an external potential in the dual Hamiltonian. The dual Hamiltonian
becomes

$\displaystyle H_{Q-{\rm U}(1)}^{d}$ $\displaystyle\sim$
$\displaystyle\sum_{\bar{\triangle}}\frac{U}{2}(\Psi_{\bar{1}}+\Psi_{\bar{2}}+\Psi_{\bar{3}})^{2}-\sum_{\bar{j}}t\cos(\phi_{\bar{j}})$
(24) $\displaystyle-$ $\displaystyle\alpha\cos(2\pi\Psi_{\bar{j}}).$ (26)

The next step is to temporarily ignore the $\alpha$-terms, and expand
$-t\cos(\phi_{j})$ to the lowest nontrivial order. After this procedure, the
dual Hamiltonian takes a Gaussian form, and it is the dual of the Gaussian
Hamiltonian Eq. 21. The goal of this analysis is to check the role of the
$\alpha$-terms at this Gaussian state. This dual Gaussian Hamiltonian can be
solved, leading to a band structure of $\Psi$. The minima of the band struture
of $\Psi$ are located at the two corners of the Brillouin zone,
$\pm\mathbf{K}=(\pm 4\pi/3,0)$. We then expand $\Psi_{\mathbf{r}}$ at
$\pm\mathbf{K}$:

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

The Lagrangian of the dual theory expanded at $\pm\mathbf{K}$ becomes

$\displaystyle\mathcal{L}_{Q-{\rm
U}(1)}^{d}=(\partial_{\tau}\vec{\psi})^{2}+\rho_{2}(\nabla\vec{\psi})^{2}-\sum_{a}\alpha\cos(\vec{e}_{a}\cdot\vec{\psi}),$
(28)

where $\vec{\psi}=(\mathrm{Re}(\psi),\mathrm{Im}(\psi))$; $a=A,B,C$ label the
three sublattics of the dual triangular lattice, and $e_{A}=2\pi(1,0)$,
$e_{B}=2\pi(-1/2,\sqrt{3}/2)$, $e_{C}=2\pi(-1/2,-\sqrt{3}/2)$.

The last three terms in Eq. 28 arise from rewriting the last term of Eq. 26 by
expanding $\Psi$ at $\pm\mathbf{K}$. After this expansion, the last term of
Eq. 26 becomes $-\alpha\cos(\vec{e}_{a}\cdot\vec{\psi}(\mathbf{r}))$ for
$\mathbf{r}$ belonging to sublattice $a$ ($a=A,B,C$) of the dual triangular
lattice. Hence at long scale a nonvanishing term would survive. The $\alpha$
term in Eq. 28 will be relevant for the Gaussian theory with nonzero
$\rho_{2}$, which implies that the compactness of $\theta$, or the discrete
nature of $n$ in Eq. 20 destabilizes the semiclassical Gaussian state, and the
spectrum of Eq. 20 should be gapped even with small $U$.

The dual description of the U(1) model studied here accurately captures the
spectrum of the gapless modes arising from spontaneously breaking the global
U(1) symmetries, though it does not reproduce the spectrum at $U=0$ that arise
due to the Pascal triangle symmetries. Nevertheless, the nature of the ground-
state of the system when the pinning potential flows to strong coupling can
still be inferred. A strong $\alpha$ would pin $\Psi$ to integer values, which
implies that a relevant $\alpha$ would drive the system into an eigenstate of
$n$ in Eq. 20, and the $U$ term will lead to a unique and gapped ground state
without any spontaneous symmetry breaking. Hence we postulate that quantum
fluctuations of Eq. 20 restores all the symmetries of model Eq. 2, and
continuously connects to the large$-U$ limit of Eq. 20. The analysis here
would be more involved if $n$ takes half integer values in Eq. 20.

One possible quantum generalization of Eq. 8 is

$\displaystyle
H_{Q-Z_{p}}=\sum_{\bigtriangledown}-t\sigma_{1}^{z}\sigma_{2}^{z}\sigma_{3}^{z}-\sum_{j}h\sigma^{x}_{j}+\text{h.c.}$
(29)

for which the clock operators $\sigma^{z}$ and $\sigma^{x}$ obey
$(\sigma^{z})^{p}=(\sigma^{x})^{p}=1$ and $\sigma^{z}\sigma^{x}=e^{2\pi
i/p}\sigma^{x}\sigma^{z}$ (one can also take
$\sigma^{z}=\exp(\mathrm{i}\theta)$ and $\sigma^{x}=\exp(\mathrm{i}2\pi n/p)$,
and restrict $\theta$ to take values in $\frac{2\pi}{p}Z_{p}$). Unlike the
quantum ${\rm U}(1)$ model, Eq. 29 is exactly self-dual with the introduction
of the dual plaquette variables

$\displaystyle\tau^{x}_{\bar{j}}=\sigma_{1}^{z}\sigma_{2}^{z}\sigma_{3}^{z},\
\ \ \ \tau^{z}_{\bar{1}}\tau^{z}_{\bar{2}}\tau_{\bar{3}}^{z}=\sigma^{x}_{j}$
(30)

for which the dual Hamiltonian takes the same form as Eq. 29 with $t$ and $h$
switched. Since the spectrum of the $Z_{p}$ models at $h=0$ is gapped, and it
takes infinite order of perturbations of $h$ to mix two different ground
states in the thermodynamics limit, the classical fractal-order of the $Z_{p}$
Pascal’s triangle models is not destroyed upon the introduction of quantum
fluctuations. Furthermore, the exact self-duality implies that there should be
one or more quantum phase transitions that separate the fractal ordered phase
$(t\gg h)$ and the disordered phase $(h\gg t)$.

Discussion: Although we demonstrated that the semiclassical order of Eq. 20 is
unstable against quantum fluctuation, some deformation of Eq. 20 can support a
stable semiclassical order. In the SM we will show that if we sum over three-
boson interactions for both upward-facing and downward-facing triangles, the
semiclassical $\sqrt{3}\times\sqrt{3}$ order becomes stable against quantum
fluctuations. Also, the system may be tuned to a multicritical point where
$\rho_{2}$ in Eq. 28 vanishes, and the low energy dynamics is controlled by
$\rho_{4}(\nabla^{2}\vec{\psi})^{2}$. The system can remain gapless for a
finite range of $\rho_{4}$, though it takes tuning multiple parameters to
reach this state Fradkin et al. (2004); Ardonne et al. (2004); Vishwanath et
al. (2004).

The nature of the quantum phase transition(s) in the quantum $Z_{p}$ model is
a challenging subject. So far there is no well established paradigm for
understanding quantum phase transitions involving spontaneous breaking of a
fractal symmetry. Any approach to study the quantum phase transition of the
$Z_{p}$ models (such as the quantum $Z_{2}$ Sierpinski-triangle model) through
the ${\rm U}(1)$ generalization would need to address the enlarged Pascal’s
triangle symmetry pointed out in the current work.

C.X. is supported by NSF Grant No. DMR-1920434, and the Simons Investigator
program.

## References

  * 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.
  * Read and Sachdev (1991) N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991), URL https://link.aps.org/doi/10.1103/PhysRevLett.66.1773.
  * Wen (1991) X. G. Wen, Phys. Rev. B 44, 2664 (1991), URL https://link.aps.org/doi/10.1103/PhysRevB.44.2664.
  * Moessner and Sondhi (2001) R. Moessner and S. L. Sondhi, Phys. Rev. Lett. 86, 1881 (2001), URL https://link.aps.org/doi/10.1103/PhysRevLett.86.1881.
  * Kitaev (2003) A. Kitaev, Annals of Physics 303, 2 (2003), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/S0003491602%000180.
  * Wen (2003) X.-G. Wen, Phys. Rev. B 68, 115413 (2003), URL https://link.aps.org/doi/10.1103/PhysRevB.68.115413.
  * Moessner and Sondhi (2003) R. Moessner and S. L. Sondhi, Phys. Rev. B 68, 184512 (2003), URL https://link.aps.org/doi/10.1103/PhysRevB.68.184512.
  * Hermele et al. (2004) M. Hermele, M. P. A. Fisher, and L. Balents, Phys. Rev. B 69, 064404 (2004), URL https://link.aps.org/doi/10.1103/PhysRevB.69.064404.
  * 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.
  * Córdova et al. (2019) C. Córdova, T. T. Dumitrescu, and K. Intriligator, Journal of High Energy Physics 2019, 184 (2019), eprint 1802.04790.
  * 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.
  * Vijay et al. (2015) S. Vijay, J. Haah, and L. Fu, Phys. Rev. B 92, 235136 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.92.235136.
  * Haah (2011) J. Haah, Phys. Rev. A 83, 042330 (2011), URL https://link.aps.org/doi/10.1103/PhysRevA.83.042330.
  * 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.
  * Yoshida (2013) B. Yoshida, Phys. Rev. B 88, 125122 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.125122.
  * 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.
  * Zhou et al. (2021) Z. Zhou, X.-F. Zhang, F. Pollmann, and Y. You (2021), eprint 2105.05851.
  * Myerson-Jain et al. (2021) N. E. Myerson-Jain, S. Yan, D. Weld, and C. Xu (2021), eprint 2108.07765.
  * 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/S0003491604%%****␣main.bbl␣Line␣250␣****000247.
  * 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.
  * Xu and Vishwanath (2006) C. Xu and A. Vishwanath (2006), eprint unpublished.

##  Appendix A Construction of the ${\rm U}(1)$ parent model

The three-body interaction of the ${\rm U}(1)$ parent model is artificial
because the nature is dominated by two-body interactions. In what follows, we
will discuss a more natural construction of the ${\rm U}(1)$ model through a
set-up with only two-body interactions. The microscopic system we start with
is a honeycomb lattice, with a spin-3/2 degree of freedom $\vec{S}_{R}$ on
each site $R$ of sublattice $\mathcal{A}$ of the honeycomb lattice, and a
spin-1/2 degree of freedom $\vec{s}$ on each site of sublattice $\mathcal{B}$.
The initial model only has two body interactions:

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\sum_{R\in\mathcal{A}}\sum_{b=1}^{3}\
J\left(S^{x}_{R}s^{x}_{R,b}+S^{y}_{R}s^{y}_{R,b}\right)-D(S^{z}_{R})^{2}$ (31)
$\displaystyle=$ $\displaystyle\sum_{R\in\mathcal{A}}\sum_{b=1}^{3}\
\frac{J}{2}\left(S^{+}_{R}s^{-}_{R,b}+S^{-}_{R}s^{+}_{R,b}\right)-D(S^{z}_{R})^{2}$
(33)

The first sum of $R$ is over all the sites on sublattice $\mathcal{A}$, the
second sum $\sum_{b=1}^{3}$ is over the three neighboring sites on sublattice
$\mathcal{B}$ around $R$. We also consider other interactions between second-
neighbor sites, or equivalently first-neighbor sites within respective
sublattices $\mathcal{A}$ and $\mathcal{B}$:

$\displaystyle
H^{\prime}=\sum_{R,R^{\prime}\in\mathcal{A}}J_{1}(S^{x}_{R}S^{x}_{R^{\prime}}+S^{y}_{R}S^{y}_{R^{\prime}})+\sum_{i,j\in\mathcal{B}}J_{2}(s^{x}_{i}s^{x}_{j}+s^{y}_{i}s^{y}_{j}).$
(34)

Figure 3: The honeycomb lattice where Eq. 33 is defined on. The sublattice
$\mathcal{B}$ of the honeycomb lattice coincides with the sites of the
triangular lattice of Eq. 1 and Eq. 2.

The $D$ term in Eq. 33 is an onsite spin anisotropy term. If $D$ is positive
and large, at low energy only the process that tunnels between $S^{z}=-3/2$
and $S^{z}=+3/2$ is allowed. At the third order perturbation of $J$, a
nontrivial term is generated within the low energy Hilbert space, which acts
as a low energy effective Hamiltonian:

$\displaystyle H_{\mathrm{eff},0}$ $\displaystyle=$
$\displaystyle\sum_{R\in\mathcal{A}}K\left((S^{+}_{R})^{3}s^{-}_{R,1}s^{-}_{R,2}s^{-}_{R,3}+(S^{-}_{R})^{3}s^{+}_{R,1}s^{+}_{R,2}s^{+}_{R,3}\right),$
(35) $\displaystyle K$ $\displaystyle\sim$ $\displaystyle\frac{J^{3}}{D^{2}}.$
(37)

Notice that at low energy, only terms with $(S^{+})^{3}$ and $(S^{-})^{3}$ can
survive, while terms like $(S^{+}_{R})^{3}(s^{-}_{R,1})^{2}s^{-}_{R,2}$ must
vanish because $(s^{-})^{2}=0$ for spin-1/2 operator $s^{\pm}$.

Two other effects need to be discussed: (1) at the second order perturbation
of $J/D$, the $J_{2}$ term in Eq. 34 will be renormalized, and shifted by
$\sim J^{2}/D$. We assume that the renormalized $J_{2}$ is very small and
negligible. (2) The 3rd order perturbation of $J_{1}$ will generate a new term

$\displaystyle
H_{\mathrm{eff}}^{\prime}\sim\sum_{R,R^{\prime}\in\mathcal{A}}\frac{J_{1}^{3}}{D^{2}}\left((S_{R}^{+})^{3}(S_{R^{\prime}}^{-})^{3}+\text{h.c.}\right)$
(38)

Now, if we combine Eq. 37 and Eq. 38 together and use the standard spin-boson
mapping: $S^{+}_{R}\sim\exp(\mathrm{i}\theta_{R})$ and
$s^{+}_{R,b}\sim\exp(\mathrm{i}\theta_{R,b})$, the entire low energy effective
Hamiltonian is mapped to

$\displaystyle H_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\sum_{R\in\mathcal{A}}-t\cos\left(3\theta_{R}-\theta_{R,1}-\theta_{R,2}-\theta_{R,3}\right)$
(39) $\displaystyle-$
$\displaystyle\sum_{R,R^{\prime}\in\mathcal{A}}t^{\prime}\cos(3\theta_{R}-3\theta_{R^{\prime}}).$
(41)

Due to the $t^{\prime}$ term, at zero temperature $3\theta_{R}$ will order on
$\mathcal{A}$ sublattice. If we replace $3\theta_{R}$ with its expectation
value, then Eq. 41 reduces to Eq. 2 on a triangular lattice, which corresponds
to the sublattice $\mathcal{B}$ of the original honeycomb lattice.

When we map the spin model to a boson rotor model, there is one extra subtlety
that the boson number takes half-integer, rather than integer values. This to
some extent complicates the analysis of the phase where $\theta$ is
disordered, due to the extra degeneracy arising from $S^{z}$.

##  Appendix B Pascal’s Triangle Symmetry of the ${\rm U}(1)$ Parent Model

In this section, we establish the presence of a generalized fractal symmetry
of the ${\rm U}(1)$ parent Hamiltonian (2) given an appropriate set of
boundary conditions. When these boundary conditions are not satisfied, we
establish that these approximate symmetries give rise to low-energy modes.
Furthermore, we show that the $Z_{p}$ generalization of the Sierpinski
triangle model exhibits a ground-state degeneracy of $D=p^{L-1}$ on an
$L\times L$ system with periodic boundary conditions, when $L=p^{k}-1$.

Figure 4: No excitations are created at the upward-facing triangles that
border an edge of region $\mathcal{R}$, as shown.

For the ${\rm U}(1)$ parent Hamiltonian, starting with the lowest-energy
configuration with $\theta_{i}=0$ everywhere (with energy $E_{0}$), we may
perform the transformation

$\displaystyle\theta_{i}\rightarrow\theta_{i}+\frac{2\pi}{p}(-1)^{i_{x}+i_{y}}\left(\begin{array}[]{c}i_{x}+i_{y}\\\
i_{y}\end{array}\right)$ (44)

within a region $\mathcal{R}_{(0,0)}$ which consists of the points contained
within an equilateral triangle of side-length $p^{k}$ with the top of the
triangle at the origin $(i_{x},i_{y})=(0,0)$, i.e.

$\displaystyle\mathcal{R}_{(0,0)}\equiv\left\\{(i_{x},i_{y})\,|\,0\leq
i_{y}\leq i_{x}\,\mathrm{and}\,i_{x}+i_{y}\in[0,p^{k}-1]\right\\}$

where $p$ is a prime number and $k\geq 0$ is an integer. This transformation
creates three excitations – downward-facing triangles at which
$\theta_{i}+\theta_{i+\hat{x}}+\theta_{i+\hat{x}-\hat{y}}\neq 0$ mod $2\pi$ –
at each corner of the triangular region $\mathcal{R}$. The energy of the final
state is

$\displaystyle E_{0}+3t\left[1-\cos\left(\frac{2\pi}{p}\right)\right]$ (45)

.

We may demonstrate this as follows. The transformation (44) creates an
excitation at the downward-facing triangle whose bottom corner sits at the
origin since
$\theta_{-\hat{x}}+\theta_{0}+\theta_{-\hat{y}}\rightarrow\theta_{-\hat{x}}+\theta_{0}+\theta_{-\hat{y}}+(2\pi/p)$
under the transformation (44). Identical excitations are created at the two
downward-facing triangles at the two other corners of $\mathcal{R}_{(0,0)}$.

We now demonstrate that the transformation (44) acts trivially everywhere away
from the corners of $\mathcal{R}_{(0,0)}$. First, this transformation does not
create any excitations at downward-facing triangles that are contained
entirely within the region $\mathcal{R}_{(0,0)}$, i.e.
$\theta_{i}+\theta_{i+\hat{x}}+\theta_{i+\hat{x}-\hat{y}}$ is unchanged by the
transformation for any triangle for which $i$, $i+\hat{x}$, and
$i+\hat{x}-\hat{y}$ are contained entirely in $\mathcal{R}_{(0,0)}$. This can
be verified explicitly by applying the transformation (44).

Now, we may consider downward-facing triangles that border an edge of
$\mathcal{R}_{(0,0)}$. Consider downward-facing triangles that border the
bottom of region $\mathcal{R}_{(0,0)}$, as shown in Fig. 4. The vertices of
these triangles are located at points $(i,i+\hat{x},i+\hat{x}-\hat{y})$ for
which $i$, $i+\hat{x}-\hat{y}\in\mathcal{R}$ while
$i+\hat{x}\notin\mathcal{R}$. It can be seen that for these downward facing
triangles, the transformation (44) leads to the shift

$\displaystyle\theta_{i}+\theta_{i+\hat{x}}+\theta_{i+\hat{x}-\hat{y}}$
$\displaystyle\rightarrow\theta_{i}+\theta_{i+\hat{x}}+\theta_{i+\hat{x}-\hat{y}}+\frac{2\pi}{p}(-1)^{i_{x}+i_{y}}\left[\left(\begin{array}[]{c}i_{x}+i_{y}\\\
i_{y}\end{array}\right)+\left(\begin{array}[]{c}i_{x}+i_{y}\\\
i_{y}-1\end{array}\right)\right]=\theta_{i}+\theta_{i+\hat{x}}+\theta_{i+\hat{x}-\hat{y}}+\frac{2\pi}{p}\left(\begin{array}[]{c}p^{k}\\\
i_{y}\end{array}\right)$ (52)

In the last equality, we’ve used the fact that $i_{x}+i_{y}=p^{k}-1$ for the
downward-facing triangles bordering the bottom edge of $\mathcal{R}_{(0,0)}$.
We now observe that

$\displaystyle\left(\begin{array}[]{c}p^{k}\\\
i_{y}\end{array}\right)=0\mod\hskip 3.61371ptp$ (55)

for any $1\leq i_{y}\leq p^{k}-1$. This may be proven by observing that

$\displaystyle i_{y}\left(\begin{array}[]{c}p^{k}\\\
i_{y}\end{array}\right)=p^{k}\left(\begin{array}[]{c}p^{k}-1\\\
i_{y}-1\end{array}\right)$ (60)

The prime $p$ divides the right-hand-side of (60) at least $k$ times.
Furthermore, since $1\leq i_{y}\leq p^{k}-1$, $p$ divides $i_{y}$ at most
$k-1$ times. As a result, $p$ must divide $\left(\begin{array}[]{c}p^{k}\\\
i_{y}\end{array}\right)$ at least once, which completes the proof. From this,
we conclude that $\theta_{i}+\theta_{i+\hat{x}}+\theta_{i+\hat{x}-\hat{y}}$ is
shifted by an integer multiple of $2\pi$ for downward-facing triangles that
border the bottom edge of $\mathcal{R}_{(0,0)}$. Therefore, when starting from
the state with $\theta_{i}=0$ everywhere, only three excitations are created
by the transformation (44) which lie at the corners of $\mathcal{R}_{(0,0)}$,
and the energy cost of the final state is indeed given by (45).

We conclude by determining the ground-state degeneracy of the $Z_{p}$ models
on an $L\times L$ triangular lattice with periodic boundary conditions, for
particular values of $L$. First, the previous discussion of the Pascal’s
symmetry implies that if we start from a configuration with $\theta_{i}=0$
everywhere and perform a symmetry transformation over the region
$\mathcal{R}_{(0,0)}$, then the configuration of $\theta_{i}$ along the row
defined by $i_{x}+i_{y}=p^{k}-1$ of the lattice is given by

$\displaystyle\theta_{i}\,\,\mathrm{mod}\,\,2\pi=\left\\{\begin{array}[]{cc}\displaystyle(-1)^{i_{y}}\frac{2\pi}{p}&\,\,\,\,\,\,i_{y}\in[0,p^{k}-1]\\\
&\\\ 0&\mathrm{otherwise}\end{array}\right.$ (64)

Now consider another symmetry transformation along the region
$\mathcal{R}_{(-1,1)}$, i.e. along the equilateral triangle of side length
$p^{k}$ whose top vertex lies at the site $(i_{x},i_{y})=(-1,1)$. Superposing
this transformation with the previous transformation along
$\mathcal{R}_{(0,0)}$, we find that along the row $i_{x}+i_{y}=p^{k}-1$, that

$\displaystyle\theta_{i}\,\,\mathrm{mod}\,\,2\pi=\left\\{\begin{array}[]{cc}\displaystyle\frac{2\pi}{p}&\,\,\,\,\,i_{y}=0\,\,\mathrm{or}\,\,p^{k}\\\
&\\\ 0&\mathrm{otherwise}\end{array}\right.$ (68)

If we now impose periodic boundary conditions so that
$(i_{x},i_{y})\sim(i_{x}+L,i_{y})\sim(i_{x},i_{y}+L)$ where $L=p^{k}-1$, we
observe that the symmetry transformations along $\mathcal{R}_{(0,0)}$ and
$\mathcal{R}_{(-1,1)}$ have generated a new ground-state of the system, since
the row $i_{x}+i_{y}=p^{k}-1$ is now identified with row $i_{x}+i_{y}=0$, and
the configuration (68) is precisely the configuration of $\theta_{i}$ along
the first row $i_{x}+i_{y}=0$.

By applying additional “pairs” of these symmetry transformations (e.g. along
$\mathcal{R}_{(i_{x},-i_{x})}$ and
$\mathcal{R}_{(i^{\prime}_{x},-i^{\prime}_{x})}$), we may generate other
ground-state configurations. Each unique ground-state corresponds to a unique
configuration of $\theta_{i}=\frac{2\pi}{p}m_{i}$ for $m_{i}\in[0,\ldots,p-1]$
at each site along the row $i_{x}+i_{y}=0$. Not all of these $p^{L}$
configurations can be reached by starting with the ground-state with
$\theta_{i}=0$ everywhere and applying the fractal symmetry transformation in
pairs. It is clear, however, that by continuing to apply these symmetry
transformations in pairs, any configuration of $\theta_{i}$ along the row
$i_{x}+i_{y}=0$ can be reached for which

$\displaystyle\sum_{i|i_{x}+i_{y}=0}(-1)^{i_{y}}m_{i}=0\hskip
7.22743pt\mathrm{mod}\,\,p$ (69)

There are precisely $p^{L-1}$ such configurations. We conclude that the
ground-state degeneracy of the $Z_{p}$ model is $p^{L-1}$ on an $L\times L$
lattice with periodic boundary conditions when $L=p^{k}-1$.

##  Appendix C Three-point correlation of the $Z_{p}$ and $Z_{N}$ models

Starting from the classical $Z_{p}$ Pascal’s triangle models in Eq. 8, since
the $Z_{p}$ clock degrees of freedom are primitive $p$-roots of unity for
which $\sigma^{\dagger}=\sigma^{p-n}$,
$\mathop{\mathrm{Re}}{\sigma}=\mathop{\mathrm{Re}}{\sigma^{p-n}}$,
$\mathop{\mathrm{Im}}{\sigma}=-\mathop{\mathrm{Im}}{\sigma^{p-n}}$ and
likewise for the dual $\tau$-variables. As a result, for the partition
function of a single $\tau$-variable $Z_{1}=\sum_{n=0}^{p-1}e^{\beta
t\cos(\frac{2\pi n}{p})}$:

$\displaystyle\langle\mathop{\mathrm{Im}}{\tau^{m}}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{1}}\sum_{n=0}^{p-1}\sin\bigg{(}\frac{2\pi
mn}{p}\bigg{)}e^{\beta t\cos(\frac{2\pi n}{p})},$ (70) $\displaystyle=$
$\displaystyle 0;$ (72)
$\displaystyle\langle\mathop{\mathrm{Re}}{\tau^{m}}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{1}}\sum_{n=0}^{p-1}\cos\bigg{(}\frac{2\pi
mn}{p}\bigg{)}e^{\beta t\cos(\frac{2\pi n}{p})},$ (73) $\displaystyle=$
$\displaystyle\frac{1}{Z_{1}}\sum_{n=0}^{p-1}\cos\bigg{(}\frac{2\pi
m(p-n)}{p}\bigg{)}e^{\beta t\cos(\frac{2\pi n}{p})},$ (75) $\displaystyle=$
$\displaystyle\langle\mathop{\mathrm{Re}}{(\tau^{m})^{\dagger}}\rangle.$ (77)

As a consequence, $\langle\tau^{m}\rangle$ =
$\langle(\tau^{m})^{\dagger}\rangle$. In the dual representation, the three-
point correlation function factors since each $\tau$ is statistically
independent

$\displaystyle\mathcal{C}_{3}(r=p^{k})$ $\displaystyle=$
$\displaystyle\langle\sigma_{0,0}\sigma_{r,0}\sigma_{0,r}\rangle,$ (78)
$\displaystyle=$ $\displaystyle\langle\text{product of }\tau^{m}\text{ in
Pascal's triangle}\rangle,$ (80) $\displaystyle=$ $\displaystyle\text{product
of }\langle\tau^{m}\rangle\text{ in Pascal's triangle},$ (82) $\displaystyle=$
$\displaystyle\prod_{m=1}^{\frac{p-1}{2}}\langle\tau^{m}\rangle^{N_{m,p-m}(k)},$
(84)

where $N_{m,p-m}(k)$ is the number of times $m$ and $p-m$ appear in the
Pascal’s triangle modulo $p$ of length $p^{k}-1$. Note that the staggering of
$\tau$-variables in the representation of $\sigma$ is inconsequential for the
calculation of the three-point correlation as
$\langle\tau^{m}\rangle=\langle(\tau^{m})^{\dagger}\rangle$.

As a result, computing the three-point correlation function reduces to
counting the number of times each unique entry in a Pascal’s triangle modulo
$p$ appears. Since the elements of a Pascal’s triangle modulo $p$ are binomial
coefficients ${j\choose i}\mod p$, one may in principle calculate
$N_{m,p-m}(k)$ by brute force according to

$\displaystyle N_{m,p-m}(k)$ $\displaystyle=$ $\displaystyle\sum_{i\leq
j}\sum_{j=0}^{p^{k}-1}\delta\bigg{(}m,{j\choose i}\text{ mod }p\bigg{)}$ (85)
$\displaystyle+$ $\displaystyle\sum_{i\leq
j}\sum_{j=0}^{p^{k}-1}\delta\bigg{(}p-m,{j\choose i}\text{ mod }p\bigg{)}.$
(87)

However, this is neither efficient analytically nor numerically, even at small
$p$. It is instead desirable to derive recurrence relations that can be solved
assuming initial conditions $N_{m,p-m}(1)$ are provided which can be computed
according to Eq. 87 with $k=1$. In order to derive such recurrence relations,
we can make use that for $\cos(m\phi),m=1,\dots,\frac{p-1}{2}$ such that

$\displaystyle\cos(m\phi)=\cos(qn\phi)\Leftrightarrow\begin{cases}m=qn\text{
mod }p\\\ m=q(p-n)\text{ mod }p\end{cases}.$ (88)

This implies that $\langle\tau^{m}\rangle=\langle\tau^{qn}\rangle$ where
$q=q(m,n)$ is the smallest element of $Z_{p}$ for which one of the two
equations in Eq. 88 is solved. If we consider an arbitrary Pascal’s triangle
modulo $p$ of length $p^{k}-1$, this triangle is self-similar and can be
constructed by tiling the same Pascal’s triangle of length $p^{k-1}-1$
according to the Pascal’s triangle of length $p-1$. The number of
$\langle\tau^{n}\rangle$ which become $\langle\tau^{m}\rangle$ upon tiling is
$N_{q(m,n),p-q(m,n)}(1)$. As a result, we arrive at a set of $\frac{p-1}{2}$
coupled recurrence relations

$\displaystyle
N_{m,p-m}(k)=\sum_{n=1}^{\frac{p-1}{2}}N_{q(m,n),p-q(m,n)}(1)N_{n,p-n}(k-1)$
(89)

For small $p$, $N_{m,p-m}(1)$ can be calculated and these recurrence relations
can be solved without the use of numerics. For example, for the $p=5$ model

$\displaystyle q=\begin{bmatrix}1&2\\\ 2&1\end{bmatrix}.$

One can compute the initial conditions by hand from Eq. 87 to find
$N_{1,4}(1)=12$ and $N_{2,3}(1)=3$, yielding coupled recurrence relations

$\displaystyle N_{1,4}(k)$ $\displaystyle=$ $\displaystyle
12N_{1,4}(k-1)+3N_{2,3}(k-1),$ (90) $\displaystyle N_{2,3}(k)$
$\displaystyle=$ $\displaystyle 12N_{2,3}(k-1)+3N_{1,4}(k-1).$ (91)

This can be solved again using the initial conditions to recover

$\displaystyle N_{1,4}(k)=\frac{3^{k-1}}{10}(3\times 5^{k+1}+5\times
3^{k+1}),$ $\displaystyle N_{2,3}(k)=\frac{3^{k-1}}{10}(3\times
5^{k+1}-5\times 3^{k+1}).$

One can verify that $N_{1,4}(k)+N_{2,3}(k)$ must be the area of the fractal of
length $p^{k}-1$, $N_{1,4}(k)+N_{2,3}(k)=15^{k}=5^{k\times\frac{\ln 15}{\ln
5}}=5^{k\times d_{H}}$. The three-point correlation can now be computed as

$\displaystyle\mathcal{C}_{3}(r=5^{k})$ $\displaystyle=$
$\displaystyle\langle\tau\rangle^{N_{1,4}(k)}\langle\tau^{2}\rangle^{N_{2,3}(k)},$
(92) $\displaystyle=$
$\displaystyle\bigg{(}\frac{1}{Z_{1}}\sum_{n=0}^{p-1}\cos\bigg{(}\frac{2\pi
n}{p}\bigg{)}e^{\beta t\cos(\frac{2\pi n}{p})}\bigg{)}^{N_{1,4}(k)}$ (94)
$\displaystyle\times$
$\displaystyle\bigg{(}\frac{1}{Z_{1}}\sum_{n=0}^{p-1}\cos\bigg{(}\frac{4\pi
n}{p}\bigg{)}e^{\beta t\cos(\frac{2\pi n}{p})}\bigg{)}^{N_{2,3}(k)}$ (96)
$\displaystyle\approx$ $\displaystyle\bigg{(}\frac{e^{\beta
t}+0.62e^{0.31\beta t}-1.62e^{-0.81\beta t}}{e^{\beta t}+2e^{0.31\beta
t}+2e^{-0.81\beta t}}\bigg{)}^{N_{1,4}(k)}$ (98) $\displaystyle\times$
$\displaystyle\bigg{(}\frac{e^{\beta t}-1.62e^{0.31\beta t}+0.62e^{-0.81\beta
t}}{e^{\beta t}+2e^{0.31\beta t}+2e^{-0.81\beta t}}\bigg{)}^{N_{2,3}(k)}.$
(100)

For sufficiently large distances, $\mathcal{C}_{3}(r=5^{k})\sim
e^{-\alpha^{\prime}r^{\gamma}}$ where $\gamma$ is an exponent that is larger
than $d_{H}$, due to the decay of $\langle\tau^{2}\rangle$ being stronger than
$\langle\tau\rangle$. As a comparison, the three-point correlation in the
$Z_{5}$ model decays much faster than the three-point correlation of the
Sierpinski-triangle model and the $Z_{3}$ Pascal’s triangle model.

The multiple fractal symmetries of the $Z_{N}$ models has important
consequences for the behavior of the three-point correlations as well. All of
the $Z_{N}$ models retain the same dual description, and $\sigma$ can be
represented as an infinite staggered product of dual $\tau$-variables on a
Pascal’s triangle modulo $N$ in the thermodynamic limit. However, since $N$ is
composite, this Pascal’s triangle is not a fractal and hence there is no
product of three $\sigma$’s that has compact support causing
$\mathcal{C}_{3}(r)$ to vanish. Instead, for each unique prime divisor $p$ of
$N$ for which $qp=N$, $\sigma^{q}$ can be represented as an infinite staggered
product of $\tau^{q}$ in the shape of a Pascal’s triangle modulo $p$.
Therefore,
$\mathcal{C}_{3}^{q}(r)=\langle\sigma^{q}_{0,0}\sigma^{q}_{r,0}\sigma^{q}_{0,r}\rangle$
is non-vanishing when $r=p^{k}$. Using $Z_{N=6}$ as an example again,
$\mathcal{C}^{3}_{3}=\langle\sigma^{3}_{0,0}\sigma^{3}_{r,0}\sigma_{0,r}^{3}\rangle$
characterizes SSB of the $Z_{2}$ fractal symmetry while
$\mathcal{C}^{2}_{3}=\langle\sigma^{2}_{0,0}\sigma^{2}_{r,0}\sigma_{0,r}^{2}\rangle$
characterizes SSB of the $Z_{3}$ fractal symmetry. Note that if $N$ is taken
to be prime, since the only prime-factor of $N$ is itself, these results
reduce down to those of the $Z_{p}$ models and
$\langle\sigma_{0,0}\sigma_{r,0}\sigma_{0,r}\rangle$ is ordered.

Figure 5: The dual honeycomb lattice where Eq. 106 is defined on. Each vertex
of the honeycomb lattice coincides with either a downward or upward facing
plaquette of the original triangular lattice.

##  Appendix D A related model with stable $\sqrt{3}\times\sqrt{3}$ order

We consider a model directly related to Eq. 20, but instead is stable against
the inclusion of quantum fluctuations. In the ${\rm U}(1)$ parent model, the
Pascal’s triangle fractal symmetries only exist because a single orientation
of triangular plaquettes is being summed over. If we include the sum over
three-boson interactions on upward facing plaquettes as well, the model
becomes

$\displaystyle
H=\sum_{\bigtriangledown,\triangle}-t\cos\left(\theta_{1}+\theta_{2}+\theta_{3}\right)+\sum_{j}\frac{U}{2}n_{j}^{2}.$
(101)

While this model does not have the Pascal’s triangle fractal symmetries of Eq.
20, it still retains the global ${\rm U}(1)_{1}\times{\rm U}(1)_{2}$ symmetry
in Eq. 5 and hence the ground state of this model at small $U$ will still have
the $\sqrt{3}\times\sqrt{3}$ order. In this case, unlike Eq. 20, the boson
ordered phase will actually be stable against quantum fluctuations. This can
be seen by analyzing the dual of Eq. 101. Now, the dual variables are defined
on the sites $\bar{i}$ of the dual honeycomb lattice (Fig. 5):

$\displaystyle\theta_{1}+\theta_{2}+\theta_{3}=\phi_{\bar{i}},\ \ \
\mathrm{site}\ 1,2,3\ \mathrm{around}\ \bar{i};$ (102) (103) $\displaystyle
n_{j}=\sum_{\bar{i}}\Psi_{\bar{i}},\ \ \ \bar{i}=1,\cdots 6\ \mathrm{around}\
j.$ (104)

The dual variable is subject to a local gauge constraint around each site $j$:

$\displaystyle\sum_{\bar{i}\ \mathrm{around}\
j}(-1)^{\bar{i}}\phi_{\bar{i}}=0\ \mathrm{mod}\ 2\pi.$ (105)

Again, $\phi_{\bar{i}}$ are compact variables, while $\Psi_{\bar{i}}$ takes
discrete values. The dual Hamiltonian of Eq. 101 reads Xu and Vishwanath
(2006)

$\displaystyle
H_{d}=\sum_{\bar{j}}-t\cos(\phi_{\bar{j}})+\sum_{j}\frac{U}{2}\left(\sum_{\bar{i}\
\mathrm{around}\ j}\Psi_{\bar{i}}\right)^{2}.$ (106)

The key difference from the model considered in the main text is that, this
dual model Eq. 106 has a gauge symmetry, i.e. the dual Hamiltonian is
invariant under the following local gauge transformation:

$\displaystyle\Psi_{\bar{i}}\rightarrow\Psi_{\bar{i}}+(-1)^{\bar{i}}\epsilon_{j},\
\ \ \bar{i}\ \mathrm{around}\ j,$ (107)

where $\epsilon_{j}$ is an arbitrary integer-valued function on the original
triangular lattice $j$.

Now we can again replace $\Psi_{\bar{i}}$ by continuous variables, and the
constraint that $\Psi$ takes discrete values will be enforced through “vertex
operators” such as $\cos(2\pi\Psi)$. If we ignore these vertex operators, the
Lagrangian of the Gaussian theory of Eq. 106 reads:

$\displaystyle\mathcal{L}_{d}=\sum_{\bar{i}}\frac{1}{2t}(\partial_{\tau}\Psi_{\bar{i}})+\sum_{j}\frac{U}{2}\left(\sum_{\bar{i}\
\mathrm{around}\ j}\Psi_{\bar{i}}\right)^{2}.$ (108)

This Gaussian theory has a continuous gauge symmetry, which corresponds to
allowing $\epsilon_{j}$ to take arbitrary continuous values in Eq. 107. Then,
the Gaussian fixed point of Eq. 108 is stable against weak vertex operators,
since the vertex operators break the continuous gauge symmetry of the Gaussian
theory Eq. 108.