# Phase diagram and critical properties of a two-dimensional associating
lattice gas

I. Ibagon<EMAIL_ADDRESS>Departamento de Física, ICEx, Universidade
Federal de Minas Gerais, C. P. 702, 30123-970, Belo Horizonte, Minas Gerais -
Brazil  A. P. Furlan<EMAIL_ADDRESS>Departamento de Física, ICEx,
Universidade Federal de Minas Gerais, C. P. 702, 30123-970, Belo Horizonte,
Minas Gerais - Brazil  T. J. Oliveira<EMAIL_ADDRESS>Departamento de Física,
Universidade Federal de Viçosa, 36570-900, Viçosa, Minas Gerais, Brazil  R.
Dickman<EMAIL_ADDRESS>Departamento de Física and National Institute
of Science and Technology for Complex Systems, ICEx, Universidade Federal de
Minas Gerais, C. P. 702, 30123-970 Belo Horizonte, Minas Gerais - Brazil

(August 27, 2024)

###### Abstract

We revisit the associating lattice gas (ALG) introduced by Henriques et al.
[PRE 71, 031504 (2005)] in its symmetric version. In this model, defined on
the triangular lattice, interaction between molecules occupying nearest-
neighbor sites depends on their relative orientation, mimicking the formation
of hydrogen bonds in network-forming fluids. Although all previous studies of
this model agree that it has a disordered fluid (DF), a low-density liquid
(LDL) and a high-density liquid (HDL) phase, quite different forms have been
reported for its phase diagram. Here, we present a thorough investigation of
its phase behavior using both transfer matrix calculations and Monte Carlo
(MC) simulations, along with finite-size scaling extrapolations. Results in
striking agreement are found using these methods. The critical point
associated with the DF-HDL transition at full occupancy, identified by Furlan
and coworkers [Phys. Rev. E 100, 022109 (2019)] is shown to be one terminus of
a critical line separating these phases. In opposition to previous simulation
studies, we find that the transition between the DF and LDL phases is always
discontinuous, similar to the LDL-HDL transition. The associated coexistence
lines meet at the point where the DF-HDL critical line ends, making it
critical-end-point. Overall, the form of the phase diagram observed in our
simulations is very similar to that found in the exact solution of the model
on a Husimi lattice. Our results confirm that, despite the existence of some
waterlike anomalies in this model, it is unable to reproduce key features of
the phase behavior of liquid water.

###### pacs:

68.35.Rh

##  I Introduction

The importance water has to life on Earth is perhaps comparable to the degree
of unusual behavior this liquid displays. In fact, when compared with other
substances, liquid water presents a number of anomalies (such as in density,
diffusivity and so on), which are believed to originate in the particular
manner in which H2O molecules interact, through hydrogen bonds forming a
local, approximately tetrahedral network Stanley (2013). The three main
scenarios proposed to understand water, viz. the stability-limit Speedy and
Angell (1976), the singularity-free Sastry  _et al._ (1996) and liquid-liquid
critical point (LLCP) Poole  _et al._ (1992) hypotheses, are all consistent
with the presence of such a network. The later scenarios are based on the fact
that just as ice becomes glassy and can exist in low- and high-density
amorphous phases Mishima, Calvert, and Whalley (1985), the continuation of
these phases ( increasing the temperature $T$ above that of spontaneous
crystallization) may give rise to a low-density liquid (LDL) and a high-
density liquid (HDL) phase Stanley (2013); Lynden-Bell  _et al._ (2010).
Although there is some indirect experimental evidence for the existence of the
LLCP Soper and Ricci (2000); Mishima and Stanley (1998); Woutersen  _et al._
(2018), including some coming from water in nanopores Liu  _et al._ (2005);
Chen  _et al._ (2006); Mallamace  _et al._ (2007); Mallamace (2009), a
definitive proof is still lacking Soper (2019), mainly because the LLCP would
exist in the metastable, deeply supercooled region, which is very difficult to
access experimentally in bulk water.

In this context, much of this debate has been based on computational studies,
with several realistic models for water confirming the LLCP hypothesis (see
Palmer  _et al._ (2018) for a recent review). A number of simplified
associating lattice gas (ALG) models, designed to possess two liquid phases,
have also been introduced in two- Bell and Lavis (1970a, b); Huckaby and Hanna
(1987); Patrykiejew, Pizio, and Sokołowski (1999); Buzano  _et al._ (2004);
Henriques and Barbosa (2005); Balladares, Henriques, and Barbosa (2007) and
three-dimensions Bell (1972); Meijer, Kikuchi, and Royen (1982); Lavis and
Southern (1984); Besseling and Lyklema (1994); Roberts and Debenedetti (1996);
Roberts, Panagiotopoulos, and Debenedetti (1996); Pretti and Buzano (2004,
2005); Sastry, Sciortino, and Stanley (1993); Borick, Debenedetti, and Sastry
(1995); Sastry  _et al._ (1996); Girardi  _et al._ (2007). These ALG models
usually feature an attractive, orientation-dependent interaction between
molecules occupying nearest neighbor (NN) sites — mimicking the hydrogen
bonding in network forming fluids —, competing with a repulsive excluded-
volume interaction. Phase diagrams displaying a LLCP have been reported for
some of these models Borick, Debenedetti, and Sastry (1995); Roberts and
Debenedetti (1996); Roberts, Panagiotopoulos, and Debenedetti (1996); Pretti
and Buzano (2004, 2005); Henriques and Barbosa (2005); Girardi  _et al._
(2007); Balladares, Henriques, and Barbosa (2007), though latter works on some
these models have demonstrated a phase behavior without a LLCP Henriques and
Barbosa (2005); Girardi  _et al._ (2007); Balladares, Henriques, and Barbosa
(2007).

Figure 1: Qualitative phase diagrams for the symmetric H&B model in the
(reduced) $\bar{\mu}-\bar{T}$ plane adapted from references Balladares,
Henriques, and Barbosa (2007); Oliveira, Stilck, and Barbosa (2010); Furlan,
Fiore, and Barbosa (2015) in panels (a), (b), and (c), respectively. Solid and
dashed lines represent continuous and discontinuous phase transitions,
respectively. The coexistence points at $T=0$ are exact, while the loci of the
critical points $c_{1}$ and $c_{2}$, of the triple point $tp$, and of the
tricritical points $tc_{1}$ and $tc_{2}$ are those found numerically in the
references cited above.

An example of this is the symmetric version of the Besseling and Lyklema model
Besseling and Lyklema (1994), firstly analyzed by Girardi et al. Girardi  _et
al._ (2007). This model is defined on the body-centered cubic lattice, where
each molecule can form up to four hydrogen bonds with NNs, similar to H2O.
While these bonds lower their energy, nonbonding molecules at NN sites
experience a repulsive “van der Waals” interaction. Beyond the GAS-liquid
coexistence line ending at a critical point, a LDL-HDL coexistence line ending
at a second critical point was also reported in Girardi  _et al._ (2007),
yielding a phase diagram qualitatively similar to the one in Fig. 1(a). (As an
aside, we remark that the phase commonly called GAS in the literature related
to these models will be regarded as a disordered fluid (DF) here, since it may
extend to large chemical potentials and so, high densities.) Subsequent
studies of this ALG Buzzano, De Stefanis, and Pretti (2008); Pretti, Buzzano,
and De Stefanis (2009); Szortyka  _et al._ (2010), notwithstanding, revealed a
much more complex phase behavior, with the LLCP [$c_{2}$ in Fig. 1(a)] being
actually a tricritical point, while $c_{1}$ is disputed as being a critical-
end-point Buzzano, De Stefanis, and Pretti (2008); Pretti, Buzzano, and De
Stefanis (2009) or a bicritical point Szortyka  _et al._ (2010). These works
also uncovered a critical line separating the HDL and DF phases Buzzano, De
Stefanis, and Pretti (2008); Pretti, Buzzano, and De Stefanis (2009); Szortyka
_et al._ (2010).

A similar sequence of events marks the history of the ALG model introduced by
Henriques and Barbosa (H&B) Henriques and Barbosa (2005) which is a
nonsymmetric two-dimensional (2D) version of the 3D model analyzed by Girardi
et al. Girardi  _et al._ (2007). The H&B model is defined on the triangular
lattice, in which each molecule has four bonding arms (two donors and two
acceptors of protons), beyond two inert arms separated by 180°. In their
original Monte Carlo analysis, H&B Henriques and Barbosa (2005) reported a
phase diagram analogous to that in Fig. 1(a), with two coexistence lines and
two critical points. Further work on this model indicated, however, a
different and richer phase diagram, with two additional continuous transition
lines: DF-LDL and DF-HDL Szortyka  _et al._ (2009). In this phase diagram, the
DF-LDL coexistence line ends at a tricritical point where it meets the DF-LDL
continuous line; and the LDL-HDL coexistence line ends at a bicritical point
where it meets the continuous DF-LDL and DF-HDL lines Szortyka  _et al._
(2009).

Several modifications of the H&B model have been proposed Balladares,
Henriques, and Barbosa (2007); Furlan, Almarza, and Barbosa (2016); Furlan,
Fiore, and Barbosa (2015). Of particular interest here is the symmetric
version introduced by Balladares et al. Balladares, Henriques, and Barbosa
(2007), in which no distinction is made between donor and acceptor arms. The
literature on this symmetric ALG is also marked by controversy. For example,
the original MC simulations of Ref. Balladares, Henriques, and Barbosa (2007)
again suggested the phase behavior of Fig. 1(a), consistent with the LLCP
hypothesis. The semi-analytical mean-field solution of this simplified version
of the H&B model on a Husimi lattice built with hexagonal plaquettes Oliveira,
Stilck, and Barbosa (2010) revealed a qualitatively different phase diagram,
with three coexistence lines (DF-HDL, DF-LDL and LDL-HDL) meeting at a triple
point ($tp$), as shown in Fig. 1(b). Later MC simulations Furlan, Fiore, and
Barbosa (2015) supported the phase behavior of Fig. 1(c), with two coexistence
lines (DF-LDL and LDL-HDL) ending at two tricritical points ($tc_{1}$ and
$tc_{2}$), which are linked by a critical DF-LDL line. No transition between
the HDL phase and the disordered fluid was reported in Ref. Furlan, Fiore, and
Barbosa (2015).

More recently, this symmetric H&B model was investigated in the limit of full
occupancy (corresponding to a chemical potential $\mu\to\infty$) using MC
simulations and semi-analytical solutions on Husimi lattices Furlan  _et al._
(2019). The simulations yielded a continuous transition between the HDL and
the DF phase at a (reduced) temperature $\bar{T}_{c}=0.9525(5)$, with critical
exponents and Binder cumulant in good agreement with those of the three-state
Potts model in 2D Wu (1982). On the other hand, the Husimi lattice solutions
furnish a discontinuous transition, as a continuation of the DF-HDL
coexistence line found for finite $\mu$Oliveira, Stilck, and Barbosa (2010) to
the limit of $\mu\rightarrow\infty$. All these results raise important
questions on the phase behavior of this model: i) Does it indeed exhibit a
LLCP as in Fig. 1(a)? ii) What is the true nature of the DF-LDL transition? Is
it discontinuous as in Fig. 1(b) or does it have a tricritical point and a
continuous part like in Fig. 1(c)? iii) Assuming that a transition line does
indeed exist separating the DF and HDL phases: Is it always discontinuous but
ends at a critical point in the limit $\mu\to\infty$? Or does it start
discontinuous and change to continuous at some finite $\mu$? Or, even, is the
entire line continuous?

In this work we address these points through an extensive numerical study of
this ALG model using both transfer-matrix analysis and MC simulations in the
grand-canonical ensemble. Consistent results are obtained from both methods,
demonstrating that the correct topology of the phase diagram is that of Fig.
1(b). However, in contrast with the mean-field result Oliveira, Stilck, and
Barbosa (2010) shown in this figure, we find that the transition between the
DF and HDL phases is always continuous and the critical DF-HDL line ends at a
critical-end-point, where it meets the DF-LDL and LDL-HDL coexistence lines.
The critical properties of the continuous transition line — an interesting
problem scarcely tackled in the literature of lattice models for water (see
Furlan  _et al._ (2019) for a discussion) — are also investigated and an
intriguing behavior is found, with some quantities presenting a close
agreement with the three-state Potts class, while others have a large
deviation from it.

The paper is organized as follows. In Sec. II, we present the investigated
model. The transfer-matrix method and the associated results are discussed in
Sec. III. A similar thing is done in Sec. IV for the Monte Carlo simulations.
The comparison of the outcomes from both methods and our concluding remarks
are presented in Sec. V.

##  II Model

We study the symmetric version of the Henriques and Barbosa (H&B) Henriques
and Barbosa (2005) associating lattice gas (ALG) model, investigated in Refs.
Balladares, Henriques, and Barbosa (2007); Oliveira, Stilck, and Barbosa
(2010); Furlan, Fiore, and Barbosa (2015); Furlan  _et al._ (2019). The model
is defined on the triangular lattice (coordination number $z=6$), where each
site $i$ can be either empty ($\sigma_{i}=0$) or occupied by a single molecule
($\sigma_{i}=1$). A chemical potential $\mu$ is associated to each particle in
the system, since we will work in the grand canonical ensemble. Whenever two
molecules occupy nearest-neighbor (NN) sites, they experience a repulsive “van
der Walls” interaction, gaining an energy $\varepsilon>0$. In addition, two NN
molecules can form a hydrogen bond, which decreases their energy by a factor
$-\gamma<0$, provided that their bonding arms point toward each other. In
fact, a key feature of the model analyzed here (and of other ALGs as well) is
the existence of orientational degrees of freedom and related interactions
introduced by bonding arms. In our case, each molecule has four bonding arms
and two inert (nonbonding) arms, as shown in Fig. 2. The four bonding arms are
all equivalent and always interact with those of NN molecules in the same
manner, i.e., there is no distinction between donors and acceptors (of
protons) here, in opposition to the original (non-symmetric) H&B model
Henriques and Barbosa (2005). The inert arms of a given molecule are in
diametrically opposed lattice edges (forming an angle 180°), which results in
three orientational states for each molecule (see Fig. 2.) Following Ref.
Furlan  _et al._ (2019), we may conveniently define an orientational variable,
${\bm{\eta}}_{i}$, pointing along one of the nonbonding directions of the
molecule at site $i$. In this fashion, it is simply given by the generators of
the triangular lattice, $\hat{e}_{k}$ with $k=1,2$ and 3, which are also
depicted in Fig. 2. Namely, ${\bm{\eta}}_{i}\equiv\hat{e}_{k}$, where
$\hat{e}_{1}={{\bm{\hat{\textnormal{\i}}}}}$,
$\hat{e}_{2}=+\frac{1}{2}{{\bm{\hat{\textnormal{\i}}}}}+\frac{\sqrt{3}}{2}{{\bm{\hat{\textnormal{\j}}}}}$
and
$\hat{e}_{3}=-\frac{1}{2}{{\bm{\hat{\textnormal{\i}}}}}+\frac{\sqrt{3}}{2}{{\bm{\hat{\textnormal{\j}}}}}$.

With these definitions, we may write the Hamiltonian of this ALG model as,

${\cal H}\left({\bm{\eta}},{\bm{r}}\right)=\sum_{\left\langle
i,j\right\rangle}\sigma_{i}\sigma_{j}\left[\varepsilon+u_{ij}\right]-\mu\sum_{i}\sigma_{i},$
(1)

where ${\bm{r}}={\bm{r}}_{i}-{\bm{r}}_{j}$ is a unit vector in the set
$\left\\{\pm\hat{e}_{1},\pm\hat{e}_{2},\pm\hat{e}_{3}\right\\}$. The second
sum gives the total number of particles, $N$, while the first sum runs over
all NN pairs $(i,j)$, whose orientational interaction is enconded in variable
$u_{ij}$, given by Furlan  _et al._ (2019)

$u_{ij}=-\gamma\left(\frac{4}{3}\right)^{2}\left[1-\left({\bm{\eta}}_{i}\cdot{\bm{r}}\right)^{2}\right]\left[1-\left({\bm{\eta}}_{j}\cdot{\bm{r}}\right)^{2}\right],$
(2)

Note that whenever ${\bm{r}}$ is parallel or antiparallel to ${\bm{\eta}}_{i}$
and/or ${\bm{\eta}}_{j}$, molecules $(i,j)$ do not form a hydrogen bond and
$u_{i,j}=0$; otherwise $u_{i,j}=-\gamma$. Following most of the previous
studies of the ALG model Balladares, Henriques, and Barbosa (2007); Oliveira,
Stilck, and Barbosa (2010); Furlan, Fiore, and Barbosa (2015), we will use the
parameter $\varepsilon$ as the energy scale and set $\gamma/\varepsilon=2$. In
this way, the net interaction between two NN molecules forming (not forming) a
hydrogen bond is attractive (repulsive). We therefore have as free parameters
the dimensionless temperature $\bar{T}=k_{B}T/\varepsilon$ and chemical
potential $\bar{\mu}=\mu/\varepsilon$. We remark that in Ref. Furlan  _et al._
(2019) the temperature (called $\tau$ there) was defined in units of $\gamma$,
so that $\tau=\bar{T}/2$.

Figure 2: Possible orientational states ${\bm{\eta}}\equiv\hat{e}_{k}$, with
$k=1,2$ and 3. Bonding arms are represented by thick black lines and the
generators of the triangular lattice $\hat{e}_{1}$, $\hat{e}_{2}$ and
$\hat{e}_{3}$ are indicated by tiny coloured arrows. For this configuration no
hydrogen bond is formed between these molecules.

As already discussed in the Introduction [Sec. I], from previous studies of
the ALG model Balladares, Henriques, and Barbosa (2007); Oliveira, Stilck, and
Barbosa (2010); Furlan, Fiore, and Barbosa (2015), one knows that its phase
diagram exhibits three phases: disordered fluid (DF), low-density liquid (LDL)
and high-density liquid (HDL). In the ground state ($\bar{T}=0$), the DF phase
has particle density $\rho=0$ (the reason it is usually called GAS) and is the
most stable phase for $\bar{\mu}<-2$. For $-2<\bar{\mu}<2$ the LDL phase is
the one with the lowest free energy. In this phase, one of the four
sublattices is empty while each of the other three is occupied by molecules in
one of the orientational states, such that each molecule participates in four
hydrogen bonds [see Fig. 3 (a)]. Hence, the LDL phase has $\rho=3/4$ and a
fourfold degeneracy. For $\bar{\mu}>2$ the HDL phase is the most stable, with
the lattice fully occupied ($\rho=1$) and all inert arms aligned along the
same direction, as illustrated in Fig. 3(b). The HDL phase is threefold
degenerate. It is quite simple to show that the DF-LDL and LDL-HDL phases
coexist at $\bar{\mu}=-2$ and $\bar{\mu}=2$, respectively, when $\bar{T}=0$
Balladares, Henriques, and Barbosa (2007); Oliveira, Stilck, and Barbosa
(2010).

Figure 3: Illustration of a ground state configuration of the (a) low-density
liquid (LDL) phase and (b) high-density liquid (HDL) phase. Different colors
represent molecules in distinct orientational states. The bonding arms are
indicated by thick lines.

##  III Transfer-matrix analysis

###  III.1 Preliminaries

We define the transfer matrix, $\tilde{M}_{L}$, of the model on strips of
infinite length and finite width $L$. As illustrated in Fig. 4, we consider
square lattice strips with diagonal edges, with periodic boundary conditions
in the finite (horizontal) direction. As usual, each possible state for a
given row of the strip will be associated with a row and a column of
$\tilde{M}_{L}$. Since each lattice site can be in four different states, we
have a total of $4^{L}$ states for a strip of width $L$, so that
$\tilde{M}_{L}$ has dimension $4^{L}\times 4^{L}$. In order to investigate the
full-occupancy case ($\bar{\mu}\rightarrow\infty$), where only three states
are possible for each site, we can work with reduced matrices of dimension
$3^{L}\times 3^{L}$. In any case, the element $(i,j)$ of $\tilde{M}_{L}$ is
obtained by setting the row $n$ in state $i$, the row $n+1$ in state $j$ and
determining the statistical weight of this configuration. Then,
$\tilde{M}_{L}(i,j)=\exp[(N_{p}\bar{\mu}+2N_{hb}-N_{vdW})/2\bar{T}]$, where
$N_{p}$ is the total number of particles on both rows (note that we are
dividing this by 2), and $N_{hb}$ and $N_{vdW}$ counts the effective number of
hydrogen bonds and van der Waals interactions (i.e, first neighbors) in these
rows, respectively. By “effective”, we mean that we are counting once (twice)
interactions along the horizontal (diagonal and vertical) direction and, then,
dividing them by 2. In this way, the “internal” interactions, along the
vertical and diagonal directions, are counted once, while each one in the
horizontal contribute with $1/2$ in each state, since they appear in two
states.

Figure 4: Illustration of the approach to the triangular lattice used in our
transfer-matrix calculations. A strip of width $L=5$ is shown.

It is hard to deal directly with $\tilde{M}_{L}$, because all of its elements
are non-null. Therefore, in order to save computer RAM memory and speed up the
calculations, we replace the procedure of adding an entire row, by the
procedure of adding a single elementary cell (formed by squares with a
diagonal edge). Working in this way, we have to deal with $L+1$ matrices,
$\tilde{T}_{L}^{(s)}$, so that
$\tilde{M}_{L}=\tilde{T}_{L}^{(L)}\cdot\tilde{T}_{L}^{(L-1)}\cdots\tilde{T}_{L}^{(0)}$.
Here, $\tilde{T}_{L}^{(0)}$ creates a new site in row $n+1$ above a complete
row $n$; the matrices $\tilde{T}_{L}^{(s)}$, with $1\leqslant s\leqslant L-1$,
add new elementary cells to row $n+1$; and finally $\tilde{T}_{L}^{(L)}$
completes the row $n+1$ by imposing the appropriated (periodic here) boundary
conditions (see, e.g, Refs. Guo and Blöte (2002); Oliveira and Stilck (2015)
for more details). The advantage of this method relies on the fact that the
matrices $\tilde{T}_{L}^{(s)}$ are very sparse. In fact, the total number of
non-null terms summing over all matrices $\tilde{T}_{L}^{(s)}$ is only
$(16L-8)4^{L}$.

$L$ | $\bar{T}_{c,1}$ | $X_{c,1}$ | $c^{(1)}$ |  | $\bar{T}_{c,2}$ | $X_{c,2}$ | $c^{(2)}$
---|---|---|---|---|---|---|---
3 | 0.9630812528 | 0.1153561805 | 0.9439481004 |  | 0.9831360600 | 0.1951137102 | 0.9400187288
4 | 0.9581778321 | 0.1122361525 | 0.9019637901 |  | 0.9689256192 | 0.1810700192 | 0.8996045252
5 | 0.9546313471 | 0.1090076877 | 0.8651449784 |  | 0.9611962232 | 0.1708572156 | 0.8639745823
6 | 0.9530663416 | 0.1071998513 | 0.8438815963 |  | 0.9575263053 | 0.1648349316 | 0.8433113197
7 | 0.9524498263 | 0.1063326346 | 0.8322698047 |  | 0.9557272984 | 0.1612950160 | 0.8319934887
8 | 0.9522229606 | 0.1059533630 | 0.8256760684 |  | 0.9547690621 | 0.1590849673 | 0.8255533747
9 | 0.9521529738 | 0.1058159572 | 0.8216678338 |  | 0.9542109914 | 0.1576018062 | 0.8216323836
10 | 0.9521470966 | 0.1058006641 | 0.8190595772 |  | 0.9538600720 | 0.1565417468 | 0.8190774634
11 | 0.9521675217 | 0.1058493262 | 0.8172623591 |  | 0.9536256433 | 0.1557459461 | 0.8173147797
12 | 0.9521981542 | 0.1059325169 | 0.8159665891 |  | 0.9534614919 | 0.1551257762 | 0.8160423453
13 | 0.9522318614 | 0.1060345131 | 0.8149983146 |  | 0.9533421966 | 0.1546283056 | 0.8150903408
14 | 0.9522654200 | 0.1061464695 | 0.8142535202 |  | 0.9532528474 | 0.1542199764 | 0.8143571526
15 | 0.9522974007 | 0.1062632003 | 0.8136666597 |  | 0.9531842311 | 0.1538784508 | 0.8137787074
Table 1: Finite-size estimates of the pseudocritical temperature
$\bar{T}_{c,k}$ and scaling dimension $X_{c,k}$ from the condition
$X_{k}(L-1,\bar{T}_{c,k})=X_{k}(L+1,\bar{T}_{c,k})$ for $k=1$ (left) and $k=2$
(right). The central charges $c^{(k)}$, calculated at $\bar{T}_{c,k}(L)$, are
also shown.

If $\Lambda_{L}^{(k)}$ denotes the eigenvalues of $\tilde{M}_{L}$, for given
values of $\bar{\mu}$ and $\bar{T}$, with
$\Lambda_{L}^{(0)}>|\Lambda_{L}^{(1)}|>|\Lambda_{L}^{(2)}|>\ldots$, then, the
(dimensionless) negative free energy density is

$f(L)=\frac{\zeta}{L}\ln\Lambda_{L}^{(0)},$ (3)

where $\zeta$ is a geometric factor given by the thickness of the added row
when its circumference is $L$, which in our case is $\zeta=2/\sqrt{3}$. The
correlation length $\xi_{k}(L)$ associated with the $k$th subleading
eigenvalue is given by

$\xi_{k}^{-1}(L)=\zeta\ln\frac{\Lambda_{L}^{(0)}}{|\Lambda_{L}^{(k)}|}.$ (4)

From these correlation lengths, we may defined the scaled gaps
$X_{k}(L)=L/[2\pi\xi_{k}(L)]$. Note that all these quantities obviously depend
also on $\bar{\mu}$ and $\bar{T}$. Close to a critical point ($\bar{T}_{c}$ or
$\bar{\mu}_{c}$), from finite-size scaling Fisher (1971), we expect that

$X_{k}(L;\bar{T},\bar{\mu})=X_{k}^{*}+a_{t}tL^{y_{t}}+a_{1}L^{y_{1}}+a_{2}L^{y_{2}}+\ldots,$
(5)

where $y_{t}>0$ is the temperature exponent, $y_{i}<0$ ($i=1,2,\ldots$) are
irrelevant exponents, and $t=\bar{T}-\bar{T}_{c}$ if $\bar{\mu}$ is kept fixed
or $t=\bar{\mu}-\bar{\mu}_{c}$ when $\bar{T}$ is the fixed parameter.

In conformally invariant systems, the critical free energy density is expected
to follow

$f(L)=f_{\infty}+\frac{\pi c}{6L^{2}}+p_{1}L^{z_{1}}+\ldots,$ (6)

where $c$ is the conformal anomaly Blöte, Cardy, and Nightingale (1986);
*affleck86 and $z_{1}$ is a negative exponent.

###  III.2 The full-occupancy limit

Let us start analyzing the limiting case $\bar{\mu}\rightarrow\infty$, where
all lattice sites are occupied by molecules. Here, we can work with smaller
transfer matrices, as noticed above, what allow us to investigate strip widths
up to $L=16$. As recently demonstrated by some of us, in such a limit the
system displays only the DF and HDL phases, separated by a critical point in
the class of the two-dimensional (2D) 3-state Potts model Furlan  _et al._
(2019).

If a critical point exists at some critical temperature $\bar{T}_{c}^{*}$,
from Eq. 5, one may expect that curves of the scaled gaps
$X_{k}(L-\ell;\bar{T})$ and $X_{k}(L+\ell;\bar{T})$ versus $\bar{T}$ (e.g.,
with $\ell=1$) cross at pseudocritical temperatures $\bar{T}_{c,k}(L)$, which
shall converge to $\bar{T}_{c}^{*}$ as $L\rightarrow\infty$. Therefore,
$\bar{T}_{c,k}(L)$ can be determined from the condition
$X_{k}(L-1,\bar{T}_{c,k})=X_{k}(L+1,\bar{T}_{c,k})\equiv X_{c,k}(L)$, from
which a pseudocritical estimate [$X_{c,k}(L)$] of the asymptotic scaling
dimension $X_{k}^{*}$ is also obtained. The values of these quantities for
$k=1$ and $k=2$ are shown in Tab. 1. The first thing to note there is that
$\bar{T}_{c,1}$ has a non-monotonic behavior, initially decreasing and then
increasing with $L$. This prevents any reliable extrapolation of this
temperature. The convergence of $\bar{T}_{c,2}$, on the other hand, is
monotonic and, so, we will focus on it to estimate $\bar{T}_{c}^{*}$.

In general, from finite-size scaling, one expects that

$\bar{T}_{c,k}(L)=\bar{T}_{c}^{*}+b_{1}L^{-\omega_{1}}+b_{2}L^{-\omega_{2}}+\ldots,$
(7)

whose exponents are related to those in Eq. 5 as $\omega_{1}=y_{t}-y_{1}$ and
$\omega_{2}=y_{t}-y_{2}$. Then, by considering $b_{i}=0$ for $i\geqslant 2$
and performing three-point (3-pt) extrapolations of the values of
$\bar{T}_{c,2}$ in Tab. 1 (with $\bar{T}_{c}^{*}$, $b_{1}$ and $\omega_{1}$ as
unknowns in Eq. 7), we find the temperatures displayed in Tab. 2. Note that
for the largest $L$’s they are varying at the fifth decimal place. A second
round of 3-pt extrapolations [now of the extrapolated values in Tab. 2] yields
$0.952806\leqslant\bar{T}_{c}^{*}\leqslant 0.952818$, for the four largest
sets of widths, without any clear tendency to increase or decrease. Hence,
this gives $\bar{T}_{c}^{*}=0.952812(6)$. This value differs by only $0.02$%
from the one estimated from MC simulations in Ref. Furlan  _et al._ (2019):
$\bar{T}_{c}^{*}\approx 0.9526$.

The resulting exponents $\omega_{1}$ from the first 3-pt extrapolations are
also displayed in Tab. 2. Three-point extrapolations of them, assuming again a
power-law correction, return values close (and converging toward) to
$\omega_{1}=2$ for the largest sets of widths. We remark that, according to
Queiroz de Queiroz (2000), the irrelevant exponents in Eq. 5 are given by
$y_{1}=-4/5$ and $y_{2}=-8/5$ for the 3-state Potts model, whose critical
thermal exponent is $y_{t}=6/5$. This means that $\omega_{1}=2$ and
$\omega_{2}=14/5=2.8$ in this case. By assuming that $\omega_{1}=2$ and
$b_{i}=0$, for $i\geqslant 3$ in Eq. 7, we estimate the exponents $\omega_{2}$
for different $L$’s from 4-pt extrapolations; and subsequent 3-pt
extrapolations of them yield values close to $14/5$ for the largest $L$’s. So,
the correction exponents of the ALG model in the full lattice regime are
consistent with the ones for the 3-state Potts model.

$L$ | $\bar{T}_{c}^{*}$ | $\omega_{1}$ | $c$
---|---|---|---
4 | 0.9397502 | 1.3793263 | –
5 | 0.9516285 | 2.6536683 | 0.7857046
6 | 0.9529237 | 3.2158880 | 0.8086578
7 | 0.9530776 | 3.3614896 | 0.8119752
8 | 0.9530316 | 3.2893957 | 0.8121525
9 | 0.9529681 | 3.1489286 | 0.8117281
10 | 0.9529211 | 3.0137856 | 0.8112940
11 | 0.9528904 | 2.9040695 | 0.8109670
12 | 0.9528705 | 2.8173407 | 0.8107358
13 | 0.9528570 | 2.7468393 | 0.8105667
14 | 0.9528474 | 2.6875514 | 0.8104331
Table 2: Extrapolated temperatures, $\bar{T}_{c}^{*}$, and corresponding
correction exponents, $\omega_{1}$, from 3-pt extrapolations of the values of
$\bar{T}_{c,2}$ displayed in Tab. 1, considering Eq. 7 and sets of sizes
($L-1,L,L+1$). The resulting central charges, $c$, from analogous
extrapolations of the values of $c^{(2)}$ in Tab. 1 are also shown.

The finite-size amplitudes $X_{c,1}$ and $X_{c,2}$ in Tab. 1, on the other
hand, do not agree with the ones expected for the 2D 3-state Potts class. In
fact, for conformally invariant models in this class, the asymptotic values of
these quantities are $X_{1}^{*}=2-y_{h}=2/15$ and $X_{2}^{*}=2-y_{t}=4/5$.
However, our values of $X_{c,1}$ for the largest $L$’s are $\sim 20$% smaller
than $2/15$, though they display a non-monotonic convergence and are
increasing for large $L$’s. So, it may be the case that for very large strip
widths they will attain the expected Potts value. For $X_{c,2}$, our values
are smaller than $4/5$ by a factor $\sim 5$ and are still decreasing, without
any indication of changing in this behavior. Actually, the values in Tab. 1
extrapolate to $X_{2}^{*}\approx 0.15$, suggesting that $y_{t}\approx 1.85$
and then $\nu=1/y_{t}\approx 0.54$. Thereby, we can be either dealing with a
system displaying very slow convergences in these quantities (what seems
strange in the case of $X_{c,2}$…) or the ALG is not a conformally invariant
model or it belongs to a different universality class.

The central charge is somewhat consistent with the 3-state Potts class, for
which $c=4/5$. In fact, the values of $c$, calculated as
$c(L)=3[f(L-1)-f(L+1)](L-1)^{2}(L+1)^{2}/2\pi L$ at the pseudocritical
temperatures $\bar{T}_{c,1}$ and $\bar{T}_{c,2}$, are shown in Tab. 1; and
they approximate $4/5$ as the system size increases. However, 3-pt
extrapolations of such estimates (assuming power-law corrections) return
$c\gtrsim 0.81$ [see Tab. 2]. Further 3-pt extrapolations of the extrapolated
values in Tab. 2, considering the largest $L$’s available, do not improve
this, once they yield $c\approx 0.809$, which differs by $\sim 1$% from the
Potts value.

$L$ | $f$ | $X_{1}$ | $X_{2}$ | $c$
---|---|---|---|---
2 | 5.029729 | 0.112107 | 0.180682 | –
3 | 4.960323 | 0.109445 | 0.168819 | 0.944204
4 | 4.937032 | 0.107661 | 0.161966 | 0.902030
5 | 4.926737 | 0.106952 | 0.158348 | 0.865111
6 | 4.921304 | 0.106774 | 0.156322 | 0.843869
7 | 4.918080 | 0.106823 | 0.155076 | 0.832298
8 | 4.916008 | 0.106965 | 0.154240 | 0.825738
9 | 4.914595 | 0.107143 | 0.153641 | 0.821754
10 | 4.913587 | 0.107334 | 0.153189 | 0.819162
11 | 4.912844 | 0.107525 | 0.152836 | 0.817377
12 | 4.912280 | 0.107713 | 0.152551 | 0.816090
13 | 4.911841 | 0.107895 | 0.152317 | 0.815052
14 | 4.911494 | 0.108069 | 0.152121 | 0.814293
15 | 4.911213 | 0.108237 | 0.151953 | 0.813804
16 | 4.910984 | 0.108397 | 0.151809 | –
Table 3: Critical free energies, $f$, and scaled gaps $X_{1}$ and $X_{2}$
calculated at the critical temperature $\bar{T}_{c}^{*}=0.952812$. The central
charges, $c$, estimated from these free energies are also shown.

Table 3 shows the scaled gaps $X_{1}$ and $X_{2}$ calculated at the critical
point $\bar{T}_{c}^{*}=0.952812$. However, these quantities are very similar
to those in Tab. 1, obtained at the pseudocritical temperatures. The critical
free energies at $\bar{T}_{c}^{*}=0.952812$, as well as the central charges
obtained from them (calculated as indicated above) are also depicted in Tab.
3. Once again, the values of $c$ present an appreciable difference from
$c=4/5$ and extrapolations of them have deviations similar to those observed
above.

###  III.3 The DF-HDL transition line

Figure 5: (Color online) a) Pseudocritical lines separating the DF and HDL
phases, obtained from $X_{2}(L-1,\bar{T}_{c}^{*})=X_{2}(L+1,\bar{T}_{c}^{*})$
for fixed values of $\bar{\mu}$, for the indicated $L$’s. The outcomes from
3-pt extrapolations of these lines considering the sets $(5,7,9)$ [Extrap. 1]
and $(7,9,11)$ [Extrap. 2] are also shown. b) Detail of the region where the
DF-HDL, DF-LDL and LDL-HDL transition lines, for $L=7$, meet. The green square
and the red circle indicate the two estimates for the critical-end-point.

Employing the same analysis from the previous subsection, but for finite
$\bar{\mu}$, we find that the critical DF-HDL transition extends to chemical
potentials down to $\bar{\mu}\sim 2$, giving rise to a DF-HDL transition line,
as found in the mean-field solutions of the model on Husimi lattices Oliveira,
Stilck, and Barbosa (2010); Furlan  _et al._ (2019). Given the non-monotonic
convergence observed above in the critical points estimated from $X_{1}$, we
will focus here only on those calculated from $X_{2}$. The pseudocritical
transition lines obtained from the condition
$X_{2}(L-1,\bar{T}_{c}^{*})=X_{2}(L+1,\bar{T}_{c}^{*})$, for fixed values of
$\bar{\mu}$, are shown in Fig. 5(a), for $L\leqslant 11$. (Here, it is quite
hard to go beyond $L=12$, once the TMs have $4^{L}\times 4^{L}$ terms.
Actually, most of our results will be limited to $L\leqslant 10$.) It is
important to notice that such lines never cross each other, indicating that
the DF-HDL transition is always continuous.

Three-point extrapolations of the pseudocritical lines — assuming the finite-
size correction of Eq. 7 for fixed $\bar{\mu}$, with $a_{i}=0$ for $i>1$, and
considering the sets of widths $(5,7,9)$ and $(7,9,11)$ — yield lines whose
values of $\bar{T}_{c}^{*}$ always differ by less than $0.2$%. This
corresponds to differences (and error bars) at most at the third decimal
place, which is accurate enough for our purposes.

Interestingly, the pseudocritical lines have a corner at $\bar{\mu}\approx 2$,
which may be seen as their end-point. An example of this is shown in Fig.
5(b), in the line for $L=7$. As it will be demonstrated below, the critical
DF-HDL line ends at a critical end-point (CEP). The coordinates
($\bar{T}_{cep,1},\bar{\mu}_{cep,1}$) of the CEP, assuming that it is located
at the corner, are depicted in Tab. 4. Since the DF-HDL transition lines
continue below these corners and finally meet the DF-LDL and LDL-HDL
transition lines [see Fig. 5(b)], one obtains another set of estimates
($\bar{T}_{cep,2},\bar{\mu}_{cep,2}$) for the CEP. As shows Tab. 4, in this
last case the temperatures oscillate around $\bar{T}_{cep,2}=0.81266$
(considering the three largest $L$’s). The temperatures of the corners seem to
converge toward this value as $L$ increases; though a 3-pt extrapolation of
them (following Eq. 7, considering the largest $L$’s) returns a slightly
larger value: $\bar{T}_{cep,1}=0.81681$. This deviation, of $\sim 0.5$%, is
certainly due to the small widths considered here, so that it is reasonable to
regard $\bar{T}_{cep}=0.815(3)$ as our estimate for the CEP temperature. This
value is close, but a bit smaller than $\bar{T}_{c/tc}=0.825$, which was found
in previous MC simulations of the ALG model and reported as a critical point
in Balladares, Henriques, and Barbosa (2007) and as a tricritical point in
Furlan, Fiore, and Barbosa (2015). Our estimate is also smaller than the value
found for the triple point ($\bar{T}_{tp}=0.835$) in the Husimi lattice
solution of the model Oliveira, Stilck, and Barbosa (2010).

$L$ | $\bar{T}_{cep,1}$ | $\bar{\mu}_{cep,1}$ | $\bar{T}_{cep,2}$ | $\bar{\mu}_{cep,2}$
---|---|---|---|---
3 | 0.884217 | 2.264032 | 0.801210 | 2.010048
5 | 0.841799 | 2.018722 | 0.812741 | 1.982581
7 | 0.828528 | 1.981947 | 0.812563 | 1.972774
9 | 0.823467 | 1.970668 | 0.812692 | 1.968213
Table 4: Estimates for the critical-end-point from the corner in the DF-HDL
transition lines (leftmost values) and from the point where such lines meet
the other transition lines (rightmost values).

The 3-pt extrapolations of the values of $\bar{\mu}_{cep,1}(L)$ and
$\bar{\mu}_{cep,2}(L)$ for the largest $L$’s, assuming power-law corrections
similarly to Eq. 7, yield $\bar{\mu}_{cep,1}=1.9589$ and
$\bar{\mu}_{cep,2}=1.9605$, respectively, indicating that
$\bar{\mu}_{cep}=1.9597(8)$. This result agrees with the mean-field value
found in Ref. Oliveira, Stilck, and Barbosa (2010) for the triple point
($\bar{\mu}_{tp}=1.959$) and differs by $\approx 3$% from the critical point
reported in Ref. Balladares, Henriques, and Barbosa (2007)
($\bar{\mu}_{c}=2.02$).

Figure 6: (Color online) Central charge $c$ versus $\bar{\mu}$ calculated at
the pseudocritical transition lines from Fig. 5(a), for the indicated widths.

The central charges calculated along the pseudocritical transition lines from
Fig. 5(a) are shown in Fig. 6 as a function $\bar{\mu}$. We may note that
$c(L)$ is almost independent of the chemical potential for
$\bar{\mu}_{c}\gtrsim 10$ (similarly to $\bar{T}_{c}$), having values very
close to those of $c^{(2)}$ in Tab. 1 for $\bar{\mu}\rightarrow\infty$. When
$\bar{\mu}$ approximates $\bar{\mu}_{cep}$, however, one sees that $c$ passes
to increase and to display stronger finite-size corrections. The large $L$ is
the close the approximately constant part of the lines get to
$\bar{\mu}_{cep}$ before deviating upward. This strongly suggests that for
very large $L$’s one must have $c\approx 0.8$ in the entire DF-HDL line,
indicating that it belongs to the 3-state Potts class.

###  III.4 The LDL-HDL transition line

The LDL-HDL transition is the only one for which all previous studies of the
ALG model agree; the phases are separated by a coexistence line, which exists
in a narrow chemical potential interval near $\bar{\mu}=2$ Balladares,
Henriques, and Barbosa (2007); Oliveira, Stilck, and Barbosa (2010); Furlan,
Fiore, and Barbosa (2015). Therefore, in this case, it is convenient to
determine the transition lines by varying $\bar{\mu}$ for fixed values of
$\bar{T}$. For the strip widths analyzed here, crossings in the curves of
$X_{k}(L;\bar{\mu})$ versus $\bar{\mu}$ for different $L$’s are only observed
for $k\geqslant 3$. The curves of $X_{1}(L;\bar{\mu})$ and
$X_{2}(L;\bar{\mu})$ only exhibit a maximum in the coexistence region. The
transition lines obtained from the condition
$X_{3}(L-1,\bar{\mu}_{c})=X_{3}(L+1,\bar{\mu}_{c})$ are depicted in Fig. 7(a).
The lines for the smallest widths present maxima in the region close to the
CEP, which are absent in those for the largest $L$’s. This turns the
extrapolation of these lines unreliable in this region, which is the most
important part of such curves, since they are always very close to
$\bar{\mu}=2$ for low temperatures.

Figure 7: (Color online) a) Estimates of the LDL-HDL coexistence line from the
condition $X_{3}(L-1,\bar{\mu}_{c})=X_{3}(L+1,\bar{\mu}_{c})$ for fixed values
of $\bar{T}$, for the indicated $L$’s. b) Inverse correlation length
$\xi_{6}^{-1}$ versus $\bar{\mu}$ for $\bar{T}=0.75$ and the indicated $L$’s.
c) LDL-HDL coexistence lines obtained from the minima in the curves in (b) and
similarly ones for other temperatures. The extrapolated curve [Extrap.], from
a 3-pt extrapolation following Eq. 8, is also shown.

Motivated by this and also to confirm the first-order nature of the LDL-HDL
transition, we obtain an alternative estimate for the LDL-HDL transition line
through the correlation length $\xi_{k}^{-1}(L)$, defined in Eq. 4, for $k=6$.
We remark that, given the four-fold (three-fold) degeneracy of the ground
state of the LDL (HDL) phase, a total of seven possible phases (divided into
two sets with three and four equivalent ones) coexist at the LDL-HDL
transition line. This means that the seven largest eigenvalues (from
$\Lambda_{0}$ to $\Lambda_{6}$) of the TMs shall be degenerated in the
thermodynamic limit ($L\rightarrow\infty$). Hence, curves of $\xi_{6}^{-1}(L)$
vs $\bar{\mu}$, for a given $\bar{T}$, are expected to display a minimum at
the coexistence Rikvold  _et al._ (1983); Bartelt, Einstein, and Roelofs
(1986); Xavier  _et al._ (1998); Jung and Kim (2017). As exemplified in Fig.
7(b) for $\bar{T}=0.75$, such behavior is indeed found in these curves,
confirming that this is a discontinuous transition. The coexistence lines
estimated from the loci of these minima are depicted in Fig. 7(c). It is
noteworthy that $\bar{\mu}_{c}(\bar{T})$ decreases with $L$ in Fig. 7(a),
while for these last estimates an increasing behavior is observed in Fig.
7(c).

In order to extrapolate these lines, we recall that finite-size scaling
predicts exponential corrections at the coexistence Fisher (1971), such that

$\bar{\mu}_{c}(L)=\bar{\mu}_{c}^{*}+ae^{-bL}+\ldots,$ (8)

for a given $\bar{T}$. For a 3-pt extrapolation this gives
$\bar{\mu}_{c}^{*}=\frac{\bar{\mu}_{c}(L+1)\bar{\mu}_{c}(L-1)-\bar{\mu}_{c}(L)^{2}}{\bar{\mu}_{c}(L+1)+\bar{\mu}_{c}(L-1)-2\bar{\mu}_{c}(L)}$.
The resulting curve from such a extrapolation, for the three largest $L$’s, is
shown in Fig. 7(c).

###  III.5 The DF-LDL transition line

The transition lines separating the disordered and LDL phases were estimated
following the same procedures from the previous subsection, but varying
$\bar{T}$ for fixed values of $\bar{\mu}$. The curves obtained from the
condition $X_{1}(L-1,\bar{T}_{c})=X_{1}(L+1,\bar{T}_{c})$ are depicted in Fig.
8(a), where one sees that $\bar{T}_{c}$, for a given $\bar{\mu}$ decreases
with $L$. Significantly, these lines never cross each other, indicating that
no tricritical point exists in the DF-LDL transition, in contrast with the
result from MC simulations reported in Ref. Furlan, Fiore, and Barbosa (2015).
Since one knows that this transition is discontinuous at $\bar{T}=0$, and for
low $\bar{T}$ as well, the absence of the tricritical point implies that the
entire transition line is of first-order nature, in agreement with the Husimi
lattice solution of the model Oliveira, Stilck, and Barbosa (2010).

Figure 8: (Color online) Finite-size estimates of the DF-LDL coexistence line:
a) from the condition $X_{1}(L-1,\bar{T}_{c})=X_{1}(L+1,\bar{T}_{c})$, for
fixed values of $\bar{\mu}$; and b) from the minima in the inverse correlation
length $\xi_{4}^{-1}$. The respective extrapolated curves [Extrap.], from a
3-pt extrapolation following Eq. 8, are shown in both panels.

Additional evidence of this is obtained from the curves of $\xi_{4}^{-1}$
versus $\bar{T}$, for a given $\bar{\mu}$, which have minima hallmarking the
coexistence. A behavior analogous to that seen in Fig. 7(b) for $\xi_{6}^{-1}$
at the LDL-HDL transition. Note that in the case of the DF phase (one ground
state) coexisting with the LDL phase (which has four equivalent ground
states), minima are indeed expected in $\xi_{4}^{-1}$. From such minima, we
estimate the coexistence lines displayed in Fig. 8(b). Interestingly, in this
case the values of $\bar{T}_{c}(\bar{\mu})$ increases with $L$, in opposition
to the behavior seen in Fig. 8(a).

Figures 8(a) and 8(b) present also the asymptotic curves, obtained from 3-pt
extrapolations assuming exponential finite-size corrections in
$\bar{T}_{c}(L)$ similarly to Eq. 8 [i.e.,
$\bar{T}_{c}(L)\simeq\bar{T}_{c}^{*}+ae^{-bL}$]. The difference between the
values of $\bar{T}_{c}^{*}(\bar{\mu})$ along these lines is always smaller
than $0.5$% and in the main part of them it is $\lesssim 0.01$%. The turning
point in the coexistence line (i.e., the point where $\bar{T}_{c}^{*}$ is
maximum) is located at $\bar{T}_{c}^{*}\approx 0.983$ and
$\bar{\mu}_{c}^{*}\approx 0.31$.

The complete phase diagram, presenting all the transition lines found in our
MT calculations, will be presented in Sec. V, where they are compared with the
results from the MC simulations discussed in the next section.

##  IV Monte Carlo simulations

###  IV.1 Simulation details and quantities of interest

The thermodynamic properties of the symmetric ALG model were also obtained via
extensive Monte Carlo (MC) simulation on the triangular lattice with periodic
boundary conditions. We study system sizes $L=32,48,64,96,128$ and $192$ using
the standard Metropolis algorithm in the grand canonical ensemble.

Following references Balladares, Henriques, and Barbosa (2007); Furlan, Fiore,
and Barbosa (2015), simulations with fixed chemical potential are used to
locate the DF-HDL line while simulations along isotherms in the
$\bar{\mu}-\bar{T}$ plane are used to locate the DF-LDL and LDL-HDL
coexistence lines. We study $\bar{\mu}=1.95,2.3,2.4,2.5,2.6,5.0$ and $10$ in
the HDL region and use a temperature interval $0.50\leq\bar{T}\leq 0.97$ in
the LDL region.

Studies of the HDL region are performed in two steps. First, broad ranges of
temperature are studied in order to identify the region of peaks in the
response functions. In this step $50$ independent realizations are used to
compute averages. After identifying the region of interest for each system
size, additional simulations are performed in a reduced range of temperatures.
The number of points around the maximum depends on the system size; for
$L=64,128$ we use $60$ points while for larger system sizes $30$ points are
sufficient. In this step 100 independent replicas are used. In all simulations
we use $10^{6}$ MC steps for equilibration, followed by $10^{5}-2\times
10^{6}$ MC steps for data production.

The transition between the HDL and the DF phases can be analyzed using the
order parameters, Szortyka  _et al._ (2009); Furlan  _et al._ (2019)

$\theta(\bar{T},\bar{\mu};L)\equiv\frac{3}{2}\left[\frac{\max(n_{1},n_{2},n_{3})}{N}-\frac{1}{3}\right],\\\
$ (9)

and

$Q(\bar{T},\bar{\mu};L)\equiv\frac{1}{N}\sqrt{n^{2}_{1}+n^{2}_{2}+n^{2}_{3}-n_{1}n_{2}-n_{1}n_{3}-n_{2}n_{3}},\\\
$ (10)

where $n_{k}$ denotes the number of molecules in orientational state $k$ and
$N=\sum_{k=1}^{3}n_{k}$ is the total number of molecules. In the DF phase, one
has $n_{1}\approx n_{2}\approx n_{3}$, so that $\theta\rightarrow 0$ and
$Q\rightarrow 0$. On the other hand, $\theta\rightarrow 1$ and $Q\rightarrow
1$ in the HDL phase, where, e.g., $n_{1}\approx N$ and $n_{2}\approx
n_{3}\approx 0$.

To analyze the transitions we also use the susceptibilities related to the
order parameters:

$\chi_{X}\left(\bar{T};L\right)=\frac{L^{2}}{\bar{T}}\left[\left\langle
X^{2}\right\rangle-\left\langle X\right\rangle^{2}\right]$ (11)

with $X=\theta$ or $Q$; the specific heat at constant volume Allen and
Tildesley (1989)

$c_{V}\left(\bar{T};L\right)=\frac{1}{V\bar{T}^{2}}\left[\left(\left\langle
E^{2}\right\rangle-\left\langle
E\right\rangle^{2}\right)-\frac{\left(\left\langle
EN\right\rangle-\left\langle E\right\rangle\left\langle
N\right\rangle\right)^{2}}{\left\langle N^{2}\right\rangle-\left\langle
N\right\rangle^{2}}\right];$ (12)

and the isothermal compressibility

$\kappa\left(\bar{T};L\right)=\frac{V}{k_{B}\bar{T}}\left(\frac{\left\langle
N^{2}\right\rangle-\left\langle N\right\rangle^{2}}{\left\langle
N\right\rangle^{2}}\right).$ (13)

In Eqs. 11-13, $\left\langle\cdot\right\rangle$ denotes grand canonical
averages.

###  IV.2 Discontinuous transitions

In this subsection we present the results for DF-LDL and LDL-HDL transitions.
As reported in the Refs. Balladares, Henriques, and Barbosa (2007); Furlan,
Fiore, and Barbosa (2015) the density plays the role of an order parameter for
the referred transitions. Figures 9(a) and (b) show the density as a function
of chemical potential for the DF-LDL transition and the LDL-HDL transition,
respectively. The density histograms for each transition are shown in the
insets of Figs. 9(a) and (b). In both cases, the histograms are bimodal
indicating a discontinuous transition and confirming the results from the
transfer matrix (TM) calculations (see Subsec. III.4 and III.5). For each
temperature, we plotted the histogram for the chemical potential that was
closest to the coexistence value within our resolution
($\Delta\bar{\mu}=10^{-4}$ or larger). We take the coexistence chemical
potential as that for which the heights of the two peaks are equal.

Figure 9: Density versus chemical potential for distinct temperatures at DF-
LDL (a) and LDL-HDL (b) transitions. Insets show histograms of the density for
DF-LDL and LDL-HDL near coexistence. Different colors represent temperatures
as indicated. System size is $L=64$ in (a) and $L=128$ in (b).

The coexistence chemical potential is estimated as follows: First we identify
the values of the maxima of the probability distributions,
$\max\left[P\left(\rho_{DF}\right)\right]$ and
$\max\left[P\left(\rho_{LDL}\right)\right]$, for each chemical potential
studied (see blue crosses and red pluses in Fig. 10). Subsequently, we plot
the $\max\left[P\left(\rho_{LDL,DF}\right)\right]$ as function of chemical
potential as shown in Fig. 10 for $\bar{T}=0.55$. The crossing point of the
curves is estimated through a polynomial fit in the vicinity of the
intersection, using the the reciprocals of the variance ($1/\sigma_{i}^{2}$)
of the averages computed over the replicas as weights for the polynomial fit
of the points. In most cases, 7 points around the intersection region are
enough to obtain a satisfactory region of crossing. The order of the
polynomial is chosen as the lowest order that yields random behavior of
residuals; a 3rd-order polynomial fit was sufficient in most cases.
Subsequently, the error of the fit was computed through the root mean square
error (RMSE) over the points used in the fit. Finally, the uncertainty in the
crossing-point estimate is computed as the square root of the RMSE resulting
from each polynomial. This procedure yields the dots in the phase diagram of
Fig. 15.

Figure 10: Main panel: thin black lines represent the probability distribution
of the density $\rho$ for distinct chemical potentials and $\bar{T}=0.55$.
Blue(red) crosses (pluses) indicate the location of the maximum of probability
distributions in the DF (LDL) phase. Inset: maxima of the probability
distributions for the disordered phase (blue crosses) and the LDL phase (red
pluses) versus chemical potential $\bar{\mu}$. The vertical line indicates the
crossing point. System size $L=64$.

In the region around $\bar{\mu}=0$, $\bar{T}>0.95$ the densities $\rho_{LDL}$
and $\rho_{DF}$ are very similar and the coexistence temperature and chemical
potential cannot been obtained using the density histogram (see Fig. 15(b)).

For the LDL-HDL transition, the data for the density are very noisy for
$\bar{T}\leq 0.75$ (see Fig. 9(b)) and it is not possible to use the procedure
described above to determine the coexistence chemical potential. For these
cases, we adopted the following procedure: for each point on the curve
$\max\left[P\left(\rho_{i}\right)\right]\times\bar{\mu}$ (see inset of Fig.
10) we sum a Gaussian random variable proportional to the error bar of that
point. We then fit the curves and estimate the crossing point. An average is
computed over $10^{4}$ realizations of this procedure and the error is
propagated as described above. The values obtained using this strategy are the
triangles in Fig. 15 and are in good agreement with the TM results of Subsec.
III.4.

The highest temperature for which we observe LDL-HDL coexistence is
$\bar{T}=0.816$, $\bar{\mu}=1.9498(7)$, for the next temperature studied,
$\bar{T}=0.817$, we observe LDL-DF coexistence and then a continuous DF-HDL
transition. We therefore take this $\bar{T}$, $\bar{\mu}$ pair as an estimate
of the critical end point (CEP) at which the LDL-HDL line meets the DF-HDL and
the DF-LDL lines (see Subsec. III.3). This value, $\bar{T}_{cep}=0.816$,
$\bar{\mu}_{cep}=1.9498(7)$, is in excellent agreement with the estimate from
TM calculations, $\bar{T}_{cep,TM}=0.815(1)$, $\bar{\mu}_{cep,TM}=1.9597(8)$.
Our estimate for $\bar{T}_{cep,MC}$ and $\bar{\mu}_{cep,TM}$ is smaller than
the values reported in previous MC simulations Balladares, Henriques, and
Barbosa (2007); Furlan, Fiore, and Barbosa (2015) for $c_{2}$ and $tc_{2}$
(see Fig. 1).

###  IV.3 Continuous transitions

To determine the critical temperatures we use finite-size scaling (FSS)Allen
and Tildesley (1989) analysis of the susceptibilities
$\chi_{X}\left(\bar{T};L\right)$ of the order parameters $X=\theta$ and $Q$
(Eq. 11), the specific heat at constant volume $c_{V}\left(\bar{T};L\right)$
(Eq. 12), and isothermal compressibility (Eq. 13).

Figure 11 shows simulation results for the order parameters, $Q$ and $\theta$;
susceptibilities, $\chi_{Q}$ and $\chi_{\theta}$, specific heat $c_{V}$, and
isothermal compressibility $\kappa$ for chemical potential $\bar{\mu}=2.3$.
For clarity, error bars are not shown in the figure. The susceptibility,
specific heat and isothermal compressibility exhibit signatures of a
continuous transition. In all cases the peaks increase systematically with
system size. Similar results are obtained for $\bar{\mu}=2.4,2.5,2.6,5.0$, and
10. As the chemical potential is increased, the maxima shift to higher
temperatures, approaching the value
$\bar{T}^{*}_{c,\bar{\mu}\to\infty}=0.9526$ reported in Ref. Furlan  _et al._
(2019). This behavior confirms that the DF-HDL line is continuous, in
agreement with the TM calculations from Subsec. III.3.

Figure 11:  (a) Order parameters $Q$ and $\theta$, (b) susceptibilities
$\chi_{Q}$ and $\chi_{\theta}$, (c) specific heat $c_{V}$, and (d) isothermal
compresibility for $\bar{\mu}=2.3$. Blue, red, magenta, black, cyan, and
yellow correspond to system sizes $L=32,48,64,96,128$ and $192$. The order
parameter $Q$ and its susceptibility $\chi_{Q}$ are represented by solid
lines, $\theta$ and $\chi_{\theta}$ by symbols.

Using the results for $\chi_{Q},\chi_{\theta},c_{V}$ and $\kappa$ we estimate
the critical temperature $\bar{T}_{c}^{*}$. Initially, we estimate the size
dependent pseudocritical temperature $\bar{T}_{c}$ through a polynomial fit to
the data near the maximum, as described in the preceding subsection.

The pseudocritical temperatures $\bar{T}_{c}$ for $\bar{\mu}=2.3$ are plotted
versus $1/L$ in Fig. 12. Similar results were obtained for
$\bar{\mu}=2.4,2.5,2.6,5.0$ and $10$. The set of pseudocritical temperatures
appear to converge as $L\to\infty$ for all the studied values of $\bar{\mu}$.
The global estimate of $\bar{T}^{*}_{c}$ was obtained though a weighted
average with weight $1/\sigma^{2}$, where $\sigma$ represents the uncertainty
of each quantity. The final estimates of $\bar{T}^{*}_{c}$ are the squares in
the phase diagram Fig. 15.

Figure 12: Pseudocritical temperatures associated with, $c_{v}$, $\chi_{Q}$,
$\chi_{\theta}$ and $\kappa$, with symbols as indicated, versus $1/L$ for
$\bar{\mu}=2.3$. Dash lines are linear fits to the data and the filled black
dot indicates the global estimate for the critical temperature
$\bar{T}^{*}_{c}$.

###  IV.4 Critical exponents

For $\bar{\mu}=2.3$ and $10$, we estimate the critical- exponent ratios
$\beta/\nu$ and $\gamma/\nu$ from fits to the data for $Q(\bar{T}^{*}_{c})$,
$\theta(\bar{T}^{*}_{c})$, $\chi_{Q}(\bar{T}^{*}_{c})$, and
$\chi_{\theta}(\bar{T}^{*}_{c})$ versus system size (see Fig. 13). From the
linear fits, and including the effect of the uncertainty in $\bar{T}^{*}_{c}$,
we obtained the estimates in Table 5. The values for these ratios for the two-
dimensional three-state Potts model are
$\beta_{Potts}/\nu_{Potts}=2/15=0.1333\cdots$ and
$\gamma_{Potts}/\nu_{Potts}=26/15=1.7333\cdots$. In the worst cases, the
discrepancy is about $6\%$ for $\mu=2.3$ and $8\%$ for $\mu=10$ .

Figure 14 shows the linear fit for $\ln c_{v}(\bar{T}_{c})$ for $\mu=10.0$,
yielding the estimate $\alpha/\nu=0.44(1)$, i.e., a discrepancy with the ratio
for the two-dimensional, three-state Potts model $\alpha/\nu_{Potts}=0.4$ of
10%. For $\bar{\mu}=2.3$, the discrepancy in this ratio is $73\%$. We believe
this huge difference is likely due to finite-size corrections that are more
important close to the CEP ($\bar{\mu}_{cep}=1.9498$) where the DF-HDL line
encounters the LDL-HDL and DF-LDL lines.

Figure 13: (a) Order parameter $Q$($\theta$) versus system size for $\bar{\mu}=2.3$. The circles (squares) are the simulation data; dashed lines are least-squares linear fits of the data. Inset shows the residuals of the fits. (b) Analogous plots for the susceptibilities $\chi_{Q}$($\chi_{\theta}$).  Figure 14: Specific heat $c_{V}$ versus system size at critical temperature for $\bar{\mu}=10.0$ The circles are the simulation data; the dashed lines is a least-squares linear fit of the data. Inset shows the residuals of the fit. $\mu$ | $\beta_{Q}/\nu$ | $\beta_{\theta}/\nu$ | $\gamma_{Q}/\nu$ | $\gamma_{\theta}/\nu$ | $\alpha/\nu$
---|---|---|---|---|---
2.3 | 0.126(7) | 0.136(8) | 1.66(4) | 1.63(4) | 0.11(1)
10.0 | 0.126(7) | 0.122(7) | 1.75(5) | 1.75(4) | 0.44(1)
Table 5: Estimates of the ratios $\beta/\nu$, $\gamma/\nu$ and $\alpha/\nu$
for chemical potentials $\bar{\mu}=2.3$ and 10.

##  V Discussion

Figure 15(a) shows the phase diagram of the symmetric ALG model in the
$\bar{T}-\bar{\mu}$ plane, summarizing the results reported above. The
predictions of the TM analysis are in very good quantitative agreement with MC
simulation. In fact, the maximum difference between the estimates from each
method is smaller than 2% along the LDL-HDL coexistence line and, notably, is
at most 0.1% for the DF-HDL and DF-LDL lines. The natures of the transitions
are also the same in both cases. This remarkable agreement confirms that we
are accessing the correct thermodynamic behavior of the model.

Figure 15: Phase diagram for the symmetric ALG model in the
$\bar{\mu}\times\bar{T}$ plane (a) and the $\bar{T}\times\rho$ plane (b).
Symbols are results from Monte Carlo (MC) simulations. Dashed lines represent
discontinuous transitions and solid ones continuous transitions. Only in (a)
Lines are results obtained using Transfer Matrix (TM) calculations. The black
triangle marks the location of the critical end point (CEP) estimated using
TM, the black square marks the location of the critical point in the limit
$\bar{\mu}\to\infty$ obtained via MC simulation in Ref. Furlan  _et al._
(2019) and the black dot marks the density at $\bar{\mu}=1.95$ and $L=128$,
very close to our estimate of $\bar{\mu}_{cep}=1.9498(6)$. The green dotted
line is a conjecture as to how to critical line meets the CEP.

Our results demonstrate that the DF-LDL transition is always discontinuous,
without any critical or multicritical point. This is confirmed by the absence
of a crossing in the pseudocritical transition lines from the TM analysis, by
the bimodal histograms in MC simulations and also in the phase diagram of
temperature versus total density ($\bar{T}-\rho$), depicted in Fig. 15(b).
This diagram makes it clear that the densities are discontinuous at the DF-LDL
transition line up to temperatures very close to the turning point (in the
$\bar{\mu}-\bar{T}$ diagram), which is located at $\bar{T}\approx 0.983$. It
is, however, difficult to estimate the densities accurately very near this
point, since $\rho_{DF}$ becomes very similar to $\rho_{LDL}$. An analogous
behavior was found in the Husimi lattice solution of this model Oliveira,
Stilck, and Barbosa (2010) and indeed our ($\bar{T}-\rho$) diagram is very
similar to the one found there [see Fig. 4(a) of Ref. Oliveira, Stilck, and
Barbosa (2010)]. Importantly, these results provide compelling evidence that
the entire DF-LDL transition line is discontinuous, ruling out a critical
point at $\bar{T}\approx 0.55$ or a tricritical point at $\bar{T}\approx
0.65$, as claimed in previous MC analyses in references Balladares, Henriques,
and Barbosa (2007) and Furlan, Fiore, and Barbosa (2015), respectively.

Also in contrast with earlier simulations, we find a critical line separating
the DF and HDL phases. This line meets the LDL-HDL and DF-LDL coexistence
lines at a critical end point (CEP). This shows that the point where the LDL-
HDL coexistence line ends — reported in literature as critical Balladares,
Henriques, and Barbosa (2007), triple Oliveira, Stilck, and Barbosa (2010) and
tricritical Furlan, Fiore, and Barbosa (2015) — is actually a CEP. Therefore,
although this ALG model has an interesting phase behavior, beyond displaying
some thermodynamic anomalies, it is not a good model for liquid water, as
initially suggested, because it has neither a liquid-gas coexistence line
ending at a critical point, nor a LLCP. In light of the discussion in the
Introduction, it seems that absence of a LLCP is a general feature of lattice
models with attractive orientational interactions competing with “van der
Waals” repulsion between first neighbor molecules.

Although determining the universality class of the continuous transitions
displayed by models for water is key to establish a full connection with the
actual behavior of water, this has not been widely investigated, perhaps
because the very nature of these transitions are the subject of controversy. A
recent account of these critical properties for realistic models can be found
in Ref. Debenedetti, Sciortino, and Zerze1 (2020). For lattice models, some
examples include the debate in the literature on the critical exponents of the
liquid-liquid transition in the classic Bell-Lavis model Fiore  _et al._
(2009); Simenas, Ibenskas, and Tornau (2014), and our previous study of the
present ALG model in the fully occupied limit, where the DF-HDL transition was
found to be in the 3-state Potts class. The central charge obtained here,
deviating by $\sim 1$% from the Potts value, gives additional evidence of
this. The scaled gaps, on the other hand, present a large deviation from those
of the 3-state Potts class. The analysis of the central charge for finite
chemical potentials indicates that it converges to the same value along the
entire DF-HDL critical line. Nevertheless, finite-size corrections become very
strong as $\bar{\mu}$ approaches $\bar{\mu}_{CEP}=1.9597(8)$. This certainly
explains why some critical exponent ratios estimated in simulations for
$\bar{\mu}=2.3$ exhibit considerable deviations from the Potts values, while
for $\bar{\mu}=10$ the deviations are much smaller. Despite this, the overall
picture suggests that this order-disorder transition is always in the 3-state
Potts class, as is expected from its three-fold ($S_{3}$) symmetry breaking.

## Acknowledgments

Authors acknowledge financial support from CAPES Brazil and the National
Institute of Science and Technology for Complex Systems of Brazil.

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Stanley (2013) H. E. Stanley, ed.,  _Liquid Polymorphism: Advances in Chemical Physics_, Vol. 152 (Wiley, NY, 2013).
  * Speedy and Angell (1976) R. J. Speedy and C. A. Angell, “Isothermal compressibility of supercooled water and evidence of a singularity at $-45\degree$C,” J. Chem. Phys. 65, 851 (1976).
  * Sastry  _et al._ (1996) S. Sastry, P. G. Debenedetti, F. Sciortino, and H. E. Stanley, “Singularity-free interpretation of the thermodynamics of supercooled water,” Phys. Rev. E 53, 6144 (1996).
  * Poole  _et al._ (1992) P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, “Phase behavior of metastable water,” Nature 360, 324 (1992).
  * Mishima, Calvert, and Whalley (1985) O. Mishima, L. D. Calvert, and E. Whalley, “An apparently 1st-order transition between 2 amorphous phases of ice induced by pressure,” Nature 314, 76 (1985).
  * Lynden-Bell  _et al._ (2010) R. M. Lynden-Bell, S. C. Morris, J. D. Barrow, J. L. Finney, and C. L. Harper Jr., eds.,  _Water and life: The unique properties of H 2O_ (CRC Press, Boca Raton, FL, 2010).
  * Soper and Ricci (2000) A. K. Soper and M. A. Ricci, “Structures of high-density and low-density water,” Phys. Rev. Lett. 84, 2881 (2000).
  * Mishima and Stanley (1998) O. Mishima and H. E. Stanley, “Decompression-induced melting of ice IV and the liquid–liquid transition in water,” Nature 392, 164 (1998).
  * Woutersen  _et al._ (2018) S. Woutersen, B. Ensing, M. Hilbers, Z. Zhao, and C. A. Angell, “A liquid-liquid transition in supercooled aqueous solution related to the HDA-LDA transition,” Science 359, 1127 (2018).
  * Liu  _et al._ (2005) L. Liu, S.-H. Chen, A. Faraone, C.-W. Yen, and C.-Y. Mou, “Pressure dependence of fragile-to-strong transition and a possible second critical point in supercooled confined water,” Phys. Rev. Lett. 95, 117802 (2005).
  * Chen  _et al._ (2006) S.-H. Chen, F. Mallamace, C.-Y. Mou, M. Broccio, C. Corsaro, A. Faraone, and L. Liu, “Violation of the Stokes-Einstein relation in supercooled water,” Proc. Natl. Acad. Sci. USA 103, 12974 (2006).
  * Mallamace  _et al._ (2007) F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, D. Majolino, V. Venuti, L. Liu, C.-Y. Mou, and S.-H. Chen, “Evidence of the existence of the low-density liquid phase in supercooled, confined water,” Proc. Natl. Acad. Sci. USA 104, 424 (2007).
  * Mallamace (2009) F. Mallamace, “The liquid water polymorphism,” Proc. Natl. Acad. Sci. USA 106, 15097 (2009).
  * Soper (2019) A. K. Soper, “Is water one liquid or two?” J. Chem. Phys. 150, 234503 (2019).
  * Palmer  _et al._ (2018) J. C. Palmer, P. H. Poole, F. Sciortino, and P. G. Debenedetti, “Advances in computational studies of the liquid–liquid transition in water and water-like models,” Chem. Rev. 118, 9129 (2018).
  * Bell and Lavis (1970a) G. M. Bell and D. A. Lavis, “Two-dimensional bonded lattice fluids. I. Interstitial model,” J. Phys. A 3, 427 (1970a).
  * Bell and Lavis (1970b) G. M. Bell and D. A. Lavis, “Two-dimensional bonded lattice fluids. II. Orientable molecule model,” J. Phys. A 3, 568 (1970b).
  * Huckaby and Hanna (1987) D. A. Huckaby and R. S. Hanna, “A two-dimensional lattice gas model for water,” J. Phys. A 20, 5311 (1987).
  * Patrykiejew, Pizio, and Sokołowski (1999) A. Patrykiejew, O. Pizio, and S. Sokołowski, “Novel phase behavior in a two-dimensional network-forming lattice fluid,” Phys. Rev. Lett. 83, 3442 (1999).
  * Buzano  _et al._ (2004) C. Buzano, E. D. Stefanis, A. Pelizzola, and M. Pretti, “Two-dimensional lattice-fluid model with waterlike anomalies,” Phys. Rev. E 69, 061502 (2004).
  * Henriques and Barbosa (2005) V. B. Henriques and M. C. Barbosa, “Liquid polymorphism and density anomaly in a lattice gas model,” Phys. Rev. E 71, 031504 (2005).
  * Balladares, Henriques, and Barbosa (2007) A. L. Balladares, V. B. Henriques, and M. C. Barbosa, “Liquid polymorphism, density anomaly and H-bond disruption in associating lattice gases,” J. Phys.: Condens. Matter 19, 116105 (2007).
  * Bell (1972) G. M. Bell, “Statistical mechanics of water: lattice model with directed bonding,” J. Phys. C 5, 889 (1972).
  * Meijer, Kikuchi, and Royen (1982) P. H. E. Meijer, R. Kikuchi, and E. V. Royen, “Phase diagram of water based on a lattice model,” Physica A 115, 124 (1982).
  * Lavis and Southern (1984) D. A. Lavis and B. W. Southern, “Renormalization group study of a three-dimensional lattice model with directional bonding,” J. Stat. Phys. 35, 489 (1984).
  * Besseling and Lyklema (1994) N. A. M. Besseling and J. Lyklema, “Equilibrium properties of water and its liquid-vapor interface,” J. Phys. Chem. 98, 11610 (1994).
  * Roberts and Debenedetti (1996) C. J. Roberts and P. G. Debenedetti, “Polyamorphism and density anomalies in network‐forming fluids: Zeroth‐ and first‐order approximations,” J. Chem. Phys. 105, 658 (1996).
  * Roberts, Panagiotopoulos, and Debenedetti (1996) C. J. Roberts, A. Z. Panagiotopoulos, and P. G. Debenedetti, “Liquid-liquid immiscibility in pure fluids: Polyamorphism in simulations of a network-forming fluid,” Phys. Rev. Lett. 77, 4386 (1996).
  * Pretti and Buzano (2004) M. Pretti and C. Buzano, “Thermodynamic anomalies in a lattice model of water,” J. Chem. Phys. 121, 11856 (2004).
  * Pretti and Buzano (2005) M. Pretti and C. Buzano, “Thermodynamic anomalies in a lattice model of water: Solvation properties,” J. Chem. Phys. 123, 024506 (2005).
  * Sastry, Sciortino, and Stanley (1993) S. Sastry, F. Sciortino, and H. E. Stanley, “Limits of stability of the liquid phase in a lattice model with water‐like properties,” J. Chem. Phys. 98, 9863 (1993).
  * Borick, Debenedetti, and Sastry (1995) S. S. Borick, P. G. Debenedetti, and S. Sastry, “A lattice model of network-forming fluids with orientation-dependent bonding: Equilibrium, stability, and implications for the phase behavior of supercooled water,” J. Phys. Chem. 99, 3781 (1995).
  * Girardi  _et al._ (2007) M. Girardi, A. L. Balladares, V. B. Henriques, and M. C. Barbosa, “Liquid polymorphism and density anomaly in a three-dimensional associating lattice gas,” J. Chem. Phys. 126, 064503 (2007).
  * Oliveira, Stilck, and Barbosa (2010) T. J. Oliveira, J. F. Stilck, and M. A. A. Barbosa, “Solution of an associating lattice-gas model with density anomaly on a husimi lattice,” Phys. Rev. E 82, 051131 (2010).
  * Furlan, Fiore, and Barbosa (2015) A. P. Furlan, C. E. Fiore, and M. C. Barbosa, “Influence of disordered porous media on the anomalous properties of a simple water model,” Phys. Rev. E 92, 032404 (2015).
  * Buzzano, De Stefanis, and Pretti (2008) C. Buzzano, E. De Stefanis, and M. Pretti, “Cluster-variation approximation for a network-forming lattice-fluid model,” J. Chem. Phys. 129, 024506 (2008).
  * Pretti, Buzzano, and De Stefanis (2009) M. Pretti, C. Buzzano, and E. De Stefanis, “Revisiting waterlike network-forming lattice models,” J. Chem. Phys. 131, 224508 (2009).
  * Szortyka  _et al._ (2010) M. M. Szortyka, M. Girardi, V. B. Henriques, and M. C. Barbosa, “Dynamic transitions in a three dimensional associating lattice gas model,” J. Chem. Phys. 132, 134904 (2010).
  * Szortyka  _et al._ (2009) M. M. Szortyka, V. B. Henriques, M. Girardi, and M. C. Barbosa, “Dynamic transitions in a two dimensional associating lattice gas model,” J. Chem. Phys. 130, 184902 (2009).
  * Furlan, Almarza, and Barbosa (2016) A. P. Furlan, N. G. Almarza, and M. C. Barbosa, “Lattice model for water-solute mixtures,” J. Chem. Phys. 145, 144501 (2016).
  * Furlan  _et al._ (2019) A. P. Furlan, T. J. Oliveira, J. F. Stilck, and R. Dickman, “Order-disorder transition in a two-dimensional associating lattice gas,” Phys. Rev. E 100, 022109 (2019).
  * Wu (1982) F. Y. Wu, “The Potts model,” Rev. Mod. Phys. 54, 235–268 (1982).
  * Guo and Blöte (2002) W. Guo and H. W. J. Blöte, “Finite-size analysis of the hard-square lattice gas,” Phys. Rev. E 66, 046140 (2002).
  * Oliveira and Stilck (2015) T. J. Oliveira and J. F. Stilck, “Transfer-matrix study of a hard-square lattice gas with two kinds of particles and density anomaly,” Phys. Rev. E 92, 032101 (2015).
  * Fisher (1971) M. E. Fisher, in  _Proceedings of International School of Physics ”Enrico Fermi”_ , edited by M. S. Green (Academic, New York, 1971).
  * Blöte, Cardy, and Nightingale (1986) H. W. J. Blöte, J. L. Cardy, and M. P. Nightingale, Phys. Rev. Lett. 56, 742 (1986).
  * Affleck (1986) I. Affleck, Phys. Rev. Lett. 56, 746 (1986).
  * de Queiroz (2000) S. L. A. de Queiroz, “Finite-size scaling corrections in two-dimensional Ising and Potts ferromagnets,” J. Phys. A: Math. Gen. 33, 721 (2000).
  * Rikvold  _et al._ (1983) P. A. Rikvold, W. Kinzel, J. D. Gunton, and K. Kaski, “Finite-size-scaling study of a two-dimensional lattice-gas model with a tricritical point,” Phys. Rev. B 28, 2686 (1983).
  * Bartelt, Einstein, and Roelofs (1986) N. C. Bartelt, T. L. Einstein, and L. D. Roelofs, “Transfer-matrix approach to estimating coverage discontinuities and multicritical-point positions in two-dimensional lattice-gas phase diagrams,” Phys. Rev. B 34, 1616 (1986).
  * Xavier  _et al._ (1998) J. C. Xavier, F. C. Alcaraz, D. P. Lara, and J. A. Plascak, “Critical behavior of the spin-3/2 Blume-Capel model in two dimensions,” Phys. Rev. B 57, 11575 (1998).
  * Jung and Kim (2017) M. Jung and D.-H. Kim, “First-order transitions and thermodynamic properties in the 2D Blume-Capel model: the transfer-matrix method revisited,” Eur. Phys. J. B 90, 245 (2017).
  * Allen and Tildesley (1989) M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids_, Oxford Science Publ (Clarendon Press, 1989).
  * Debenedetti, Sciortino, and Zerze1 (2020) P. G. Debenedetti, F. Sciortino, and G. H. Zerze1, “Second critical point in two realistic models of water,” Science 369, 289 (2020).
  * Fiore  _et al._ (2009) C. E. Fiore, M. M. Szortyka, M. C. Barbosa, and V. B. Henriques, “Liquid polymorphism, order-disorder transitions and anomalous behavior: A monte carlo study of the Bell-Lavis model for water,” J. Chem. Phys. 131, 164506 (2009).
  * Simenas, Ibenskas, and Tornau (2014) M. Simenas, A. Ibenskas, and E. E. Tornau, “Phase transition properties of the Bell-Lavis model,” Phys. Rev. E 90, 042124 (2014).


# On the source contribution to the Galactic diffuse gamma rays above 398 TeV
detected by the Tibet AS$\gamma$ experiment

S. Kato  Corresponding author  Institute for Cosmic Ray Research, University
of Tokyo, Kashiwa 277-8582, Japan<EMAIL_ADDRESS>D. Chen
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China  J. Huang  Key Laboratory of Particle Astrophysics, Institute of
High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China  T.
Kawashima  Institute for Cosmic Ray Research, University of Tokyo, Kashiwa
277-8582, Japan<EMAIL_ADDRESS>K. Kawata  Institute for Cosmic
Ray Research, University of Tokyo, Kashiwa 277-8582, Japan;
<EMAIL_ADDRESS>A. Mizuno  Institute for Cosmic Ray Research,
University of Tokyo, Kashiwa 277-8582, Japan<EMAIL_ADDRESS>M.
Ohnishi  Institute for Cosmic Ray Research, University of Tokyo, Kashiwa
277-8582, Japan<EMAIL_ADDRESS>T. Sako  Institute for Cosmic Ray
Research, University of Tokyo, Kashiwa 277-8582, Japan;
<EMAIL_ADDRESS>T. K. Sako  Institute for Cosmic Ray Research,
University of Tokyo, Kashiwa 277-8582, Japan<EMAIL_ADDRESS>M.
Takita  Institute for Cosmic Ray Research, University of Tokyo, Kashiwa
277-8582, Japan<EMAIL_ADDRESS>Y. Yokoe  Institute for Cosmic
Ray Research, University of Tokyo, Kashiwa 277-8582, Japan;
<EMAIL_ADDRESS>
###### Abstract

Potential contribution from gamma-ray sources to the Galactic diffuse gamma
rays observed above 100 TeV (sub-PeV energy range) by the Tibet AS$\gamma$
experiment is an important key to interpreting recent multi-messenger
observations. This paper reveals a surprising fact: none of the 23 Tibet
AS$\gamma$ diffuse gamma-ray events above $398\,{\rm TeV}$ within the Galactic
latitudinal range of $|b|<10^{\circ}$ come from the 43 sub-PeV gamma-ray
sources reported in the 1LHAASO catalog, which proves that these sources are
not the origins of the Tibet AS$\gamma$ diffuse gamma-ray events. No
positional overlap between the Tibet AS$\gamma$ diffuse gamma-ray events and
the sub-PeV LHAASO sources currently supports the diffusive nature of the
Tibet AS$\gamma$ diffuse gamma-ray events, although their potential origin in
the gamma-ray sources yet unresolved in the sub-PeV energy range cannot be
ruled out.

Galactic cosmic rays (567) — Gamma-ray astronomy (628)

##  1 Introduction

The Tibet AS$\gamma$ experiment has performed the first detection of Galactic
diffuse gamma rays above $100\,{\rm TeV}$ (sub-PeV range) well concentrated on
the Galactic Plane (Amenomori et al., 2021). The hadronic origin of the Tibet
AS$\gamma$ diffuse gamma rays is supported by the observation of Galactic
neutrinos by IceCube (IceCube Collaboration et al., 2023); the energy spectrum
of the Tibet AS$\gamma$ diffuse gamma rays converted into the neutrino flux
assuming the hadronic interaction smoothly connects at $\simeq 60\,{\rm TeV}$
with the all-sky Galactic-neutrino flux scaled by the intensity of the diffuse
emission template proposed by Ackermann et al. (2012) in the field of view of
the Tibet AS$\gamma$ experiment. However, it is controversial whether they
have a totally diffusive nature or include some contributions from yet-
unresolved, dim gamma-ray sources. Fang & Murase (2021) suggest that some very
energetic sources such as hypernova remnants could largely contribute to the
Tibet AS$\gamma$ diffuse gamma-ray events. Vecchiotti et al. (2022) also
demonstrate that the Tibet AS$\gamma$ diffuse gamma-ray events can be
explained with the truly diffuse emission plus the emission from unresolved
pulsar-wind nebulae too dim to be detected, the latter of which could make a
significant contribution.

The possible pictures above can be tested experimentally with the LHAASO
observatory which is now steadily operating and surveying gamma rays from the
northern sky in the sub-PeV energy range with a sensitivity one order of
magnitude better than the other existing gamma-ray observatories (Cao et al.,
2021a, b). Considering the excellent sensitivity of LHAASO, one can naturally
expect their detection of many gamma-ray sources that the Tibet AS$\gamma$
experiment could not detect in the directions the Tibet AS$\gamma$ diffuse
gamma-ray events come from. Furthermore, the LHAASO collaboration recently
reported their first catalog (1LHAASO catalog) of gamma-ray sources detected
in $1\,{\rm TeV}<E<25\,{\rm TeV}$ and/or $E>25\,{\rm TeV}$, and 43 out of the
sources found in the latter energy range are also detected above $100\,{\rm
TeV}$ with more than $4\,{\sigma}$ levels (Cao et al., 2023). The detection of
the sub-PeV gamma-ray sources gives us a good opportunity to consider whether
these sources contribute to the Tibet AS$\gamma$ diffuse gamma-ray events.

##  2 Results and discussions

###  2.1 Tibet AS$\gamma$ Galactic diffuse gamma-ray emission events

This study focuses on the 23 gamma-ray-like diffuse events above $398\,{\rm
TeV}$ detected by the Tibet AS$\gamma$ experiment in the Galactic latitudinal
range of $|b|<10^{\circ}$, as reported in Amenomori et al. (2021); hereafter
let us simply call these events the Tibet AS$\gamma$ Galactic diffuse gamma-
ray emission (T-GDE) events. These events are found in the data taken by the
observatory during a live time of $719$ days from February 2014 to May 2017.
The Tibet AS$\gamma$ experiment consists of a surface air shower array
covering a geometrical area of $65\,700\,{\rm m}^{2}$ and an underground muon
detector array with a $3400\,{\rm m}^{2}$ area. The latter counts the number
of muons in air showers, leading to powerful discrimination of primary gamma
rays and cosmic rays. Amenomori et al. (2021) found the T-GDE events after the
event selection with the tight muon cut which realizes the cosmic-ray survival
ratio of $10^{-6}$ and the gamma-ray survival ratio of $30\%$ above $398\,{\rm
TeV}$. Note that the 23 events include $2.73$ background cosmic-ray events and
four gamma-ray-like events detected within $4^{\circ}$ from the center of the
Cygnus Cocoon; see Amenomori et al. (2021).

###  2.2 Positions of the T-GDE events and the sub-PeV LHAASO sources

Vecchiotti et al. (2022) propose that $\simeq 50\%$ of the diffuse gamma-ray
flux above $398\,{\rm TeV}$ measured by the Tibet AS$\gamma$ experiment could
be the contributions from gamma-ray sources unresolved by the Tibet AS$\gamma$
experiment. This means that $\sim 50\%$ of the T-GDE events come from
unresolved sources. If the 1LHAASO catalog sources detected above 100 TeV
(hereafter simply called sub-PeV LHAASO sources) are the origins of the T-GDE
events, about 10 sources should overlap the T-GDE events under the picture
proposed by Vecchiotti et al. (2022). Is that the case?

Figure 1 shows the incoming directions of the T-GDE events (red points) and
the positions and extensions of the sub-PeV LHAASO sources (blue circles; Cao
et al., 2023). Assuming a two-dimensional Gaussian morphology, the extensions
of $95\%$ containment radii are calculated from the published extensions of
the $39\%$ containment radii (Cao et al., 2023). For point-like sources, the
$95\%$ upper limits on the extensions are shown. Surprisingly, no T-GDE events
have their arrival directions within the extensions of the sub-PeV LHAASO
sources. This means that none of the sub-PeV LHAASO sources is the origin of
the T-GDE events.

###  2.3 Statistical consistency of the observations by the Tibet AS$\gamma$
experiment and LHAASO

One may think of the accidental positional overlap of the T-GDE events with
the sub-PeV LHAASO sources. The sum of the extensions (in terms of $95\%$
containment radii) of the sub-PeV LHAASO sources located in the area of
interest shown in Figure 1 covers $\simeq 3.7\%$ of the area of interest, and
the expected number of the T-GDE events that accidentally overlap the sub-PeV
LHAASO sources is estimated to be $\simeq 0.037\times 23=0.86$. Therefore, it
is natural to see no accidental overlap between the T-GDE events and the sub-
PeV LHAASO sources, as observed.

In the Tibet AS$\gamma$ diffuse-gamma-ray analysis performed by Amenomori et
al. (2021), they masked the gamma-ray sources listed in the TeVCat catalog as
of 2021. Twenty out of the 43 sub-PeV LHAASO sources are associated with the
TeVCat catalog sources masked in the analysis of Amenomori et al. (2021),
while the remaining 22 sources are not; hereafter such sources are called
newly-reported sub-PeV LHAASO sources, listed in Table 1. The statistical
consistency of the non-detection of gamma-ray events by the Tibet AS$\gamma$
experiment from these sources should be studied. The number of gamma-ray
events $n$ above 398 TeV from each of the newly-reported sub-PeV LHAASO
sources can be estimated as

$n=F(>398\,{\rm
TeV})\times\Bigg{(}\frac{\displaystyle\int_{\theta<40^{\circ}}\,D\,{\rm
cos}\,\theta(t)\,{\rm d}t}{\displaystyle\int_{\theta<40^{\circ}}\,{\rm
d}t}\Bigg{)}\times S_{\rm MD}\times T_{\theta<40^{\circ}}.$ (1)

Here $F(>398\,{\rm TeV})$ is the integral gamma-ray flux above $398\,{\rm
TeV}$ of the source. $F(>398\,{\rm TeV})$ is estimated from differential flux
and spectral index determined from the LHAASO KM2A measurement above 25 TeV,
because Cao et al. (2023) do not present the integral gamma-ray fluxes above
$398\,{\rm TeV}$ for the sub-PeV LHAASO sources. $D$ $(=65\,700\,{\rm m}^{2})$
is the geometrical area of the Tibet AS$\gamma$ experiment and $S_{\rm MD}$
$(=0.3)$ is the survival ratio of gamma rays above 398 TeV due to the event
selection using the underground muon detector array (Amenomori et al., 2021).
The time integrals are performed for the period in a sidereal day in which a
source of interest resides in the zenith-angle ($\theta$) range within
$40^{\circ}$ employed in the Tibet AS$\gamma$ data analysis. The term enclosed
with the parenthesis thus calculates the effective detection area for gamma-
ray events by integrating the geometrical area of the AS array projected to
the time-dependent source direction. $T_{\theta<40^{\circ}}$ is the total
period in which the source resides in $\theta<40^{\circ}$ over the live time
(719 days) of the dataset used in the Tibet AS$\gamma$ data analysis. The
results of the $n$’s for the newly-reported sub-PeV LHAASO sources are shown
in Figure 2 as a function of declination $\delta$. Since the newly-reported
sub-PeV LHAASO sources with $\delta<-10^{\circ}$ or $\delta>70^{\circ}$ do not
culminate in $\theta<40^{\circ}$ at the site of the Tibet AS$\gamma$
experiment ($30\fdg 102{\rm N}$), they have $n=0$ and are not shown in the
figure. The sum of the $n$’s of the newly-reported sub-PeV LHAASO sources is
$1.2$, so the non-detection of gamma-ray events by the Tibet AS$\gamma$
experiment from the newly-reported sub-PeV LHAASO sources is statistically
consistent. Therefore, our result currently supports the diffusive nature of
the T-GDE events, although the potential origin of the T-GDE events in the
gamma-ray sources yet unresolved in the sub-PeV energy range cannot be ruled
out.

Figure 1: Incoming directions of the Galactic diffuse gamma-ray events above
$398\,{\rm TeV}$ observed by the Tibet AS$\gamma$ experiment (red points) and
the positions and extensions of the sub-PeV LHAASO sources: the gamma-ray
sources detected by LHAASO above $100\,{\rm TeV}$ (blue circles; Cao et al.
(2023)) in the Galactic coordinates. The extensions of the sub-PeV LHAASO
sources have $95\%$ containment radii assuming a two-dimensional Gaussian
morphology. For point-like sources, the $95\%$ upper limits on the extensions
are shown. The green dashed circle with a $6^{\circ}$ radius encloses the
Cygnus region centered at $(l,\,b)=(79\fdg 62,\,0\fdg 96)$, the position of
${\rm HAWCJ}2030+409$ (Abeysekara et al., 2021). Note that some of the sub-PeV
LHAASO sources have $|b|>10^{\circ}$ and thus are not shown in the plot.
Figure 2: Expected number of gamma-ray events above 398 TeV that would be
observed by the Tibet AS$\gamma$ experiment from the 22 newly-reported sub-PeV
LHAASO sources. For more details, see the main text. Table 1: List of the
newly-reported sub-PeV LHAASO sources

Source Name | ${\delta}(^{\circ})$ | $n$
---|---|---
1LHAASO J0007+5659u | $57.00$ | $2.5\times 10^{-2}$
1LHAASO J0007+7303u | $73.07$ | $0$
1LHAASO J0056+6346u | $63.77$ | $4.7\times 10^{-2}$
1LHAASO J0206+4302u | $43.05$ | $7.5\times 10^{-2}$
1LHAASO J0212+4254u | $42.91$ | $6.0\times 10^{-2}$
1LHAASO J0216+4237u | $42.63$ | $6.3\times 10^{-2}$
1LHAASO J0343+5254u | $52.91$ | $3.0\times 10^{-2}$
1LHAASO J1740+0948u | $9.81$ | $2.6\times 10^{-2}$
1LHAASO J1809-1918u | $-19.30$ | $0$
1LHAASO J1814-1719u | $-17.89$ | $0$
1LHAASO J1814-1636u | $-16.62$ | $0$
1LHAASO J1825-1256u | $-12.94$ | $0$
1LHAASO J1825-1337u | $-13.63$ | $0$
1LHAASO J1839-0548u | $-5.8$ | $6.4\times 10^{-2}$
1LHAASO J1959+2846u | $28.78$ | $0.12$
1LHAASO J1959+1129u | $11.49$ | $5.5\times 10^{-2}$
1LHAASO J2002+3244u | $32.64$ | $3.8\times 10^{-2}$
1LHAASO J2020+3649u | $36.82$ | $0.12$
1LHAASO J2031+4052u | $40.88$ | $0.10$
1LHAASO J2108+5153u | $51.90$ | $0.16$
1LHAASO J2200+5643u | $56.73$ | $5.5\times 10^{-2}$
1LHAASO J2229+5927u | $59.55$ | $0.11$

Note. — $\delta$, declination of the source; $n$, the expected number of
gamma-ray events that would be observed by the Tibet AS$\gamma$ experiment
from the source during the live time, calculated from Equation (1). The
sources outside the field of view of the Tibet AS$\gamma$ experiment have
$n=0$; see the main text.

###  2.4 Positional correlation between the T-GDE events and gamma-ray sources
other than the sub-PeV LHAASO sources

This study checks overlaps between the T-GDE events and existing gamma-ray
sources other than the sub-PeV LHAASO sources to study the potential source
origin of the T-GDE events. First, the positional correlation between the
T-GDE events and the total 47 1LHAASO catalog sources detected only below 100
TeV is studied. The circular region with a radius of $0.5^{\circ}$ centered at
each of the 47 1LHAASO catalog sources is considered and the T-GDE events
within these circular regions are seached. The radius of the circular regions
is the same as that employed for the masked regions in the Tibet AS$\gamma$
diffuse-gamma-ray analysis (Amenomori et al., 2021). As a result, it is found
that only one T-GDE event named TASG-D01-025 with the arrival direction of
$(\alpha_{\rm J2000},\,\delta_{\rm J2000})=(286.96^{\circ},\,7.96^{\circ})$
(see the Supplemental Material of Amenomori et al. (2021)) is within the
circular region of 1LHAASO J1907$+$0826 which is detected only in $1\,{\rm
TeV}<E<25\,{\rm TeV}$ by the LHAASO Water Cherenkov Detector Array. The fact
of the one overlap is statistically consistent with an expected accidental
overlap (0.19 events).

Furthermore, overlaps between the T-GDE events and the gamma-ray sources
listed in the latest TeVCat catalog (Wakely & Horan, 2008) are checked. Out of
the total 273 sources listed in the latest TeVCat catalog 111As of 7th
December 2023, there are 92 sources in the sky region displayed in Figure 1
after removing extragalactic sources. Note that the 92 TeVCat catalog sources
include the sources masked in the Tibet AS$\gamma$ diffuse-gamma-ray analysis
(Amenomori et al., 2021), while they do not include the 1LHAASO catalog
sources. The circular region with a radius of $0.5^{\circ}$ centered at each
of the 92 TeVCat catalog sources is considered and the T-GDE events within
these circular regions are searched. As a result, there is no overlap between
the T-GDE events and the TeVCat catalog sources, and it is statistically
consistent with an expected accidental overlap (0.38 events). The results
presented in this section further support the diffusive nature of the T-GDE
events.

##  3 Conclusion

Our study shows that the 43 1LHAASO catalog sources detected above $100\,{\rm
TeV}$ (sub-PeV LHAASO sources) are not the origins of the 23 T-GDE events
above $398\,{\rm TeV}$ from the fact that no T-GDE events have their arrival
directions within the extensions of the sources. The number of accidental
overlaps between the T-GDE events and the sub-PeV LHAASO sources is estimated
to be $0.86$, and the expected number of gamma-ray events from the 22 newly-
reported sub-PeV LHAASO sources is $1.2$, both ensuring the statistical
consistency of no overlap between the T-GDE events and the sub-PeV LHAASO
sources. No overlap between the T-GDE events and the sub-PeV LHAASO sources
currently supports the diffusive nature of the T-GDE events, although the
potential origin of the T-GDE events in the gamma-ray sources yet unresolved
in the sub-PeV energy range cannot be ruled out.

This work is supported in part by Grants-in-Aid for Scientific Research from
the Japan Society for the Promotion of Science in Japan, the joint research
program of the Institute for Cosmic Ray Research (ICRR), the University of
Tokyo, and the use of the computer system of ICRR. This work is also supported
by the National Natural Science Foundation of China under Grants No. 12227804,
and the Key Laboratory of Particle Astrophysics, Institute of High Energy
Physics, CAS.

## References

  * Abeysekara et al. (2021) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2021, Nature Astronomy, 5, 465, doi: 10.1038/s41550-021-01318-y
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012, The Astrophysical Journal, 750, 3, doi: 10.1088/0004-637X/750/1/3
  * Amenomori et al. (2021) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2021, Phys. Rev. Lett., 126, 141101, doi: 10.1103/PhysRevLett.126.141101
  * Cao et al. (2023) Cao, Z., Aharonian, F., An, Q., et al. 2023, arXiv:2305.17030v2.  https://arxiv.org/pdf/2305.17030.pdf
  * Cao et al. (2021a) Cao, Z., Aharonian, F. A., An, Q., et al. 2021a, Nature, 594, 33, doi: 10.1038/s41586-021-03498-z
  * Cao et al. (2021b) Cao, Z., Aharonian, F., An, Q., et al. 2021b, Science, 373, 425, doi: 10.1126/science.abg5137
  * Fang & Murase (2021) Fang, K., & Murase, K. 2021, The Astrophysical Journal, 919, 93, doi: 10.3847/1538-4357/ac11f0
  * IceCube Collaboration et al. (2023) IceCube Collaboration, Abbasi, R., Ackermann, M., et al. 2023, Science, 380, 1338, doi: 10.1126/science.adc9818
  * Vecchiotti et al. (2022) Vecchiotti, V., Zuccarini, F., Villante, F. L., & Pagliaroli, G. 2022, The Astrophysical Journal, 928, 19, doi: 10.3847/1538-4357/ac4df4
  * Wakely & Horan (2008) Wakely, S. P., & Horan, D. 2008, 34th International Cosmic Ray Conference, 3, 1341\.  http://TeVCat.uchicago.edu