Ashkin-Teller phase transition and multicritical behavior in a classical monomer-dimer model

The usual dimer worm algorithm [27, 26, 17] provides a rejection-free nonlocal update scheme for interacting dimer models at full-packing. Here, we build on ideas developed in Ref. 32 to generalize this dimer worm algorithm and obtain an efficient grand-canonical algorithm for the monomer-dimer model at nonzero monomer fugacity. In the first step of our grand-canonical scheme, one chooses at random a site $j_{\rm init}$ of the lattice. There are two possibilities at this first step: either the initially chosen site $j_{\rm init}$ has a monomer on it, or it is covered by a dimer. Let us consider each in turn.

If $j_{\rm init}$ has a monomer on it, we have five options at our disposal. The first four options consist of placing a dimer connecting the initial site to one of its four neighbors. The fifth option is to exit without doing anything. Each of these possibilities is assigned a probability from a probability table. We will discuss the construction of this probability table in some detail below. For now we simply introduce some language that will subsequently be useful in describing the construction of this probability table: the initial site is our first “pivot” site $\pi_{0}$, which we have “entered” from the “entry” $\sigma_{0}=0$, i.e. from “outside the lattice” (Fig. 3). Aborting our attempted worm move at this step itself without doing anything corresponds to “exiting the pivot” $\pi_{0}$ via “exit” $\sigma^{\prime}_{0}=0$. On the other hand, if we opt to place a dimer connecting the pivot $\pi_{0}$ to its $k$-th neighbor, this option corresponds to exiting the pivot via exit $\sigma^{\prime}_{0}=k$ (so $k$ can take on values from 1 to 4). If the chosen exit is $\sigma^{\prime}_{0}\neq 0$ , we now move to the site corresponding to the chosen exit and continue the construction.

Before we describe what is done next, we need to specify the procedure to be used if the initially chosen site $j_{\rm init}$ has a dimer covering it. In this case, one walks to the other end of this dimer; the site covered by this other end becomes our first pivot $\pi_{0}$, which we have “entered” from the entry $\sigma_{0}$ corresponding to $j_{\rm init}$. Now, the choices available are again five in number: One can delete this dimer that connects the pivot site $\pi_{0}$ to the entrance site $j_{0}$. This introduces two monomers in the system and concludes the worm move. As before, this corresponds to “exiting” the first pivot $\pi_{0}$ via exit $\sigma^{\prime}_{0}=0$. Or, one can pivot the dimer covering $\pi_{0}$ so that it now connects $\pi_{0}$ to its $k$-th neighbor. If one of these latter four options is chosen, we say the pivot $\pi_{0}$ is exited via exit number $\sigma^{\prime}_{0}=k$, and we move to the site corresponding to the chosen exit to proceed further as described below.

At this stage of our worm construction, we are at the site corresponding to exit $\sigma^{\prime}_{0}$ of the previous pivot, having arrived there because we chose to place a dimer connecting the previous pivot point $\pi_{0}$ to this exit site. If this site does not already have another dimer covering it, we have reached an allowed configuration and the worm move ends. On the other hand, if this exit site does have another dimer already covering it, this site becomes the current “overlap site” $o_{0}$. We now walk along this pre-existing dimer from $o_{0}$ to its other end. The site at this other end becomes our next pivot site $\pi_{1}$, which has been “entered” via entry number $\sigma_{1}$ that corresponds to the overlap site $o_{0}$.

At this step, there are again five choices for $\sigma^{\prime}_{1}$, the exit to be used to exit the current pivot site $\pi_{1}$. As before, exit $\sigma^{\prime}_{1}=0$ corresponds to deleting the dimer covering the current pivot site $\pi_{1}$. If this is chosen, the worm move ends. On the other hand, exits numbered $\sigma^{\prime}_{1}=1$ through $\sigma^{\prime}_{1}=4$ correspond to pivoting the dimer covering $\pi_{1}$ so that it now connects $\pi_{1}$ to the site corresponding to $\sigma^{\prime}_{1}$. If this exit site does not already have a dimer covering it, we have reached an allowed configuration and the move ends. Otherwise, this exit site becomes the next overlap site $o_{1}$, and the procedure is repeated.

It is easy to see that this worm construction yields a valid rejection-free algorithm if we choose the probability table $P_{\sigma\rightarrow\sigma^{\prime}}$ for transition probabilities in a way that it satisfies local detailed balance at each step. This amounts to requiring that the probability obeys the constraint equations:

Here, $P_{\sigma\rightarrow\sigma^{\prime}}$ is the conditional probability for exiting a pivot via exit $\sigma^{\prime}$ given that we have entered it via entrance $\sigma$, $P_{\sigma^{\prime}\rightarrow\sigma}$ is the conditional probability for the reverse process, and the weights $\omega_{\sigma}$ and $\omega_{\sigma^{\prime}}$ represent the Boltzmann weights of the configurations corresponding to the choices $\sigma$ and $\sigma^{\prime}$, respectively. These Boltzmann weights are to be calculated ignoring the violation of the hard-core constraint on dimers in the configurations that arise during the worm construction.

The simplest choice of solution is the heat-bath solution (sometimes called the Gibbs sampler) given as $P_{\sigma\rightarrow\sigma^{\prime}}=\omega_{\sigma^{\prime}}/\sum_{\sigma^{\prime}}\omega_{\sigma^{\prime}}$. In practice, we use the iterative Metropolized Gibbs sampler to reduce the bounce process [33], i.e., reduce the magnitude of the diagonal elements of the probability table. Note also that the computation of the weights $\omega_{\sigma}$ and $\omega_{\sigma^{\prime}}$ is simplified by the fact that they only differ due to factors arising from the contribution of the immediate neighborhood of the pivot. Since the equation set is homogenous, we can cancel all common factors to define reduced weights that only depend on the local environment of the pivot, and use these in Eq. (2).

These reduced weights can be written as

where $n$ ($m$) denotes the number of nearest neighbor dimers parallel in the transverse (longitudinal) direction to the dimer that covers the pivot when the configuration corresponds to entrances/exits $\sigma\neq 0$. These numbers $n,m\in\{0,1,2\}$ can be calculated by checking the direction of dimers on eight sites around $\bm{r}_{\text{h}}$, that is, $\bm{r}_{\text{h}}\pm\bm{e}_{x}$, $\bm{r}_{\text{h}}\pm\bm{e}_{y}$, $\bm{r}_{\text{h}}\pm 2\bm{e}_{x}$, $\bm{r}_{\text{h}}\pm 2\bm{e}_{y}$. The number of valid configurations in these eight sites is $65089$, much smaller than $5^{8}=390625$.

The factor of $z^{2}$ in the first term of Eq. (3) reflects the fact that configuration $\sigma=0$ has one fewer dimer, i.e., two additional holes (monomers) in comparison with the configurations with $\sigma\neq 0$.

Finally, we note for completeness that this worm construction and its detailed balance property generalizes straightforwardly to lattices with arbitrary coordination number. The “out-of-plane” entrance/exit in the general case is numbered $0$, and the other entrances/exits are numbered from $1$ to $n_{c}$, where $n_{c}$ is the coordination number of the pivot site in question ($n_{c}$ can be different for different sites, and there is thus no restriction of regularity for this algorithm to remain valid)

Using this algorithm, we perform the MC simulations of $N=L\times L$ systems with periodic boundary conditions for system size $L$ up to $L=512$. The number of the worm updates $n_{w}$ used per Monte Carlo step is chosen for each set of control parameters to be such that $n_{w}\langle l_{w}\rangle=N$, $\langle l_{w}\rangle$ is the mean number of sites visited during the construction of a single worm. We perform $10^{3}\times N$ worm updates to estimate $\langle l_{w}\rangle$ and to thermalize the system before measuring physical quantities. With this convention defining a MC step, we ensure that we obtain at least $2\times 10^{6}$ MC configurations of the system from which we can calculate equilibrium properties.

The tensor network representation of our model is based on the singular value decomposition of the local Boltzmann weight on a bond. The partition function is rewritten as the contraction of the tensor network,

The tensor $A$ is located on sites of the square lattice and has four indices representing links to the nearest neighbor sites. An element of $A$ has the form

where $s$ denotes the local configuration at a site as shown in Fig. 3. The matrices $X$’s are determined by the singular value decomposition of the local Boltzmann weight on a bond. The Boltzmann weight on a horizontal bond is represented as a $5\times 5$ matrix,

where the row (column) index of $W_{h}$ corresponds to a state at $\bm{r}$ ($\bm{r}+\bm{e}_{x}$), respectively. Elements with a value of zero indicate a configuration prohibited by the hard-core constraint. The singular value decomposition, $W_{h}=U_{h}S_{h}V_{h}^{T}$, defines $X_{r}\equiv U_{h}S_{h}^{1/2}$ and $X_{l}\equiv V_{h}S_{h}^{1/2}$. We note that $U_{h}$ and $V_{h}$ can be chosen to be real orthogonal $5\times 5$ matrices since $W_{h}$ is a real square matrix. Similarly we obtain the Boltzmann weight on a vertical bond between $\bm{r}$ (row) and $\bm{r}+\bm{e}_{y}$ (column) as

and we define $X_{t}\equiv U_{v}S_{v}^{1/2}$ and $X_{b}\equiv V_{v}S_{v}^{1/2}$. The chemical potential of a monomer acts as the external field and gives the corresponding on-site factor as

The first term has the factor of $z$ in contrast to Eq. (3). The former corresponds to the Boltzmann weight for valid configurations of our monomer-dimer model, while tha latter is for extended configurations containing one doubly occupied site.

The row-to-row transfer matrix with infinite width is represented as a uniform matrix product operator with a local tensor $A$. Using the variational uniform matrix product state algorithm (VUMPS) [28, 29, 30, 31], we calculate the eigenvector corresponding to the largest eigenvalue of the transfer matrix. The uniform matrix product state (uMPS) obtained in this way approximates this eigenvector with accuracy that is controlled by the bond dimension $\chi$ of the uMPS. We increase $\chi$ up to $128$ to ensure sufficient accuracy. In practice, we assume that the uMPS has a $2\times 2$ unit-cell structure [31] as is appropriate for a description of the columnar ordered state. After calculating the horizontal uMPS in this way, we also calculate the vertical uMPS, which approximates the corresponding eigenvector of the column-to-column transfer matrix. A good initial guess for the vertical uMPS can be given by the fixed point tensor of the horizontal uMPS [28]. We have confirmed that the physical quantities calculated from the horizontal and vertical uMPS agree with each other to machine precision even in the ordered phases.

We detect columnar order using a complex order parameter constructed from the following local order parameter field defined at each site $\bm{r}=(r_{x},r_{y})$ as

The corresponding order parameter,

takes $\pm 1$ or $\pm i$ when the state has the complete columnar order.

Nematic order, which breaks the symmetry of $\pi/2$ rotations, can be detected by comparing the number of horizontal and vertical dimers. With this motivation, we define the local nematic order parameter field as

The corresponding order parameter,

takes on values $\pm 1$ both in the nematic and columnar states.

The corresponding Binder ratios [34] are defined in the usual way:

As is well-known, the Binder ratios converge in the ordered phase to $1$ as the denominator and numerator take on the same limiting value in the thermodynamic limit. On the other hand, in a phase without symmetry breaking, the limiting value depends on the nature of fluctuations of the order parameter. $U_{\text{col}}$ in the nematic and $U_{\text{nem}}$ in the disordered fluid phase are both expected to converge to $3$ in the thermodynamic limit because the fluctuations of the corresponding order parameters obey a one-dimensional Gaussian distribution in these regimes. On the other hand, $U_{\text{col}}\rightarrow 2$ in the disordered fluid phase because the fluctuations of $m_{\text{col}}$ obey a two-dimensional Gaussian distribution. We note that the conventional and Ising definition of the Binder parameter for the frustrated Ising model discussed in Ref. [35] correspond $U_{\text{col}}$ and $U_{\text{nem}}$, respectively. Since the Binder ratio is a dimensionless quantity in the sense of the renormalization group, curves that represent the dependence of a Binder ratio on a control parameter ($z$ or $T$) for systems of different sizes $L$ are all expected to cross at the critical value of the control parameter. This allow us to locate the phase transitions involving loss of order in a convenient way.

The correlation length is obtained as

where $\lambda_{0}$ ($\lambda_{1}$) is the (second) largest eigenvalue of the transfer matrix defined by the uMPS. We note that the uMPS is always normalized such that $\lambda_{0}=1$, and the factor $2$ comes from the unit-cell size of the uMPS. Since we calculate the uMPS in both the horizontal and vertical directions, two kinds of correlation length exist in a phase that breaks the symmetry of lattice rotations. The “transverse” correlation length $\xi_{\perp}$ is measured along the direction perpendicular to the dimers, while the “longitudinal” correlation length $\xi_{\parallel}$ is in the direction parallel to the dimers. In the disordered fluid phase, both correlation lengths are the same as expected. In the nematic phase, a level crossing of $\lambda_{1}$ occurs, and $\xi_{\perp}$ takes a smaller value than $\xi_{\parallel}$ below a certain temperature. In the columnar ordered phase, we always have $\xi_{\perp}<\xi_{\parallel}$. In other words, the correlation length scale along the dimers become larger than in the perpendicular direction. Thus, the transverse correlation length seems to be suitable for studying the phase transition between the columnar ordered and disordered fluid phases. We also note that $\xi$ corresponds to the truncated correlation function and takes a finite value even in the ordered phase.

The entanglement entropy is defined as

where $\sigma_{i}$ denotes the singular value of the core matrix in the mixed canonical form of the uMPS. Since we calculate the horizontal and vertical uMPS with a $2\times 2$ unit-cell structure, $S_{\text{EE}}$ may depend on direction and position in the unit cell. In the disordered and nematic phases, we find that the all $S_{\text{EE}}$ are equal to each other. On the other hand, in the columnar ordered phase, we find that $S_{\text{EE}}$ takes on two values, depending on whether the core matrix of the horizontal uMPS is on a dimer or between dimers. The former always yields the larger $S_{\text{EE}}$. In our analysis below, we use the largest entanglement entropy thus obtained.

It is instructive to think in terms of a coarse-grained version $\psi$ of our local complex columnar order parameter field $\Psi_{\text{col}}$, and write

Clearly, the corresponding coarse-grained version $\phi$ of the local nematic order parameter field $\Psi_{\text{nem}}$ satisfies

The $\tau$ defined in the above are two coarse-grained Ising fields. In the columnar ordered phase, both $\tau_{1}$ and $\tau_{2}$ are ordered; this correctly accounts for the four-fold symmetry breaking in the columnar ordered state. Nematic order corresponds to the product $\tau_{1}\tau_{2}$ being ordered, without any long range order in the individual $\tau$. From the symmetries of the original problem, we see that interchanging the $\tau$ is a symmetry of the theory. Thus, the natural description is in terms of a symmetric Ashkin-Teller theory with two Ising fields $\tau_{1}$ and $\tau_{2}$.

This connection to the physics of the Ashkin-Teller model [19, 20, 21] yields a wealth of information. For instance, along a line of continuous transitions from the columnar ordered state to the disordered fluid state, we expect the critical behavior to controlled by the critical properties of the corresponding fixed line in the Ashkin-Teller model. Along this line, both the Ising fields $\tau_{1}$ and $\tau_{2}$ have a fixed anomalous dimension of $\eta=1/4$ [19, 20, 21]. Since the columnar order parameter is linear in the Ising fields $\tau_{1/2}$, we expect it to also scale with an anomalous exponent $\eta=1/4$ all along the line of continuous transitions between columnar ordered and disordered fluid phases.

Along the fixed line of the Ashkin-Teller theory, $\tau_{1}\tau_{2}$ scales with an anomalous dimension $\eta_{2}$ that varies continuously [19, 20, 21] and is related to the continuously varying correlation length exponent by the Ashkin-Teller relation $\eta_{2}=1-1/2\nu$. Since the nematic order parameter $\phi\sim\tau_{1}\tau_{2}$, we expect it to have an anomalous dimension $\eta_{2}$ [7] given by this relation all along the line of continuous transitions between columnar ordered and disordered fluid phases.

In the Ashkin-Teller model, the point at which Ashkin-Teller line splits into two lines of Ising transitions is known to have the symmetries for the four-state Potts model [19, 20, 21, 25] ; in the phase between these two Ising transition lines, $\tau_{1}\tau_{2}$ is ordered although $\tau_{1}$ and $\tau_{2}$ remain individually disordered . The enhanced Potts symmetry at this multicritical point implies that $\tau_{1}$, $\tau_{2}$, and $\tau_{1}\tau_{2}$ all have the same anomalous exponent. Thus $\eta_{2}=\eta=1/4$ at this point, and the Ashkin-Teller relation implies $\nu=2/3$. Given the correspondence made above, this implies that the nematic order parameter is expected to have an anomalous exponent of $1/4$ at the multicritical point at which the Ising phase boundaries of the nematic phase meet the line of continuous transitions between columnar ordered and disordered fluid phases.

From this perspective, it is clear that the value of $\nu$ (or equivalently $\eta_{2}$) serves as a universal coordinate for the line of continuous transitions from the columnar ordered phase to the disordered fluid phase. At full-packing, i.e., $z=0$, the system has a description in terms of a coarse-grained Gaussian height action for a scalar height $h$, and the transition from the power-law ordered high temperature state to the low temperature state is expected to be governed by a Kosterlitz-Thouless transition at which the leading cosine nonlinearity $\cos(8\pi h)$ becomes relevant. As a result, one expects $\nu\to\infty$ as $z\to 0$ along the line of continuous transitions between the columnar ordered state and the disordered fluid state. From this it is clear that the multicritical point at which the two Ising lines meet cannot be at $z=0$, since this multicritical point corresponds to a value of $\nu=2/3$.

This theoretical perspective and the resulting expectations informs much of the data analysis we present in the next section.

To this end, we first consider a fixed $z=0.2$ slice and display the computed two-dimensional phase diagram in the $T$–$v$ plane (with $u=1-v$). This is shown in Fig. 4.

For $v\leq 0.7$, there is a direct temperature driven phase transition between the columnar ordered and disordered fluid phases. As $v$ is increased further, this phase boundary splits slightly below $v=0.8$ into two transition lines and the nematic phase appears as an intermediate phase beyond this multicritical point at which three transition lines meet. We have also checked that a corresponding slice at somewhat larger $z$ reduces the temperature scale of the transitions. The transition temperatures shown in Fig. 4 is determined from the peak position of the correlation length estimated by the TN method. The result of the MC method agrees with it within errors that are smaller than the symbol sizes used.

Next we consider slices with fixed $v=0.9$ and $v=0.8$ (with $u=1-v$) and display the computed two-dimensional phase diagram in the $T$–$z$ plane. This is shown in Fig. 5. The transition points are determined in the same way as Fig. 4. Since the fully packed $z=0$ system is particularly challenging for TN computations, we do not extend our study all the way to $z=0$. Nevertheless, we are able to go to low enough $z$ to demarcate the essential features of the phase diagram. It would be difficult to obtain this phase diagram using the MC method because the phase boundaries cannot be calculated with such precision due to the finite size effect.

At $v=0.8$, the low-temperature nematic phase is quite narrow and disappears below the multicritical value of $z$ which is close to $z=0.1$. On the other hand, for $v=0.9$, the nematic phase is very broad and seems to exist even at very low values of $z$. However, our detailed computations reveal that there is no nematic state below a multicritical threshold value of $z$ which is close to $z=0.002$ (Fig. 5). The actual monomer densities associated with these multicritical points are extremely small: For $v=0.8$, the multicritical monomer density is about $\delta=0.005$. At $v=0.9$, the multicritical monomer density is a much lower value of $5\times 10^{-5}$.

This conclusion is contrary to that of Ref. [22]. However, it is entirely consistent with our understanding of the phase boundaries based on the coarse-grained effective field theory. Indeed, as we have already reviewed earlier, the multicritical point is expected to have rather different universal behavior from the transition at full packing since the former is expected to correspond to a $\eta_{2}=1/4$ and the latter corresponds to $\eta_{2}=1$ (where $\eta_{2}$ is the anomalous exponent associated with the nematic order parameter). As a result, for $v/u<\infty$, there is no consistent scenario in which the multicritical point coincides with the transition at full packing. The resolution is of course that the multicritical value of $z$ approaches $z=0$ extremely rapidly as $v/u$ is increased, but does not reach $z=0$ at any finite $v/u$.

In the rest of this section, we display our results for a few representative values of $v$ and $z$, and provide a detailed account of the analysis that leads to these phase diagrams and this overall conclusion.

We use both tensor network (TN) and Monte Carlo (MC) methods to obtain the nematic and columnar order parameters, since these two complementary methods provide a nontrivial consistency check on each other. Figure 6 shows the temperature dependence of the order parameters at $v=0.6$ and $v=0.9$ with $z=0.2$.

At $v=0.6$, both the order parameters take on nonzero values below $T=0.7714$, signaling the onset of columnar order. On the other hand, at $v=0.9$, we find that a nematic phase exists between $T_{1}=0.722$ and $T_{2}=0.790$ , as is clear from the fact that the columnar order parameter vanishes but the nematic order parameter takes on a nonzero value. The transition temperatures are estimated by identifying the crossing points of the Binder ratios, as shown in Fig. 7 for the Binder ratios calculated by the MC simulations at $v=0.9$ and $z=0.2$. These crossing points are found to be consistent with the transition temperature estimated by the TN method.

We obtain the scaling dimensions of various quantities at these transitions by performing a finite-size scaling (FSS) analysis, using the following FSS form for the columnar and nematic order parameters,

where $x_{a}$ denotes the scaling dimension of the corresponding operator $\Psi_{a}$ ($a=\text{col},\text{nem}$); these scaling dimensions are related to the anomalous exponents introduced earlier via $x_{\text{col}}=\eta/2$, $x_{\text{nem}}=\eta_{2}/2$. Likewise, the scaling dimension of the energy operator is related to the correlation length exponent introduced earlier: $x_{t}=2-1/\nu$.

We use the kernel method [36] to infer the critical exponents and the critical temperature and estimate their confidential intervals. At the best fit value of the scaling dimensions and the transition point, all data collapse reasonably well onto a single curve as shown in Fig. 8.

The scaling dimensions are plotted in Fig. 9. Below the multicritical point, i.e., along the line of continuous phase transitions between the columnar ordered and disordered fluid states, the scaling dimensions $x_{\text{nem}}$ and $x_{t}$ are seen to continuously change with $v$, while $x_{\text{col}}$ remains constant at a value consistent with the theoretical expectation of $x_{\text{col}}=1/8$. We have checked that the corresponding exponents $\eta_{2}$ and $\nu$ satisfy the Ashkin-Teller relation $\eta_{2}=1-1/2\nu$, i.e., $x_{t}/4=x_{\text{nem}}$ within our errors. On the other hand, we have $x_{t}/4\neq x_{\text{nem}}$ for $v\geq 0.8$. This discontinuous change in the critical index ratio strongly suggests the change from the single transition to the two separate transitions, i.e., the existence of the multi-critical point. The expected $\eta=\eta_{2}=1/4$ at the multicritical point is actually realized at about $v=0.7$. Parenthetically, we note that we find a decoupled Ising point at about $v=0.3$ in the $z=0.2$ plane, where we have $x_{\text{nem}}=0.25$ and $x_{t}=1.0$.

The scaling dimension of the energy operator can also be estimated from the FSS analysis of the Binder ratio as shown in Fig. 8. Although the nonmonotonic behavior of the Binder ratio for columnar order and the relatively closely spaced successive transitions makes this difficult, we obtain almost the same results for $x_{t}$ from this analysis, as we do using the earlier analysis in terms of the order parameters. We emphasize that nonmonotonicity has to do with the differing character of the columnar order parameter fluctuations in the nematic and the disordered fluid phases, and the related presence of a proximate multicritical point. Similar nonmonotonicity has been noted in the $J_{1}$-$J_{2}$ Ising model earlier [37]. In contrast to other examples of such behavior in frustrated Ising models, which is associated with a proximate weakly first order transition, we do not find evidence of any first-order transitions for the values of $v$ studied here.

Although it is difficult to completely exclude possibility of the weakly first-order phase transition around the multicritical point in our simulations, the weight of evidence suggests that the presence of a sizable $v$ term replaces the first order transition found in Refs. [16, 17] by an intermediate nematic phase flanked by two second-order Ising phase boundaries. This is consistent with the fact that generalized four-state clock models, which serve as a discrete hard-spin analog of the coarse-grained description in terms of order parameter fields $\psi$ and $\phi$, are known for some parameter values to have such an intermediate phase flanked by two Ising lines that meet at a multicritical point with enhanced four-state Potts symmetry at the end of a line of Ashkin-Teller transitions [25].

The multicritical point at $v=v_{c}(z)$ is expected to have the four-state Potts universality. It is difficult to determine accurate location of the multicritical point because of the finite-size or finite bond-dimension effect. To make matters worse, one also expects that a logarithmic correction appears at the four-state Potts model. Our data, however, support the $S_{4}$ symmetry at the multicritical point. The order parameters satisfy $m_{\text{nem}}>m_{\text{col}}$ below $v_{c}$, while $m_{\text{nem}}>m_{\text{col}}$ for $v>v_{c}$ (Fig. 6). Thus we expect that $m_{\text{col}}\simeq m_{\text{nem}}$ at $v=v_{c}$, which indicates the emergent $S_{4}$ symmetry.

By definition, the scaling dimension also appears in the corresponding critical two-point correlation function as $C_{a}(r)\propto r^{-2x_{a}}$. We consider the correlation function along an axis because the uMPS can easily calculate it. The correlation function between the local quantities is defined as

In our TN simulation, $\bm{\rho}$ runs over sites in the $2\times 2$ unit cell, which corresponds to a value of $N=4$. For correlation between monomers, the truncated correlation function,

is expected to scale as $r^{-2x_{t}}$. As shown in Fig. 10, the scaling dimensions extracted by the linear fitting of the correlation functions agree with the results obtained by the FSS analysis.

The central charge is another important universal property of a critical point. According to the conformal field theory, the correlation length $\xi$ and the entanglement entropy $S_{\text{EE}}$ are related by the Calabrese-Cardy formula at criticality:

where $c$ denotes the central charge.

One of the advantages of the TN method is that these two quantities can be calculated naturally. Figure 11 clearly shows that the values obtained from the TN method are consistent with the Cardy-Calabrese formula (21). Based on the theoretical framework outlined in the previous section, the continuous phase transition between the columnar ordered and disordered fluid phases is expected to have a central charge of $c=1$. On the other hand, the continuous Ising phase boundaries of the nematic phase are expected to have a central charge of $c=1/2$. From the results displayed in Fig. 11, we see that our results do indeed conform to both these expectations.

Above the multicritical point ($v>v_{c}$), there are two phase transitions, the columnar-nematic and nematic-disorder ones. Although the FSS result of $x_{t}$ deviates from the expected value (Fig. 9), we believe that it is due to the correction to scaling and effects from another critical point. The TN result at $v=0.9$ shows that the eighth power of the order parameters and the inverse of the correlation length are linear to the temperature near the criticality (Fig. 12). This fact strongly indicates that these critical exponents satisfy $\beta=1/8$ and $\nu=1$, which is consistent with the Ising universality.

Finally, we comment on the approach of the multicritical point to $z=0$ as $v/u$ is increased. This happens very rapidly, and simultaneously, the phase boundary between the columnar and nematic phases, shown by orange squares in Fig. 5, approaches the vertical temperature axis as $v/u$ increases. This is clear from monitoring the relative values of $m_{\text{col}}$ and $m_{\text{nem}}$ as in our earlier discussion about the symmetry of the multicritical point. Eventually, in the limit $v/u\rightarrow\infty$, the columnar phase vanishes and only the phase boundary between the nematic and disordered phases remains.