Multi-particle states inverstigation with tensor renormalization group method

Let us start with a brief review of the spectroscopy scheme using transfer matrix and tensor network that we proposed in Az-zahra et al. (2024). The scheme is explained in the framework of the two dimensional scalar field theory with nearest-neighbor interactions on lattice. However, its extension to higher dimensional systems is straightforward, and in principle, the scheme is also applicable to fermionic or gauge systems. The lattice action of the (1+1)d scalar field theory in Euclidean space-time is given by

$$S[\phi]=\sum_{{\bm{r}}\in\Gamma}\left[\sum_{\mu=0}^{1}\frac{1}{2}\left(\phi({\bm{r}}+\hat{\mu})-\phi({\bm{r}})\right)^{2}+V(\phi({\bm{r}}))\right].$$

(1)

Here, $\phi(\bm{r})$ is the scalar field resides on two dimensional square lattice ${\bm{r}}=(t,x)\in\Gamma$ and $\hat{\mu}$ is the unit vector for the $\mu=0,1$ direction, where the $0$ direction ($1$ direction) is considered as the time (space) direction. The lattice $\Gamma$ has periodic boundary in spatial direction, and can be defined as

$$\Gamma=\{(t,x)|t=0,1,2,\ldots,L_{\rm t}-1\,\,{\rm and}\,\,x=0,1,2,\ldots,L_{\rm s}-1\},$$

(2)

where $L_{\rm t}$ and $L_{\rm s}$ denote the lattice size in time and space direction, respectively. Note that the mass term and self-interaction term are already included in the potential $V(\phi({\bm{r}}))$. Accordingly, the partition function of the system is given by

$$Z=\int\prod_{{\bm{r}}\in\Gamma}d\phi({\bm{r}})e^{-S[\phi]}.$$

(3)

This partition function can be reformulated in terms of transfer matrix

$$Z={\rm Tr}\left[{\cal T}^{L_{\rm t}}\right],$$

(4)

where the transfer matrix ${\cal T}$ (see Fig. 1(a) for the diagram) is given by Montvay and Münster (1994)

	$\displaystyle{\cal T}_{\Phi^{\prime}\Phi}$	$\displaystyle=$	$\displaystyle\left(\prod_{x=0}^{L_{\rm s}-1}\exp\left[-\frac{1}{2}(\phi_{x}^{\prime}-\phi_{x})^{2}-\frac{1}{4}V(\phi^{\prime}_{x})-\frac{1}{4}V(\phi_{x})\right]\right)$		(5)
		$\displaystyle\times$	$\displaystyle\left(\prod_{x=0}^{L_{\rm s}-1}\exp\left[-\frac{1}{4}(\phi_{x+1}^{\prime}-\phi_{x}^{\prime})^{2}-\frac{1}{8}V(\phi^{\prime}_{x+1})-\frac{1}{8}V(\phi^{\prime}_{x})\right]\right)$
		$\displaystyle\times$	$\displaystyle\left(\prod_{x=0}^{L_{\rm s}-1}\exp\left[-\frac{1}{4}(\phi_{x+1}-\phi_{x})^{2}-\frac{1}{8}V(\phi_{x+1})-\frac{1}{8}V(\phi_{x})\right]\right),$

where

	$\displaystyle\Phi^{\prime}$	$\displaystyle=$	$\displaystyle\{\phi^{\prime}_{x}=\phi(t+1$	,	$\displaystyle x)$	$\displaystyle|$	$\displaystyle x=0,1,2,\ldots,L_{\rm s}-1\},$		(6)
	$\displaystyle\Phi$	$\displaystyle=$	$\displaystyle\{\phi_{x}=\phi(t$	,	$\displaystyle x)$	$\displaystyle|$	$\displaystyle x=0,1,2,\ldots,L_{\rm s}-1\}$		(7)

are field configurations on the Euclidean time slice at $t+1$ and $t$. In this case, ${\cal T}_{\Phi^{\prime}\Phi}$ is treated as a usual matrix where $\Phi$ is treated as integer-valued index for notational convenience. The diagonalization of ${\cal T}_{\Phi^{\prime}\Phi}$ is given by

$${\cal T}_{\Phi^{\prime}\Phi}=\sum_{a=0}^{\infty}U_{\Phi^{\prime}a}\lambda_{a}(U^{\dagger})_{a\Phi}\,\,\Longleftrightarrow\,\,\hat{\cal T}|a\rangle=\lambda_{a}|a\rangle.$$

(8)

Here $U_{\Phi^{\prime}a}=\langle\Phi^{\prime}|a\rangle$ is the field representation of the eigenstate $|a\rangle$. The eigenvalues $\lambda_{a}$ give us the information of the energy eigenvalues $E_{a}$ of the system

$$\lambda_{a}=e^{-E_{a}}\hskip 28.45274pt\mbox{ for }a=0,1,2,3,\ldots.$$

(9)

where $\lambda_{0}\geq\lambda_{1}\geq\ldots$. Instead of energy $E_{a}$, the energy gaps $\omega_{a}$

$$\omega_{a}=E_{a}-E_{0}\hskip 28.45274pt\mbox{ for }a=1,2,3,\ldots$$

(10)

where $E_{0}$ is the ground state energy, are more useful. So that, hereafter, the energy gap spectrum $\omega_{a}$ will be mentioned as energy spectrum for simplicity.

The quantum number of each eigenstate $|a\rangle$ is not a priori known. Therefore, we identify it using matrix elements

$$\langle b|\hat{\mathcal{O}}_{q}|a\rangle=\left(U^{\dagger}{\mathcal{O}_{q}}U\right)_{ba},$$

(11)

where $U$ is the unitary matrix from Eq. (8), $\hat{\mathcal{O}}_{q}$ is an interpolating operator with quantum number $q$ and $\mathcal{O}_{q}$ is the field representation of the interpolating operator

$$(\mathcal{O}_{q})_{\Phi^{\prime}\Phi}=\langle\Phi^{\prime}|\hat{\mathcal{O}}_{q}|\Phi\rangle.$$

(12)

A selection rule derived from the symmetry determines the quantum number of the eigenstate $|a\rangle$. As derived in Az-zahra et al. (2024), the selection rule for the system with discrete symmetry is given by

$$\langle b|\hat{\mathcal{O}}_{q}|a\rangle\neq 0\rightarrow q_{b}qq_{a}=1,$$

(13)

where $q_{a},q_{b}$ are the quantum numbers for eigenstates $|a\rangle$ and $|b\rangle$, respectively. Solving Eq. (13) for known $q_{b}$ and $q$, yields the quantum number of the state $|a\rangle$. The quantum number in systems with continuous symmetries can also be classified using a similar approach, see Az-zahra et al. (2024).

After presenting the formalism for computing the energy spectrum and identifying the corresponding quantum numbers, we now discuss how to compute them numerically. Because the dimensionality of transfer matrix is extremely large, we apply an approximation method to compute its eigenvalues and extract the energy spectrum, which in our case is based on tensor network techniques. To construct the tensor network representation of the transfer matrix in Eq. (5), we first decompose the transfer matrix into

$${\cal T}_{\Phi^{\prime}\Phi}=\sum_{k}Y_{\Phi^{\prime}k}Y^{\dagger}_{k\Phi},$$

(14)

where $k=\{k_{x}|x=0,1,2,\ldots,L_{\rm s}-1\}$ is a newly defined integrated index and the matrix $Y$ is defined as

	$\displaystyle Y_{\Phi k}$	$\displaystyle=$	$\displaystyle\left(\prod_{x=0}^{L_{\rm s}-1}\sum_{k_{x}=0}^{\infty}u_{\phi_{x}{k_{x}}}\sqrt{\sigma_{k_{x}}}\right)$		(15)
		$\displaystyle\times$	$\displaystyle\left(\prod_{x=0}^{L_{\rm s}-1}\exp\left[-\frac{1}{4}(\phi_{x+1}-\phi_{x})^{2}-\frac{1}{8}V(\phi_{x+1})-\frac{1}{8}V(\phi_{x})\right]\right).$

The matrix $u_{\phi_{x}k_{x}}$ and $\sigma_{k_{x}}$ are obtained from the eigenvalue decomposition (EVD) of the Boltzmann weight in Eq. (5) for the fields $(\phi,\phi^{\prime}\in\mathbb{R})$, specifically from the temporal hopping term, namely

$$\exp\left[-\frac{1}{2}(\phi^{\prime}-\phi)^{2}-\frac{1}{4}V(\phi^{\prime})-\frac{1}{4}V(\phi)\right]=\sum_{k=0}^{\infty}u_{\phi^{\prime}{k}}\sigma_{k}(u^{\dagger})_{k\phi}.$$

(16)

Using the matrix $Y$ in Eq. (15), the partition function in Eq. (4) can be rewritten as

$$Z=\operatorname{Tr}_{\Phi}{\left[{\cal T}^{L_{\rm t}}\right]}=\operatorname{Tr}_{\Phi}{\left[(YY^{\dagger})^{L_{\rm t}}\right]}=\operatorname{Tr}_{k}{\left[(Y^{\dagger}Y)^{L_{\rm t}}\right]}=\operatorname{Tr}_{k}{\left[{\cal A}^{L_{\rm t}}\right]}.$$

(17)

Here, we have defined

$${\cal A}\coloneqq Y^{\dagger}Y.$$

(18)

As seen in Fig. 1(b), $\mathcal{A}$ can be written as a product of rank-four tensors

$${\cal A}_{kj}=\sum_{\{l\}}\prod_{x=0}^{L_{\rm s}-1}A_{k_{x}l_{x}j_{x}l_{x-1}}$$

(19)

where the tensor $A_{k_{x}l_{x}j_{x}l_{x-1}}$ is computed by

$$A_{k_{x}l_{x}j_{x}l_{x-1}}\coloneqq\sqrt{\sigma_{k_{x}}\sigma_{l_{x}}\sigma_{j_{x}}\sigma_{l_{x-1}}}\sum_{\phi_{x}}\left(u^{\dagger}\right)_{k_{x}\phi_{x}}\left(u^{\dagger}\right)_{l_{x}\phi_{x}}u_{\phi_{x}j_{x}}u_{\phi_{x}l_{x-1}}.$$

(20)

The EVD of this single-time slice tensor network $\mathcal{A}_{kj}$ is given by

$${\cal A}_{kj}=W_{ka}\lambda_{a}W^{\dagger}_{aj},$$

(21)

where $\lambda_{a}$ are the same eigenvalues as in Eq. (8) and $W$ is a diagonalization matrix.

Next, we consider the matrix elements $\langle b|\hat{\mathcal{O}}_{q}|a\rangle$. This matrix also has very large dimensionality, so we apply the tensor network methods to approximate it as well. For this purpose, we first rewrite $\langle b|\hat{\mathcal{O}}_{q}|a\rangle$ in the tensor network language. The complete procedure for formulating matrix elements in terms of tensor network is given in Az-zahra et al. (2024), where the final expression is given by

$$\langle b|\hat{\mathcal{O}}_{q}|a\rangle=\left(\lambda^{-m+1/2}W^{\dagger}\mathcal{A}^{m-1}\mathcal{A}^{\prime}\mathcal{A}^{m}W\lambda^{-m-1/2}\right)_{ba}.$$

(22)

Here, ${\cal A}^{m-1}{\cal A}^{\prime}{\cal A}^{m}$ represents an impurity tensor network, and $1\leq m\leq L_{\rm t}/2$ specifies the position of $\cal A^{\prime}$ within the network (see Fig. 4). Note that, ${\cal A}^{\prime}$ is the impurity version of $\mathcal{A}$, that is $\mathcal{A}^{\prime}\coloneqq Y^{\dagger}\mathcal{O}_{q}Y$, where $Y$ is the matrix given in Eq. (15). For a single field $\mathcal{O}_{q}=\phi_{x}$ at lattice site $x$, an impurity tensor $A^{\prime}$ is defined as

$$A_{k_{x}l_{x}j_{x}l_{x-1}}^{\prime}\coloneqq\sqrt{\sigma_{k_{x}}\sigma_{l_{x}}\sigma_{j_{x}}\sigma_{l_{x-1}}}\sum_{\phi_{x}}\phi_{x}(u^{\dagger})_{k_{x}\phi_{x}}(u^{\dagger})_{l_{x}\phi_{x}}u_{\phi_{x}j_{x}}u_{\phi_{x}l_{x-1}}.$$

(23)

In Fig. 1(c) we show an image of $\mathcal{A}^{\prime}$, which is composed of several pure tensors $A$ and the single impurity tensor $A^{\prime}$ located at $x=0$. The mathematical expression of the single-time slice impurity tensor network $\mathcal{A}^{\prime}_{kj}$ shown in terms of $A$ and $A^{\prime}$ is given by

$$\mathcal{A}^{\prime}_{kj}=\sum_{\{l\}}A^{\prime}_{k_{0}l_{0}j_{0}l_{L_{\rm s}-1}}\prod_{x=1}^{L_{\rm s}-1}A_{k_{x}l_{x}j_{x}l_{x-1}}.$$

(24)

Here, we explain how to coarse grain the tensor network representation of transfer matrix and matrix elements derived in the previous section using higher order tensor renormalization group (HOTRG) algorithm Xie et al. (2012). In the following, the bond dimension of the initial tensor is denoted by $\chi$.

It is well-known that single coarse-graining iteration of HOTRG naturally reduces the tensor network size by a factor of two, thereby restricting the accessible system sizes to $L_{\rm s}=2^{n}$, where $n=0,1,2,\ldots$ is the number of coarse-graining step. To access a wider range of spatial sizes, we prepare an initial single-time slice tensor network ¹¹1The idea of constructing the initial tensor network for transfer matrix was first proposed in Ref. Huang et al. (2023). with spatial size $L_{\rm s}^{(\rm init)}$ and then embed it into the main tensor network such that the total spatial size becomes $L_{\rm s}=L_{\rm s}^{(\rm init)}\times 2^{n}$. Note that, throughout this work, we use initial tensor sizes $L_{\rm s}^{(\rm init)}\in\{1,3,5,7\}$.

For this purpose, we denote the tensor $A$ in Eq. (20) as a bare tensor $A^{[0,0]}$. The size of initial tensor network is set as $L_{\rm s}^{(\rm init)}=r+1$ for $r\in\{0,2,4,6\}$. As shown in Fig. 2(a), contractions are applied to the bare tensors $A^{[0,0]}$ using HOTRG-like algorithm repeatedly until only a single tensor $A^{[r,0]}$ remains. We refer to this as HOTRG-like algorithm because the coarse-graining is applied to two different tensors (see Fig. 2), although the procedure is essentially the same as in the standard HOTRG algorithm. To be concrete, a single coarse-graining step that produces the new tensor $A^{[s,0]}$ from the previously coarse-grained tensor $A^{[s-1,0]}$ is given by

$\displaystyle A^{[s,0]}_{c_{0}a_{1}b_{0}a_{0}}$

$\displaystyle=\vbox{\hbox{\includegraphics[height=85.35826pt]{picture/init_contraction_pure_eqn.png}}}\hskip 28.45274pt(s=1,2,\ldots,r).$

(25)

For $s=1$, $A^{[s-1,0]}$ corresponds to the bare tensor. The tensor $U^{[s,0]}$ in Eq. (25) is an isometry that is computed from the EVD of the following self-adjoint matrix

$$M_{(k_{0}k_{1},k_{0}^{\prime}k_{1}^{\prime})}^{[s,0]}=\vbox{\hbox{\includegraphics[height=85.35826pt]{picture/init_contraction_M_eqn.png}}}.$$

(26)

Applying EVD and truncating the eigenvalues of $M_{(k_{0}k_{1},k_{0}^{\prime}k_{1}^{\prime})}^{[s,0]}$ yield

$$M_{(k_{0}k_{1},k_{0}^{\prime}k_{1}^{\prime})}^{[s,0]}\approx\sum_{c_{0}=0}^{\chi}\vbox{\hbox{\includegraphics[height=56.9055pt]{picture/decomposition.png}}}$$

(27)

where $U^{[s,0]}$ is the isometry and $\sigma_{c_{0}}$ are the eigenvalues of $M$ that are truncated up to the cut-off bond dimension $\chi$. See Xie et al. (2012) for a complete description of the isometry.

After completing the coarse-graining of the initial tensor network with size $L_{\rm s}^{(\rm init)}$, we obtain the tensor $A^{[r,0]}$. As mentioned earlier, we embed this tensor into the main tensor network, which then consists of $2^{p}\times 2^{n}$ tensors $A^{[r,0]}$, as illustrated in Fig. 3. Accordingly, the total size of the main tensor network is $L_{\rm t}\times L_{\rm s}=2^{p}\times(L_{\rm s}^{(\rm init)}\times 2^{n})$. The coarse-graining of this network is performed using the standard HOTRG algorithm Xie et al. (2012), where the new tensor $A^{[r,i]}$ for $i=1,2,\ldots,n$ is obtained from previously coarse-grained tensors $A^{[r,i-1]}$ by coarse-graining for the time and spaces direction.

An important point to note is that, $p$ and $n$ do not only determine the size of tensor network, but also represent the number of coarse graining steps in the time and spatial directions, respectively. These values are chosen such that

$$n\geq p.$$

(28)

As illustrated in Fig. 3, when $n>p$, the coarse graining procedure is as follows: 1.) Apply $p$ coarse-graining iterations in both directions using two-dimensional HOTRG. 2.) Continue with $(n-p-1)$ coarse-graining iterations in the spatial direction using the one-dimensional HOTRG. 3.) Perform an exact contraction to the last two remained tensors $A^{[r,n-1]}$ with the periodic boundary condition (PBC), obtaining

$$\left({\cal A}^{L_{\rm t}}\right)_{kj}\approx\sum_{l_{0},l_{1}}{A}^{[r,n-1]}_{k_{0}l_{0}j_{0}l_{1}}{A}^{[r,n-1]}_{k_{1}l_{1}j_{1}l_{0}}\eqqcolon{\cal A}^{[r,n]}_{k_{0}k_{1},j_{0}j_{1}}.$$

(29)

Meanwhile, the coarse graining procedure when $n=p$ is as follows: a.) Apply $(n-1)$ coarse graining iterations in both directions using the two-dimensional HOTRG. b.) Perform an exact contraction to the last four remained tensors with PBC, resulting in

$$\left({\cal A}^{L_{\rm t}}\right)_{kj}\approx\sum_{l_{0},l_{1},m_{0},m_{1},n_{0},n_{1}}{A}^{[r,n-1]}_{k_{0}l_{0}n_{0}l_{1}}{A}^{[r,n-1]}_{k_{1}l_{1}n_{1}l_{0}}{A}^{[r,n-1]}_{n_{0}m_{0}j_{0}m_{1}}{A}^{[r,n-1]}_{n_{1}m_{1}j_{1}m_{0}}\eqqcolon{\cal A}^{[r,n]}_{k_{0}k_{1},j_{0}j_{1}}.$$

(30)

Subsequently, we apply the EVD to ${\cal A}^{[r,n]}$ as

$${\cal A}^{[r,n]}=W^{[r,n]}\lambda^{[r,n]}W^{[r,n]\dagger},$$

(31)

where $W^{[r,n]}$ is a unitary matrix that approximates $W$ in Eq. (21) and $\lambda^{[r,n]}$ denotes the numerical eigenvalues. These are related to the eigenvalues of $\cal A$ in Eq. (21) by

$$\lambda_{a}\approx\left(\lambda_{a}^{[r,n]}\right)^{1/L_{\rm t}}.$$

(32)

Finally, the estimated energy spectrum is given by

$$\omega_{a}\approx\frac{1}{L_{\rm t}}\log\left(\frac{\lambda_{0}^{[r,n]}}{\lambda_{a}^{[r,n]}}\right)\eqqcolon\omega_{a}^{[{\rm hotrg}]}\,\,\,\,\mbox{ for }a=1,2,3,\ldots.$$

(33)

Next, we explain the coarse graining procedure for the impurity tensor network ${\cal A}^{m-1}{\cal A}^{\prime}{\cal A}^{m}$ in Eq. (22). Similar to the pure tensor case, we begin by preparing the initial impurity tensor network. The bare impurity tensor is denoted by $A^{{}^{\prime}[0,0]}=A^{\prime}$. To obtain the $s$-th coarse-grained impurity tensor $A^{{}^{\prime}[s,0]}$, we perform

$$A^{\prime[s,0]}_{c_{0}a_{1}b_{0}a_{0}}=\vbox{\hbox{\includegraphics[height=85.35826pt]{picture/init_contraction_impure_eqn.png}}},\hskip 28.45274pt(s=1,2,\ldots,r)$$

(34)

where $A^{{}^{\prime}[s-1,0]}$ is the previously coarse grained impurity tensor, and $U^{[s,0]}$ is the same isometry with the pure tensor case obtained from Eqs. (26-27). After $r$ coarse-graining iterations, we obtain $A^{\prime[r,0]}$, which is then inserted to the main impurity tensor network, as shown in Fig. 4. The new tensor $A^{\prime[r,i]}$ in this main impurity tensor network is obtained from the previous tensors $A^{\prime[r,i-1]}$ and $A^{[r,i-1]}$. As in the pure tensor case, when $n>p$, the coarse graining of the impurity tensor network is carried out according to the procedure 0.)-3.) shown in Fig. 4. In the last step, we perform an exact contraction for the remaining one pure tensor and one impurity tensor along with PBC, and obtain the estimate of ${\cal A}^{m-1}{\cal A}^{\prime}{\cal A}^{m}$, that is

$$({\cal A}^{m-1}{\cal A}^{\prime}{\cal A}^{m})_{kj}\approx\sum_{l_{0},l_{1}}{A}^{\prime[r,n-1]}_{k_{0}l_{0}j_{0}l_{1}}{A}^{[r,n-1]}_{k_{1}l_{1}j_{1}l_{0}}\eqqcolon{\cal A}^{\prime[r,n]}_{k_{0}k_{1},j_{0}j_{1}}.$$

(35)

On the other hand, when $n=p$, we follow the procedure 0.)-b.). In the final step, an exact contraction of one impurity tensor and three pure tensors is performed to obtain

$$({\cal A}^{m-1}{\cal A}^{\prime}{\cal A}^{m})_{kj}\approx\sum_{l_{0},l_{1},m_{0},m_{1},n_{0},n_{1}}{A}^{{}^{\prime}[r,n-1]}_{k_{0}l_{0}n_{0}l_{1}}{A}^{[r,n-1]}_{k_{1}l_{1}n_{1}l_{0}}{A}^{[r,n-1]}_{n_{0}m_{0}j_{0}m_{1}}{A}^{[r,n-1]}_{n_{1}m_{1}j_{1}m_{0}}\eqqcolon{\cal A}^{\prime[r,n]}_{k_{0}k_{1},j_{0}j_{1}}.$$

(36)

After the building blocks of matrix elements have been estimated, we compute the approximate matrix elements given in Eq. (22), using the following expression

$$\langle b|\hat{\mathcal{O}}_{q}|a\rangle\approx\left(\left(\lambda^{[r,n]}\right)^{-(m-1/2)/L_{\rm t}}W^{[r,n]\dagger}{\cal A}^{\prime[r,n]}W^{[r,n]}\left(\lambda^{[r,n]}\right)^{-(m+1/2)/L_{\rm t}}\right)_{ba}\eqqcolon\langle b|\hat{\mathcal{O}}_{q}|a\rangle^{[\rm hotrg]}.$$

(37)

In Az-zahra et al. (2024), a square tensor network was employed for all calculations. Using this setup, we were able to investigate the spectrum up to system size $L_{\rm s}=64$, where the highest excited state for which the quantum numbers are correctly identified corresponds to eigenstate number $a=20$, and the error in the energy spectrum is of order $O(10^{-3})$. In the present work, we slightly modify the coarse-graining strategy to improve the accuracy of the calculations, thereby enabling access to higher eigenstates and larger system sizes. Specifically, we compute the spectrum at fixed $L_{\rm s}$ while varying system size in time direction $L_{\rm t}=2,4,\ldots,L_{\rm s}$, where $L_{\rm t}=L_{\rm s}$ correspond to a square tensor network. The case $L_{\rm t}=1$ is excluded, as it was demonstrated in Ref. Az-zahra et al. (2024) that this choice leads to large numerical error.

We compute $\omega_{a}^{[\rm hotrg]}$ using Eq. (33) for $L_{\rm s}=64$ and $L_{\rm t}=2,4,8,16,64$, following the procedure introduced in Sec. II.2. As shown in Fig. 5(a), the tensor network with $L_{\rm t}=4,8$ yields smaller relative error, $\delta\omega_{a}=\left|\frac{\omega_{a}^{[\rm hotrg]}-\omega_{a}^{[\rm exact]}}{\omega_{a}^{[\rm exact]}}\right|$, compared to other choices $L_{\rm t}=2,16,64$. Here, the exact spectrum $\omega_{a}^{[\rm exact]}$ is computed using Kaufman’s formula Kaufman (1949), see the appendix of Ref. Az-zahra et al. (2024) for the concise summary. There are two main reasons why $L_{\rm t}=4,8$ produce relatively smaller errors compared to the other choices $L_{\rm t}=2,16,64$. First, the number of the coarse-graining steps required for $L_{\rm t}=4,8$ is smaller than for $L_{\rm t}=16,64$ resulting in smaller accumulated coarse-graining error. Second, the eigenvalues of the transfer matrix tend to be more degenerate for smaller $L_{\rm t}$, particularly when $L_{\rm s}$ is large, which leads to large errors after each coarse-graining step. This explains why $L_{\rm t}=4,8$ yield smaller errors compared to $L_{\rm t}=2$. Furthermore, when comparing the results for $L_{\rm t}=4$ and $L_{\rm t}=8$, we observe that the low-lying eigenstates are determined more accurately for $L_{\rm t}=8$, whereas the higher eigenstates exhibit smaller error for $L_{\rm t}=4$. Since each choice, $L_{\rm t}=4,8$ has its own advantages and disadvantages, we employ both in the computations presented in this paper. In Fig. 5(b), we present the energy spectrum $\omega_{a}^{[\rm hotrg]}$ of the system for $L_{\rm s}=64$ computed using $\chi=80$ and $L_{\rm t}=4$. In principle, our scheme allows the extraction of $\omega_{a}^{[\rm hotrg]}$ up to $a=\chi^{2}$. However, we only present the result for which the numerical errors are in order $O(10^{-3})$. For the case $L_{\rm s}=64$, $L_{\rm t}=4$ and $\chi=80$, these correspond to the eigenstates $a=0-57$.

Next, we examine how the relative error of the energy spectrum changes over different cut-off bond dimension. We compute the energy spectrum for $L_{\rm s}=64$ and $L_{\rm t}=4,8$, using $\chi=70,80,90$, and present the results in Fig. 6. From the figure, we clearly observe that the relative error decreases as $\chi$ increases for both $L_{\rm t}=4,8$, as expected.

One of the quantum numbers for (1+1)d Ising model is $q=\{+1,-1\}$ which arises from the internal $\mathbb{Z}_{2}$ symmetry of the spin field. To classify the eigenstates according to this quantum number, we compute matrix elements of a spin field operator $\hat{s}_{x}$ (we use $\hat{s}_{0}$ at $x=0$ in practice), namely $\langle\Omega|\hat{s}_{x}|a\rangle$ where $\langle\Omega|$ is the ground state whose quantum number is $q_{\Omega}=+1$. From the selection rule of the discrete symmetry given in Eq. (13), if $\langle\Omega|\hat{s}_{0}|a\rangle\neq 0$, then the eigenstate $|a\rangle$ has a quantum number $q_{a}=-1$. To estimate the matrix elements, we coarse grain the corresponding impurity tensor network for $L_{\rm s}=64$ and $L_{\rm t}=4,8$ (with the impurity placed at $m=2,4$, respectively), using $\chi=80$, and subtitute the result into Eq. (37).

We present the energy spectrum $\omega_{a}^{[\rm hotrg]}$ previously shown in Fig. 5(b) with the additional corresponding quantum number in the upper panel of Fig. 7. The quantum number classification is based on $|\langle\Omega|\hat{s}_{0}|a\rangle^{[\rm hotrg]}|$ shown in the middle panel. Meanwhile, the bottom panel displays the exact quantum number for comparison, computed following Ref. Kaufman (1949). From the figures, using the matrix elements obtained with $L_{\rm t}=4$, we can correctly identify the quantum number for $L_{\rm s}=64$ up to $a=57$, whereas those obtained with $L_{\rm t}=8$ are accurate only up to $a=42$. In general, these results indicate that, by using tensor network with $L_{\rm t}=4,8$, we are able to identify more higher-energy eigenstates correctly compared to the previous results reported in Refs. Az-zahra et al. (2024, 2025) using square lattice.

The momentum of the eigenstate $|a\rangle$ in $q=-1$ sector can be identified by computing the matrix elements of an operator with momentum $P$

$$\langle\Omega|\hat{\mathcal{O}}_{1}(P)|a\rangle,$$

(42)

where

$$\hat{\mathcal{O}}_{1}(P)=\frac{1}{L_{\rm s}}\sum_{x=0}^{L_{\rm s}-1}\hat{s}_{x}e^{-iPx}$$

(43)

and $P=2\pi n/L_{\rm s}$ is the discrete momentum with $n=0,1,2,\ldots,L_{\rm s}-1$. The selection rule for momentum states that for a fixed $P$, if the following condition

$$\langle\Omega|\hat{\mathcal{O}}_{1}(P)|a\rangle\approx\langle\Omega|\hat{\mathcal{O}}_{1}(P)|a\rangle^{[\rm hotrg]}\neq 0$$

(44)

is satisfied, then the eigenstate $|a\rangle$ carries the momentum $P$. In Fig. 8, we show the impurity tensor network diagram corresponding to the operator $\hat{\mathcal{O}}_{1}(P)$. These networks are then coarse-grained to obtain $\langle\Omega|\hat{\mathcal{O}}_{1}(P)|a\rangle^{[\rm hotrg]}$.