Weak first-order phase transitions in the frustrated square lattice $J_{1}$-$J_{2}$ classical Ising model

The study of minimal models is well known to provide valuable fundamental insights. A seminal example is the understanding of thermal phase transitions gleaned from the classical Ising model on the square lattice with only nearest-neighbor interactions. Adding next-nearest-neighbor (i.e., diagonal) interactions to that model provides us with a minimal model of frustrated magnetism. Its Hamiltonian is given by

where $\langle.\rangle$ and $\langle\langle.\rangle\rangle$ respectively denote nearest and next-nearest-neighbors, and $\sigma_{i}=\pm 1$. Frustration arises when $J_{2}>0$: the nearest-neighbor and next-nearest-neighbor terms can not be simultaneously minimized on any minimal plaquette of the lattice (regardless of the sign of $J_{1}$). While the phases of this model are well understood, a complete understanding of the order of its phase boundaries has eluded decades of effort [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33] (see Refs. 29, 30 for recent summaries).

The phase diagram is illustrated in Fig. 1). The phase transition at $g=0.5$ occurs at zero temperature [30, 34]. The low temperature phase is either ferromagnetic ($J_{1}<0$) or antiferromagnetic ($J_{1}>0$) when $0\leq g<0.5$, and stripe-ordered when $0.5<g<\infty$, where $g=J_{2}/|J_{1}|$. Both low temperature phases transition directly to a paramagnetic phase. Thus, this model provides an opportunity to study a phase transition between striped and paramagnetic phases.

Regarding $g>0.5$, the phase boundary must go to Ising universality when $g\rightarrow\infty$, where the model becomes two decoupled copies of the nearest-neighbor square lattice Ising model. The current consensus, based on Refs. [23, 24, 25, 26], is that of second-order transitions at $g\gtrsim 0.67$. The order of the transitions at $0.5\leq g\lesssim 0.67$ remains uncertain; there is suggestive evidence that at least part of it is weak first-order [11, 12, 14, 17, 18, 19, 20, 24, 27, 30, 32], but also suggestive evidence that it is entirely second order [15, 22].

Regarding $g<0.5$, the phase boundary must go to Ising universality when $g\rightarrow 0^{+}$. A few studies find suggestive evidence for first-order transitions in a contiguous range of $g$ up to $0.5$ [30, 24, 27, 31, 33], but strong results are lacking.

The difficulty in determining the order of the phase transitions in this model is essentially due to the difficulty in numerically distinguishing between weak first-order and second order phase transitions. We illuminate the salient features of weak first-order phase transitions that give rise to this difficulty via a comparison between strong and weak first-order phase transitions: Both strong and weak first-order phase transitions necessarily entail a discontinuity in the first derivative of the Helmholtz free energy and a discontinuity of the order parameter (see Fig. 2). Both also entail metastable branches of finite extent (terminating at spinodal points) in the phase diagram. As a spinodal point is approached along a metastable branch, the metastable state approaches criticality as if there were a critical point (termed a pseudo critical point [36, 37] ¹¹1This usage of “pseudo critical point” is different from the usage in Ref. [58]. or pseudospinodal point [39]) located beyond (or, in the mean-field limit, precisely at) the spinodal point [39]. Around a weak first-order phase transition, these pseudo critical points lie very close to the phase transition, which makes the coexistence region very small and imparts a large (but still finite) correlation length to the phase transition. Around a strong first-order phase transition, the pseudo critical points are far from the phase transition, and the correlation lengths are therefore small. For a weak first-order phase transition, the smallness of the coexistence region means that the hysteresis is difficult to resolve, and the large correlation lengths can make weak first-order phase transitions numerically appear as second-order phase transitions [40].

For the region $0<g\lesssim 0.67$ in the present model, all existing results that show first-order signatures leave open the possibility that those first-order signatures are only numerical artifacts (due to, for example, finite cluster size), and all existing results that show second-order signatures leave open the possibility that those second-order signatures are in reality due to weak first-order phase transitions. For the region $g\gtrsim 0.67$, the same is true for all results that appear prior to Ref. [23]. The authors of Ref. [23] (seemingly) make a breakthrough by combining numerical data with theoretical arguments instead of relying upon numerical data alone. Specifically, they (seemingly) rule out the possibility of first-order transitions at $g\gtrsim 0.67$ in the following way: because the present model has a theoretically known second-order phase transition at $g=\infty$, they [23] assume that there is a line of second-order transitions that connects to it, and they combine numerical data with theoretical arguments to conclude that this line of second-order transitions in the region $g\gtrsim 0.67$ maps on to a phase boundary of the Ashkin-Teller model that exhibits continuously varying critical exponents. Subsequent numerics [24, 25, 26] confirm signatures of this “Ashkin-Teller criticality” at $g\gtrsim 0.67$, hence the current consensus.

However, this consensus overlooks the fact that the theoretical arguments in Ref. [23] do not rule out an alternative possibility that is also consistent with the numerical data: a line of weak first-order phase transitions at $0.67\lesssim g<\infty$ that asymptotically turns into second-order transitions as $g\rightarrow\infty$. This is possible because a line of weak first-order phase transitions necessarily has a line of pseudo critical points on each side, and these three lines can asymptotically merge as $g\rightarrow\infty$; see Fig. 3. The mapping onto the Ashkin-Teller model would then apply to one or both of the pseudo critical lines instead of the phase boundary itself. Since weak first-order phase transitions can numerically appear as second-order phase transitions [40], the results [24, 25, 26] indicating Ashkin-Teller criticality at $g\gtrsim 0.67$ do not contradict this scenario. When $g\rightarrow\infty$, the two pseudo critical lines in this scenario become a single critical point, whose central charge must be $c=1$ for consistency with the Ashkin-Teller criticality of one or both of the pseudo critical lines. The central charge is a measure of the number of local degrees of freedom, and thus the two central charges of the pseudo critical lines do not add together at $g=\infty$ because the number of degrees of freedom does not change when $g\rightarrow\infty$. This is consistent with the fact that at $g=\infty$ the model is two decoupled copies of the ordinary Ising model, each of which has central charge $c=1/2$, and each of which also has only half the number of degrees of freedom of the full model.

Below we present strong evidence that this alternative scenario is the actual case. Further, we present strong evidence that the phase transitions in $0.5<g\lesssim 0.67$ are also first-order. These results comprise strong evidence for first-order phase transitions in this model; the evidence is strong instead of only suggestive because we take care to show that it is highly unlikely that the first-order signatures that we find are mere numerical artifacts.

Additionally, from the discussion above and from Fig. 2 it seems plausible that there can be a first-order phase transition where the correlation length is very large in one phase but very small in the other; we find suggestive evidence that this is what occurs as $g\rightarrow 0.5^{+}$ in this model. Such a first-order phase transition can neither be categorized as a weak first-order phase transition nor as a strong first-order phase transition. We are not aware of any previous suggestions in the literature of such a scenario.

In Sec. II we explain the numerical method that we employ. In Sec. III we present numerical data that reveals the signatures of weak and anomalous first-order phase transitions. In Sec. IV we present further data and theoretical arguments to show that the first-order signatures are not numerical artifacts. In Sec. V we conclude with a discussion.

The class of computational methods designed around tensor networks [41, 42, 43, 44, 45, 46] is useful for studying both classical and quantum spin lattices; the following results from that field provide the basis for the numerical method of the present work: Ref. [47] demonstrates with a quantum 2D spin lattice that adiabatic evolution of 2D tensor networks is able to reveal ground state energy density level crossings and order parameter discontinuities directly in the thermodynamic limit. Ref. [48] demonstrates how to use a 1D tensor network, known as a matrix product state (MPS), to compute local observables directly in the thermodynamic limit of a classical 2D spin lattice. This MPS approach amounts to a power method wherein the row-to-row transfer operator of the partition function is repeatedly applied to the MPS until the MPS is sufficiently converged toward the dominant eigenvector of the transfer operator; direct access to the thermodynamic limit is provided by presuming translation invariance up to a chosen unit cell size and computing with only a single unit cell. The converged MPS may be interpreted as the groundstate wavefunction of a quantum 1D system that is dual to the classical 2D system [49, 50], and the correlations in the MPS may be interpreted as the bipartite entanglement in that wavefunction. The entanglement capacity of an MPS is limited by the size of its matrices, and this approach therefore suffers from finite-entanglement effects instead of finite-size effects.

We use the tensor network representation of the partition function for the present model that is given in the Appendix of Ref. [51]. We compute the Helmholtz free energy density and order parameter with an MPS method similar to that in Ref. [48].

The converged MPS gives ready access to the Schmidt decomposition of the one-dimensional wavefunction that it represents:

where $|\phi_{i}^{L}\rangle$ ($|\phi_{i}^{R}\rangle$) are states of the left (right) half of a bipartition of the system, $\lambda_{i}$ are real scalars in descending order, $\sum_{i}\lambda_{i}^{2}=1$, $\langle\phi_{i}^{L}|\phi_{j}^{L}\rangle=\langle\phi_{i}^{R}|\phi_{j}^{R}\rangle=\delta_{ij}$, and $\chi\in\mathbb{Z}^{+}$ (known as the bond dimension) limits the entanglement capacity of the MPS and also determines the computational cost of its storage and manipulation. A single application of the transfer operator to the MPS results in an increase of the bond dimension to $4\chi$, after which the Schmidt decomposition must be truncated (and renormalized) to keep the first $\chi$ terms only; the sum $\sum_{i=\chi+1}^{4\chi}\lambda_{i}^{2}$ is known as the truncation error.

The converged MPS also gives ready access to the correlation length [44]. While a more accurate method for computing the correlation length from the MPS exists [52], the simple method in Ref. [44] is sufficient for our goal here of arriving at only a qualitative understanding.

For each value of $g$, we perform independent adiabatic evolutions with fixed values of $\chi$ using a two-site unit cell (represented by an MPS with two tensors, one for each physical site). For each adiabatic evolution in temperature at a given value of $g$, we first converge a randomly initialized MPS with the transfer operator at the starting temperature, compute the physical quantities with the converged MPS, then use that converged MPS as the initial MPS for convergence at the subsequent temperature; the rest of the adiabatic evolution is completed by using the converged MPS at each temperature as the initial MPS for the subsequent temperature. The entanglement entropy, given by

depends on both the local and nonlocal correlations in the MPS, and we therefore consider the MPS converged at a given temperature when $S$ is converged to within a chosen tolerance ($\epsilon$). Adiabatic evolutions in $g$ at fixed temperature are done in an analogous way. This method becomes numerically exact in the limit of vanishing increments of $T$ or $g$, $\chi\rightarrow\infty$, and $\epsilon\rightarrow 0$. We note that a similar adiabatic evolution method for probing phase boundaries with MPSs is used in Refs. [53, 54], the difference being that their approach uses a variational method instead of a power method to converge the MPSs.

The order parameter $(m_{x},m_{y})$ is defined by

where $N_{\textrm{x}}$ ($N_{\textrm{y}}$) is the number of lattice sites in the $\hat{x}$ ($\hat{y}$) direction, and ($x,y$) are the coordinates of a spin. The stripe-ordered phase can have either $~{}m_{x}>0$ and $~{}m_{y}=0$ or $~{}m_{x}=0$ and $~{}m_{y}>0$.

We check for Helmholtz free energy density crossings at $g=0.49$ but find none. This by itself does not rule out the possibility of extremely weak first-order transitions at $g=0.49$; we return to this in Sec. V.

In Fig. 4) we illustrate the results of high-to-low $T$ and low-to-high $T$ adiabatic evolutions at $g=0.51$ with $\chi=64$ and $\epsilon=1$e-6. The defining features of a first-order phase transition are apparent: a Helmholtz free energy density crossing occurs and the order parameter for the equilibrium state is discontinuous at the crossing.

The apparent energy density crossing in Fig. 4) can not correspond to an energy density gap closing of a second-order transition for the following reasons: (1) As the gap-closing of a second-order transition is adiabatically approached from either side, the MPS must approach the same state and therefore the entanglement entropy ($S$) must approach the same value, whereas in Fig. 4) the two values of $S$ at the crossing differ greatly. (2) If the crossing is really a gap closing of a second-order transition, then the (quasi-)adiabatic evolution jumps from the ground state manifold to the excited state manifold as the crossing is traversed, which must result in a discontinuity in $S$ along either adiabatic path and a discontinuity in $(m_{x},m_{y})$ over the low-to-high $T$ path, but we find that both $S$ and $(m_{x},m_{y})$ are continuous over the energy density crossing along either adiabatic path. (3) A finite bond dimension imposes a finite correlation length [44] and therefore a finite gap, but there is no sign of a gap at the crossing point.

We perform adiabatic evolutions with different values of $\chi$ and $\epsilon$ at a series of $g$. With increasing $g$ at fixed $\chi$, the order parameter discontinuity (see Fig. 5) and the angle of incidence between the two free energy curves (see Fig. 6) both decrease. This is consistent with the requirement that the order parameter discontinuity and Helmholtz free energy crossing both disappear when $g\rightarrow\infty$. Though the discontinuity at $g=0.7$ is not very small, this is still compatible with the possibility that previous numerics mistakenly find second-order transitions there; the weak first-order phase transition of the classical square lattice 5-state Potts model has a large order parameter discontinuity (about 0.49) [35], yet appears as a second-order transition in certain finite-size numerics [40] ²²2A weak first-order transition may have a large order parameter discontinuity if the order parameter curve drops very steeply such that the energy crossing occurs at a large value of the order parameter yet is still very close to the pseudo critical point that lies at the zero point of the order parameter curve..

We also compute the correlation length in each phase near the phase boundary, as shown in Figs. (7), (8), and (9). We obtain the data in Fig. 9) with $\epsilon=1$e-4. Converging $S$ to within this value of $\epsilon$ is sufficient for the correlation length data near the phase boundary because near criticality, $S\approx(c/6)\textrm{log}(\xi)$ [56], where $c$ is the central charge and $\xi$ is the correlation length; $\epsilon=1$e-4 yields an error in $\xi$ of at most $\sim 10$.

The large correlation lengths and their hysteresis that we find near the phase boundary are consistent with weak first-order transitions in the region $g\geq 0.501$, and are also consistent with the numerical detection of criticality in the region $0.67\lesssim g\lesssim 1$ in Refs. [25, 26] because, as mentioned above, weak first-order phase transitions can numerically appear as second-order phase transitions. The hysteresis (assuming it is not a numerical artifact), however, is not consistent with actual second-order transitions. Interestingly, the data in Fig. 7) suggests that as $g\rightarrow 0.5^{+}$, a first-order transition arises that is neither strong nor weak but part of both: the correlation length is very large in one phase but very small in the other. Further investigation of this is left for future work.

Thus, we find numerical evidence for first-order phase transitions in the region $0.5<g\leq 10$. This is at odds with the current consensus of second-order transitions in the region $g\gtrsim 0.67$, but is consistent with the alternative possibility explained in Sec. I of first-order transitions in the region $0.5<g<\infty$. However, our numerical evidence necessarily contains contributions from numerical imperfections, such as finite bond dimension. In the next section we show that such contributions are not responsible for the first-order signatures that we find, and thereby that our evidence for first-order phase transitions in the region $0.5<g\leq 10$ is strong. We also explain why this leads to the conclusion that the first-order phase boundary only asymptotically becomes second-order as $g\rightarrow\infty$.

First-order phase transitions have long been suspected to occur in this model, yet strong evidence has remained absent. We addressed this in the following ways: We pointed out that the numerical data and theoretical arguments that form the basis for the current consensus of Ashkin-Teller criticality in the region $g\gtrsim 0.67$ leave open the possibility of an alternative scenario: weak first-order phase transitions at $g\gtrsim 0.67$ that asymptotically become second-order as $g\rightarrow\infty$. We presented data and theoretical arguments to provide strong evidence that this alternative scenario is the true one, and further that the phase boundary is also first-order in the region $0.5<g\lesssim 0.67$. In contrast to previous studies that presented only suggestive evidence for first-order transitions in this model, we established the first strong evidence by taking care to account for contributions from numerical imperfections. This was facilitated in part by the theory of finite entanglement scaling, which allowed us to make conclusions about the limit of infinite bond dimension ($\chi$) from data obtained with only modest $\chi$. Additionally, we provided suggestive evidence that as $g\rightarrow 0.5^{+}$, the first-order transitions become of a previously unobserved type that can be categorized as neither strong first-order nor weak first-order: the correlation length is very large in one coexisting phase and very small in the other.

The numerical correlation lengths at the phase boundary at $g>0.5$ are determined in part by the proximity of the pseudo critical lines to the phase boundary and in part by the value of $\chi$. The locations of the pseudo critical lines shift with $\chi$, as required by Eq. (9). Therefore, in order to confirm the anomalous first-order transitions at $g\rightarrow 0.5^{+}$ with the present approach of adiabatic evolutions of MPSs, values of $g$ in the region $0.5<g<0.501$ must be investigated showing convergence with increasing $\chi$ of both the correlation length data and the location of the pseudo critical lines. If the anomalous first-order transitions are confirmed in this model, it would be important to consider in which other models, both classical and quantum, they may also appear.

Regarding the pseudo critical lines themselves, while it is clear that they must merge with the phase boundary when $g\rightarrow\infty$, the following remains unclear: where the pseudo critical lines go when $g\rightarrow 0.5^{+}$, which pseudo critical line exhibits the Ashkin-Teller criticality at $g\gtrsim 0.67$, and what the criticality of either pseudo critical line is in the region $0.5<g\lesssim 0.67$. Recent results [32] indicate that at least one of the pseudo critical lines continues to have continuously varying critical exponents in the region $0.5<g\lesssim 0.67$. It may be possible to use the results in Refs. [58, 59, 60] to compute the critical exponents and central charge along the pseudo critical lines with MPSs.

We left open the possibility of an extremely weak first-order transition at $g=0.49$. The fate of the pseudo critical lines at $g\rightarrow 0.5^{+}$ has a bearing on this: If the pseudo critical line for the disordered phase (i.e., the pseudo critical line at lower $T$) merges with the phase boundary as $g\rightarrow 0.5^{+}$, it would give an infinite correlation length to the phase transition at $g=0.5$ from the high-$T$ side, and would thereby support the idea of a second-order transition at $g=0.49$.

Finally, Ref. [66] showed good results for numerically detecting weak first-order phase transitions in other models with a different approach than the one used here. Future work may be able to use that approach as an independent check of the results that we have presented here.

As mentioned in Sec. II, the converged MPS at a given point in the phase diagram of the classical 2D system is an approximate representation of the ground state of the quantum 1D system that is dual to the classical 2D system at that point in the phase diagram. For an infinite quantum 1D system near a conformally invariant critical point, Ref. [56] shows that the half chain entanglement entropy is given by

where $c$ is the central charge and $\xi_{q}$ is the correlation length of the quantum 1D system. Because of the duality, $c$ is the same between the quantum 1D system and its dual classical 2D system. Near criticality we also have $\xi_{q}=\delta^{-\nu}$ and $\xi_{cl}=\tau^{-\nu}$, where $\xi_{cl}$ is the correlation length of the classical 2D system, $\delta$ and $\tau$ are small parameters, and the critical exponent $\nu$ is the same in both systems due to the duality. More precisely, $\tau=|T-T_{c}|/T_{c}$, where $T_{c}$ is the temperature of the critical point, and the duality tells us that $\delta$ is a differentiable function of $\tau$. Using a Taylor expansion around $\tau=0$ we can write $\delta(\tau)=\delta(0)+a\tau+\mathcal{O}(\tau^{2})$, where $a$ is a constant. Noting that $\xi_{q}=\xi_{cl}=\infty$ at the critical point (i.e., where $\tau=0$), we find $\delta(0)=0$. Close to the critical point we can therefore write

Since the half-chain entanglement entropy $S$ can be computed efficiently from the converged MPS, we may fit Eq. (11) to the numerical data obtained by adiabatically approaching the phase boundary from both high and low temperature at fixed values of $g$. Fig. 16 shows an example. For a continuous phase transition, both $\frac{\nu c}{6}$ and $T_{c}$ must be the same between the high- and low-temperature sides of the phase transition (though an exception can occur in anisotropic models, where it is possible that $\nu$ can be different between the two sides [67]). For a weak first-order phase transition, $T_{c}$ lies at the pseudo critical points and must therefore be different between the high- and low-temperature sides of the transition (see Fig. 2). It is also possible for $\frac{\nu c}{6}$ to be different between the two sides of a weak first-order transition.

We use the known second-order phase transition at $g=0$ ($J_{2}=0$, $J_{1}=1$) as a benchmark case for the fitting of $T_{c}$ and $\frac{\nu c}{6}$. We converge $S$ to within $1$e-$6$ at each temperature step. We do not find good agreement between $\frac{\nu c}{6}$ (not shown) as extracted from fits of Eq. (11) on the high- and low-temperature sides of the phase transition at $g=0$, but we do find agreement of $T_{c}$ to high precision (see Table 1 and Fig. 17). We therefore only use fitted values of $T_{c}$ to draw conclusions at other values of $g$ (see Sec. IV).

At $g=\infty$ the present system is two copies of the nearest neighbor classical Ising model on the square lattice, each of which has the Hamiltonian

where $J_{cl}$ is a real scalar and the other variables are as defined in the Introduction, and whose critical temperature is given by [68]

where $k_{B}$ is Boltzmann’s constant.

If we write the nearest neighbor quantum 1D Ising model with a transverse magnetic field as

then it has a quantum phase transition at $\lambda=\lambda^{*}=1$ [69].

From Ref. [70] we know that there is a mapping between $H_{cl}$ and $H_{q}$ such that

Thus, the value of $\lambda$ in Eq. (14) fixes the value of $J_{cl}$ in Eq. (12). Further, the mapping is such that both Hamiltonians become critical simultaneously, which fixes the critical value of $J_{cl}$ at a value that we denote as $J_{cl}^{*}$:

As mentioned in Sec. IV, in a finite-$\chi$ MPS representation with sufficiently large $\chi$, the quantum 1D critical point is shifted to a value $\lambda^{*}_{\chi}$ so that:

Thus, combining Eqs. (13) and (17), we find that a monotonic shift in $\lambda_{\chi}^{*}$ as a function of $\chi$ results in a monotonic shift in $T_{c,\chi}$ as a function of $\chi$.

This establishes that the phase boundary at $g=\infty$ shifts smoothly and monotonically as a function of $\chi$ when $\chi$ is sufficiently large (i.e., within the finite-$\chi$ scaling limit). The same must be true for the entire phase boundary at $g>0.5$, otherwise it would be unphysically crooked when $\chi\rightarrow\infty$.

The correlation length data in Fig. 9 shows that $\chi\geq 64$ falls within the finite-$\chi$ scaling limit at $g\geq 10$. This is consistent with the fact that the model at $g=10$ is only a perturbation of the model at $g=\infty$, and at $g=\infty$ the monotonic shift of the critical point with $\chi$ that is required by Eq. (9) is valid all the way down to $\chi=2$ [58]. It is also consistent with the numerical data in Figs. 18 and 19, which shows that the phase boundary around $g=10$ monotonically shifts, in agreement with Eq. (9), to larger values of $g$ and smaller values of $T$ as $\chi$ is increased from $16$. Also, it is consistent with Eq. (7) in the following sense: At $g=\infty$ we have exactly $\xi_{\chi}=\chi^{\kappa}$ [58], so we expect the prefactor $a$ in $\xi_{\chi}=a\chi^{\kappa}$ at $g=10$ to be close to unity since the model at $g=10$ is only a perturbation of the model at $g=\infty$. If $g=10$ has a second-order phase transition, then previous results [23, 24, 25, 26] establish that it has central charge $c=1$, which yields $\kappa\approx 1.344$. Then at $g=10$ we easily satisfy $\xi_{\chi}\gg 1$ when $\chi\geq 64$.

We therefore conclude that the entire phase boundary at $g\geq 10$ monotonically shifts to larger values of $g$ as $\chi$ is increased above $64$.

Eq. (9) estalishes that conformally invariant quantum 1D critical points that appear within the finite-$\chi$ scaling limit (i.e., when $\chi$ is large enough so that $\xi_{\chi}\gg 1$) can not split into two pseudo critical points as $\chi$ is varied within the finite-$\chi$ scaling limit. The results in Ref. [49] establish that the critical points of a given classical 2D spin lattice are equivalent to the critical points of some quantum 1D spin lattice. Therefore, conformally invariant critical points that appear in the present classical 2D model within the finite-$\chi$ scaling limit also can not split into two pseudo critical points as $\chi$ is varied within the finite-$\chi$ scaling limit.

If the present classical 2D model has a line of critical points at $g\geq 10$ when $\chi=\infty$, then Refs. [23, 24, 25, 26] establish that it must have conformal invariance with $c=1$. Also, the correlation length data in Fig. 9 reveals $\xi_{\chi}\gg 1$ at $g=10$ when $\chi\geq 64$, so $\chi\geq 64$ lies within the scaling limit at $g=10$. Therefore, Eq. (9) must apply to any line of critical points in the region $g\geq 10$ when $\chi\geq 64$. This is consistent with the fact that $g=10$ is only a perturbation of the model at $g=\infty$, the critical point of the model at $g=\infty$ maps, via the classical-quantum correspondence [49], on to the critical point of the quantum 1D Ising model, and Eq. (9) is valid for the latter all the way down to $\chi=2$ [58]. But Eq. (9) forbids a critical point from splitting into multiple pseudo critical points. It is therefore forbidden for a line of critical points that appear at $g\geq 10$ when $\chi=\infty$ to split into two pseudo critical lines as $\chi$ is reduced down to $64$.