Hierarchy of the low-lying excitations for the $(2+1)$-dimensional $q=3$ Potts model in the ordered phase

The $(2+1)$-dimensional $q=3$ Potts model was simulated with the exact diagonalization method. In the ordered phase, the elementary excitations (magnons) are attractive, forming a series of bound states in the low-energy spectrum. We investigate the low-lying spectrum through a dynamical susceptibility, which is readily tractable with the exact diagonalization method via the continued-fraction expansion. As a result, we estimate the series of (scaled) mass gaps, $m_{2,3,4}/m_1$ ($m_1$: single-magnon mass), in proximity to the transition point.

According to the Monte Carlo simulations in (2 + 1) dimensions [18,19], the hierarchy m 1,2,... of the Potts model and that of the pure gauge theory look alike.
Actually, for the Z 2 case in (2 + 1) dimensions, a duality relation [20,21,22] does hold, ensuring the correspondence between them. Generically [23], the SU(N ) gauge theory displays a global Z N symmetry (center of SU(N )), which immediately establishes a relationship between them; for N ≥ 3, the transition would not be critical, and the universality idea does not apply nonetheless.
Meanwhile, an extensive lattice-gauge-theory simulation reveals the "weak N dependence" [24] of the SU(N ) theory, suggesting a robustness of the hierarchy, m 2,3,... /m 1 . On the one hand, the q = 3 Potts model exhibits an "approximate universality" [25], even though the phase transition is of first order. Hence, it is expected that the hierarchy should display a model-independent character to some extent.
In this paper, we investigate the (2 + 1)-dimensional q = 3 Potts model [26,27] by means of the numerical diagonalization method. The method allows us to evaluate the dynamical susceptibilities, Eqs. (12) and (17), via the continuedfraction-expansion method [28]. In Fig. 1, we present a schematic drawing for a spectral function.
To be specific, we present the Hamiltonian for the (2 + 1)-dimensional q = 3 Potts model Here, the operator is placed at each square-lattice point i = 1, 2, . . . , N ; namely, the base |l i (l i = 0, 1, 2) satisfies L i |l i = l i |l i . The summation ij runs over all possible nearest neighbor pairs ij . The operator R ± i induces the transition R ± i |l i = |l i ±1 mod 3 , and the parameter λ denotes the corresponding coupling constant.
The rest of this paper is organized as follows. In the next section, we present the numerical results. The simulation algorithm is presented as well. In Sec. 3, we address the summary and discussions.

Numerical results
In this section, we present the numerical results for the Potts model (8). To begin with, we explain the simulation algorithm.

Numerical algorithm
We employed the exact diagonalization method to simulate the Potts model (8) for a rectangular cluster with N ≤ 22 spins. In order to treat a variety of N = 16, 18, . . . systematically, we implemented the screw-boundary condition [33]. According to Ref. [33], an alignment of spins {l i } with the first-and √ N -th-neighbor interactions reduces to a rectangular cluster under the screwboundary condition. Based on this idea, we express the Hamiltonian as Here, the diagonal matrix H D (v) denotes the v-th-neighbor interaction for an The above formulae are mathematically closed; however, for an efficient simulation, Eqs. (9) and (10) of Ref. [34] may be of use.
We performed the numerical diagonalization for the Hamiltonian matrix (10) by means of the Lanczos method. The single-magnon mass gap m 1 is given by with the ground-state (E 0 ) and first-excited (E 1 ) energy levels within the zeromomentum sector. Because the N spins form a rectangular cluster, the linear dimension is given by L = √ N , which sets a foundermental length scale in the subsequent scaling analyses.

Single-magnon mass gap m 1
In this section, we investigate the single-magnon mass gap m 1 (11) with the scaling theory [26,31]. The first-order phase transition obeys the properly remedied scaling theory.
The data in Fig. 2 seem to collapse into a scaling curve, indicating that the simulation data already enter the scaling regime. Encouraged by this finding, we turn to the analysis of the spectral properties.

Hierarchical spectral peaks
Based on the finite-size scaling [26,31] demonstrated in the preceding section, we analyze the dynamical susceptibility with the ground-state energy (vector) E 0 (|0 ) and the energy-resolution parameter η. Here, the perturbation operator is set to with and the projection operator P = 1 − |0 0|. We calculated the dynamical susceptibility (12) with the continued-fraction expansion [28]. The dynamical susceptibility (spectral function) obeys the scaling formula with a certain scaling function f [35,12].
Each signal in Fig. 3 is interpreted by the diagram in Fig. 1. That is, the peaks around ω/m 1 ≈ 1.8, 2.5 and 3 correspond to the m 2,3,4 excitations, respectively. The shoulder peak around ω/m 1 ≈ 2 should be the two-magnonspectrum threshold. The signal ω/m 1 ≈ 3.5 may be either the m 5 particle or a composite one consisting of m 1 and m 2 . The ratios m 2,3,4 /m 1 are estimated in the next section more in detail.
Last, we mention the choice of the perturbation operator A (13). In a preliminary stage, we surveyed various types of the perturbation operators, aiming to create the m 2,3,... particles effectively. Actually, neither the first or the second term of the Hamiltonian (8) commutes with A; otherwise, the susceptibility reduces to a mere specific heat. A key ingredient is that the exact diagonalization method permits us to treat any off-diagonal operators. Last, we address a remark. Because the phase transition is discontinuous, the continuum limit cannot be taken properly. The above estimates such as the life time are specific to a lattice realization, although seemingly preferable scaling behavior was observed. However, in an approximate sense, the simulation results seem to be comparable with the related ones, as claimed by the preceeding studies [18,19,25].

Continuum-threshold peak via χ ′′ B (ω)
As a comparison, we present a simulation result for χ ′′ B (17), aiming to see to what extent the spectral weight is affected by the choice of the perturbation operator. The dynamical susceptibility χ ′′ B is defined by with a perturbation operator Note that the operator B coincides with the second term of the Hamiltonian (8).
The result indicates that a naive external disturbance such as the specificheat-type perturbation B does not create the bound states higher than m 1,2 very efficiently. It is significant to set up the perturbation operator, which does not commute with any terms of the Hamiltonian. Note that the exact diagonalization method allows us to survey various types of the (off-diagonal) perturbation operators so as to observe m 3,4,... clearly.
As a reference, we calculated χ ′′ B (17); here, the operator B coincides with the second term of the Hamiltonian (8), and hence, it would be relevant to the experimental study. It turned out that the probe χ ′′ B is insensitive to the hierarchy m 3,4,... , indicating that the choice of the perturbation operator is vital to observe m 3,4,... . In this sense, the exact diagonalization method has an advantage in that we are able to treat various perturbation operators so as to observe the hierarchy m 3,4,... clearly. here, we adopted the scaling theory [26,31] for the first-order phase transition.
[25] and [26,31], respectively; namely, there is no adjustable parameter involved in the scaling analysis.  λc and ν, are the same as those of Fig. 2. A faint signal around ω/m 1 ≈ 2 is attributed to the two-magnon-spectrum threshold. The spectral weights for χ ′′ B differ significantly from those for χ ′′ A .