text:
# Phase diagram of the 3D Axial-Next-Nearest-Neighbor Ising model

## I Introduction

Periodically modulated magnetic structures have attracted scientific interest for several decades both experimentally and theoretically. A non-trivial phase diagram obtained by experimental measurements in cerium antimonide (CeSb) shows a variety of different commensurately ordered magnetic structures with the underlying lattice CeSb1 ; CeSb2 . The three-dimensional (3D) $`S=\frac{5}{2}`$ axial next-nearest-neighbor Ising (ANNNI) model has been considered as a theoretical candidate for CeSb since it exhibits a rich structure when it is treated by mean-field approximation Boehm . The 3D $`S=\frac{1}{2}`$ ANNNI model is another example that shows a non-trivial spin modulated phase — the so-called devil’s stairs. This model has been analyzed theoretically by various approaches, including high-temperature series expansions Redner ; Otmaa , low-temperature series expansions FischerSelke , mean-field approximations MFA , Monte-Carlo simulations MCS , an effective-field approximation Surda , free-fermion methods, a phenomenological renormalization, and other methods reviewed in Ref. RepBak ; SelkeRev . The Monte Carlo simulations have also been applied to the $`S=\frac{1}{2}`$ ANNNI model with a finite number of spin layers Selke . Recently, Henkel and Pleimling considered an anisotropic scaling at the Lifshitz point using the Wolff cluster algorithm and critical exponents have been calculated Henkel .

The purpose of this paper is to clarify the phase structure of the 3D $`S=\frac{1}{2}`$ ANNNI model. Our interest is to study the spin modulated phases at intermediate temperatures, particularly, the stability of commensurate phases. For this purpose we apply a numerical variational method, the tensor product variational method (TPVA), to the model. In Section II we introduce the 3D ANNNI model and briefly discuss the variational background of the TPVA applied to the system. We present the numerical results in Section III where we construct the global phase diagram of the model and analyze the spin modulated phases. We summarize the obtained results in Section IV. In Appendix, a numerical self-consistent optimizing process is reviewed and efficiency of the modified TPVA is discussed.

## II Model and Non-uniform Product Variational State

We study the $`S=\frac{1}{2}`$ ANNNI model on a simple cubic lattice with the size $`L\times \mathrm{}\times \mathrm{}`$ along the $`x`$, $`y`$, and $`z`$ directions, respectively. The model is described by the lattice Hamiltonian

$``$ $`=`$ $`J_1{\displaystyle \underset{i,j,k}{}}\sigma _{i,j,k}\left(\sigma _{i+1,j,k}+\sigma _{i,j+1,k}+\sigma _{i,j,k+1}\right)`$ (1)
$`+J_2{\displaystyle \underset{i,j,k}{}}\sigma _{i,j,k}\sigma _{i+2,j,k},`$
where the subscripts $`i`$, $`j`$, and $`k`$ of the Ising spin $`\sigma =\pm 1`$ refer to the $`x`$, $`y`$, and $`z`$ coordinates, respectively. The ferromagnetic interaction $`J_1>0`$ acts between the nearest-neighbors and $`J_2>0`$ is the competing antiferromagnetic interaction between the next-nearest-neighbors imposed only in the $`x`$ direction.

Figure 1 shows the layer-to-layer transfer matrix $`𝒯`$ to the $`z`$ direction which connects two adjacent spin layers $`[\sigma ]`$ and $`[\overline{\sigma }]`$ (each of the size $`L\times \mathrm{}`$ in the $`x`$ and $`y`$ directions). The transfer matrix can be exactly expressed as the product of partially overlapped local Boltzmann weights (cf. Fig. 1)

$$𝒯[\sigma |\overline{\sigma }]=\underset{i=1}{\overset{L2}{}}\underset{j=\mathrm{}}{\overset{+\mathrm{}}{}}W_{i,j}^\mathrm{B}\{\sigma |\overline{\sigma }\}.$$
(2)
We simplify the notations using a group of 6 spins

$$\{\sigma \}(\sigma _{i,j}\sigma _{i^{},j}\sigma _{i^{}{}_{}{}^{},j}\sigma _{i,j^{}}\sigma _{i^{},j^{}}\sigma _{i^{}{}_{}{}^{},j^{}}),$$
(3)
with the index rule $`i^{}=i+1`$, $`i^{}{}_{}{}^{}=i+2`$, and $`j^{}=j+1`$. The local Boltzmann weight $`W_{i,j}^\mathrm{B}`$ of the Hamiltonian in Eq. (1) has the following form

$`W_{i,j}^\mathrm{B}\{\sigma |\overline{\sigma }\}=\mathrm{exp}\{{\displaystyle \frac{1}{k_\mathrm{B}T}}[{\displaystyle \frac{J_1}{6}}(\sigma _{i,j}\overline{\sigma }_{i,j}+\sigma _{i^{},j}\overline{\sigma }_{i^{},j}`$ (4)
$`+\sigma _{i^{}{}_{}{}^{},j}\overline{\sigma }_{i^{}{}_{}{}^{},j}+\sigma _{i,j^{}}\overline{\sigma }_{i,j^{}}+\sigma _{i^{},j^{}}\overline{\sigma }_{i^{},j^{}}+\sigma _{i^{}{}_{}{}^{},j^{}}\overline{\sigma }_{i^{}{}_{}{}^{},j^{}})`$
$`+{\displaystyle \frac{J_1}{8}}(\sigma _{i,j}\sigma _{i^{},j}+\sigma _{i^{},j}\sigma _{i^{}{}_{}{}^{},j}+\sigma _{i,j^{}}\sigma _{i^{},j^{}}+\sigma _{i^{},j^{}}\sigma _{i^{}{}_{}{}^{},j^{}}`$
$`+\overline{\sigma }_{i,j}\overline{\sigma }_{i^{},j}+\overline{\sigma }_{i^{},j}\overline{\sigma }_{i^{}{}_{}{}^{},j}+\overline{\sigma }_{i,j^{}}\overline{\sigma }_{i^{},j^{}}+\overline{\sigma }_{i^{},j^{}}\overline{\sigma }_{i^{}{}_{}{}^{},j^{}})`$
$`{\displaystyle \frac{J_2}{4}}(\sigma _{i,j}\sigma _{i^{}{}_{}{}^{},j}+\sigma _{i,j^{}}\sigma _{i^{}{}_{}{}^{},j}+\sigma _{i,j}\sigma _{i^{}{}_{}{}^{},j}+\sigma _{i,j^{}}\sigma _{i^{}{}_{}{}^{},j})]\}`$
with $`k_\mathrm{B}`$ being the Boltzmann constant and the temperature $`T`$. For reasons of simplicity and brevity, we consider $`J_1`$=1 and $`k_\mathrm{B}`$=1 throughout all calculations rem1 .

We consider a variational problem for the transfer matrix $`𝒯[\sigma |\overline{\sigma }]`$. For a given trial state $`|\mathrm{\Psi }`$, the variational partition function per layer is given by

$$\lambda _{\mathrm{var}}(\mathrm{\Psi })=\frac{\mathrm{\Psi }|𝒯|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}=\frac{\underset{[\sigma ],[\overline{\sigma }]}{}\mathrm{\Psi }[\sigma ]𝒯[\sigma |\overline{\sigma }]\mathrm{\Psi }[\overline{\sigma }]}{\underset{[\sigma ]}{}(\mathrm{\Psi }[\sigma ])^2}.$$
(5)
The TPVA is a numerical variational method that assumes a trial function written by the product of local weights $`V`$. For the ANNNI model, $`\mathrm{\Psi }`$ is written in the product form of mutually overlapped local weights (cf. Fig. 2)

$$\mathrm{\Psi }[\sigma ]=\underset{i=1}{\overset{L2}{}}\underset{j=\mathrm{}}{\overset{+\mathrm{}}{}}V_{i,j}\{\sigma \},$$
(6)
where we have used the simplified notation in Eq. (3).

In order to study non-uniform spin modulated phases, the local variational weights $`V_{i,j}\{\sigma \}`$ must be position dependent along the $`x`$ direction. Each $`V`$ thus contains $`2^6=64`$ adjustable parameters. Since we have written the trial function $`\mathrm{\Psi }`$ as well as the transfer matrix $`𝒯`$ in the product forms, both the numerator of Eq. (5)

$$\mathrm{\Psi }|𝒯|\mathrm{\Psi }=\underset{[\sigma ],[\sigma ^{}]}{}\underset{i=1}{\overset{L2}{}}\underset{j=\mathrm{}}{\overset{+\mathrm{}}{}}V_{i,j}\{\sigma \}W_{i,j}^\mathrm{B}\{\sigma |\overline{\sigma }\}V_{i,j}\{\overline{\sigma }\}$$
(7)
and its denominator

$$\mathrm{\Psi }|\mathrm{\Psi }=\underset{[\sigma ]}{}\underset{i=1}{\overset{L2}{}}\underset{j=\mathrm{}}{\overset{+\mathrm{}}{}}(V_{i,j}\{\sigma \})^2$$
(8)
have also the product forms. Thus these quantities can be accurately calculated by means of renormalization techniques, particularly, we used the density matrix renormalization group (DMRG) White ; Nishino .

## III Results

Figure 3 shows the global phase diagram of the ANNNI model obtained by the TPVA. It consists of

* a paramagnetic (disordered) phase,
* a uniformly ordered ferromagnetic phase,
* an antiphase with the periodic spin alignment ($`\mathrm{}\mathrm{}`$) for which we use the notation $`2`$ in the following, and
* a rich area of spin modulated phases.

Region of the spin modulated phases separates the antiphase from the paramagnetic phase. The paramagnetic, ferromagnetic, and the modulated phases meet at the Lifshitz point $`P_\mathrm{L}`$. In our calculations, it is located at $`J_2^\mathrm{L}/J_1`$=0.26 and $`k_\mathrm{B}T_\mathrm{L}/J_1`$=3.83 and is in agreement with the latest Monte Carlo calculations carried out by Pleimling and Henkel $`J_2^\mathrm{L}/J_1`$=0.270(4) and $`k_\mathrm{B}T_\mathrm{L}/J_1`$=3.7475(50) (from Ref. Henkel ).

The resulting phase diagram does not contradict to previous knowledge of the model. The phase boundary lines separating the ferromagnetic phase, the paramagnetic phase, the antiphase, and the spin modulated phases coincide with those obtained using the Monte Carlo calculations MCS . We found new features of the model in the region of the modulated phases, where the Monte Carlo simulations have not yielded a satisfactory answer. Our results are thought of as a supplement to the achievements computed by the mean-field approximations MFA ; Surda  at higher temperatures and by the low-temperature series expansions valid at low temperatures $`k_\mathrm{B}T/J_14`$ FischerSelke .

In the rest of this section, we focus on the region of the modulated phases that contains a multitude of various commensurate and incommensurate phases. For example, in Fig. 3 we plotted a few narrow areas of typical commensurate phases such as $`3,2=(\mathrm{}\mathrm{})`$, $`3=(\mathrm{}\mathrm{})`$, $`4=(\mathrm{}\mathrm{})`$, etc. all enclosed by the dotted lines. Note that the widths of these phases are substantially narrower compared to the mean-field approximation MFA  and the effective-field approximation Surda .

### III.1 Wavelength analysis

We first explain relation between the conventional notation and and modulation wave length $`\lambda `$. The antiphase $`2=`$ has periodicity of 4 lattice sites, thus $`\lambda =4`$. Another examples is the high-order commensurate phase $`3,(3,2)^2`$ which represents the periodic spin sequence $`()`$ and yields $`\lambda =26/5`$.

We calculate the spin modulations along five representative lines as depicted in Fig. 4 with their ending points listed in Table I. When we obtain the spontaneous magnetization, we compute the corresponding wavelengths by means of the Fourier transform.

Figure 5 shows the spin polarizations $`\sigma _i`$ at the two parameter points: $`C_1`$ on the line $`A_1B_1`$ and $`C_2`$ on $`A_2B_2`$. These two points are chosen near the phase boundaries. The spin polarization at $`C_1`$ exhibits the commensurate phase $`3^5,2(\mathrm{}\mathrm{})`$ with the wavelength $`\lambda =17/3`$ (the upper graph). The lower graph shows the commensurate phase at $`C_2`$ with $`\lambda 16.7`$ on the same region along the $`x`$ direction.

Now, we give a brief discussion on influence of boundary conditions imposed to the system on the resulting spin polarization. In Fig. 6 we plot $`\sigma _i`$ for three different types of the boundary conditions. We consider a lattice with the size $`401\times \mathrm{}\times \mathrm{}`$ and analyze the data at $`C_2`$. On the upper graph, the spin polarization is calculated for the fixed boundaries on the left end (the spins are aligned to the ’up’ direction) and the free boundaries on the right end. The Fourier transform applied to the whole region $`i=0,1,\mathrm{},400`$ yields $`\lambda =16.7\pm 0.53`$. On the intermediate graph, the parallel fixed boundary conditions on both sides are imposed (the spins are aligned ’up’ at the ends). It gives $`\lambda =16.7\pm 0.41`$. Finally, the lower graph shows the anti-parallel fixed boundaries (’up’ on the left end and ’down’ on the right end) with $`\lambda =16.7\pm 0.41`$. The choice of the boundary conditions does not affect the numerical results significantly. The larger lattice size is considered, the less influence of the boundaries is obtained, especially, off of the phase boundaries.

In Fig. 7 we plot the wavelength with respect to $`J_2/J_1`$ calculated on the three lines $`A_1B_1`$, $`A_2B_2`$, and $`A_3B_3`$. The dotted line is a guide for the eye to point out this structure. Near the boundary with the ferromagnetic phase, the wavelength rapidly increases. This is contradictory to the known results coming from the mean-field approximation. The mean-field approximation yields the first-order transitions between ferromagnetic phase and the individual commensurate phases on the boundary line (in details, see Ref. MFA ).

In the inset of Fig. 7 we plot details of the wavelength $`\lambda `$ in the vicinity of the commensurate phase $`3`$. We observed that the commensurate phases $`3`$, $`3^2,2`$, and $`3,2`$ ’lock-in’ at small regions of $`J_2/J_1`$ (on the line $`A_3B_3`$) and the so-called ”devil’s stairs” behavior is observed RepBak . On the contrary, the stairs-like structure is not visible at higher temperatures, as seen on the line $`A_1B_1`$ near the paramagnetic boundary.

Figure 8 shows $`\lambda `$ on the line from the point $`B_0`$ to the Lifshitz point $`P_\mathrm{L}`$. The wavelength diverges toward the calculated Lifshitz point at $`J_2/J_1=0.26`$. The stairs-like structure reveals if we zoom in on the phase diagram. For this reason, we selected an area depicted by the dotted rectangle therein. The inset shows the magnified area with the commensurate phases $`5`$ and $`5,4`$ where the wavelength locks-in.

### III.2 Low temperature behavior

To compare our results with analytical predictions, particularly, with the low-temperature series expansions (LTSE) FischerSelke , we have selected the line $`A_4B_4`$ which corresponds to the temperature $`k_\mathrm{B}T/J_1=1.5`$. The computed data are plotted in Fig. 9. The commensurate phase $`3`$ ($`\lambda =6`$) locks-in and forms a well-visible plateau. Note that the mean-field calculations MFA  do not exhibit any phases with $`\lambda >6`$ at $`k_\mathrm{B}T/J_12`$. It should be also noted that both the mean-field approximation (at higher temperatures) and the LTSE (at very low temperatures) do not result the spring of phases with $`\lambda >6`$ from the multi-phase point $`J_2/J_1=0.5`$.

Figure 10 illustrates the stairs-like structure of $`\lambda `$ and corresponds to the magnified area shown by the dotted rectangle in Fig. 9. The commensurate phases lock-in at rational values of $`\lambda `$ and are separated by the high-order commensurate phases and (possibly) the truly incommensurate ones. Several commensurate phases are denoted above the stairs-like curve in Fig. 10.

The LTSE yields a spring of an infinity of the commensurate phases, such as $`3`$ and $`3,2^n`$ for $`n=1,2,3,\mathrm{}`$, which separate the ferromagnetic phase and the antiphase. The transition from the region of the spin modulated phases ($`\lambda >4`$) to the antiphase ($`\lambda =4`$) does not contradict to LTSE. The existence of additional intermediate phases, $`3,2^n,3,2^{n+1}`$ $`n=1,2,3,\mathrm{}`$, at a higher temperatures was later reported by the LTSE. Our results contain all these commensurate phases. Moreover, we found unpredicted phases. For example, the transition between the phases $`3`$ and $`3,2`$ is not of the first order as reported by LTSE. We found that these two phases are separated by many commensurate phases, e. g., $`3^n,2`$ and $`3,(3,2)^{n+1}`$ for $`n=1,2,3,\mathrm{}`$ and the others of higher-orders as reported in Ref. Surda . We obtained such rich spin modulated structure also for $`\lambda >6`$, see Fig. 10.

In Fig. 11 we depict our numerical results obtained at $`k_\mathrm{B}T/J_1=1.5`$ and compare them with the results obtained by LTSE. The notation $`q_2`$ corresponds to the antiphase wave vector $`q=\pi /2`$. Note that while the LTSE gives the first-order transition among individual commensurate phases, our calculations yield subsequent stairs-like structures among them. Moreover, the LTSE calculations do not yield the commensurate phases with $`\lambda >6`$.

We, therefore, conjecture that the ’complete’ devil’s stairs structure exists at intermediate temperatures. The complete devil’s stairs structure suggests there are no first-order transition RepBak .

Here, we summarize those commensurate phases which were obtained by the numerical analysis of this model. Between two main commensurate phases $`p`$ and $`p+1`$, where $`p=2,3,4,\mathrm{}`$, new high-order commensurate phases are present, such as $`p^{n1},p+1`$, $`p,(p+1)^n`$, $`p,(p,p+1)^{n1}`$, and $`(p,p+1)^n,p+1`$, with $`n=2,3,4,\mathrm{}`$ Subsequently, the following higher-order commensurate phases were found $`(p^{n+1},p+1)^m,p^n,p+1`$ and $`(p,(p,p+1)^{n+1})^m,p,(p,p+1)^n`$, for $`m=2,3,4,\mathrm{}`$ etc.

## IV Summary

We applied the modified TPVA to the 3D ANNNI model and obtained the global phase diagram. The location of the Lifshitz point agrees with calculations performed by the high-temperature series expansions Redner ; Otmaa  and the recent Monte Carlo calculations Henkel . The modulated phase exhibits very complex structures. We found that (1) the commensurate phases are substantially narrower than those reported so far, (2) the wavelength of the spin modulated (commensurate) phases diverges at the boundary with the ferromagnetic phase, (3) the commensurate phases merge at low temperatures tending toward the multi-phase point $`J_2/J_1=0.5`$ and at low temperatures, the wavelengths with $`\lambda >6`$ are obtained, and (4) many (possibly infinity) phases have been found within the modulated phase, which have not yet been reported.

## Acknowledgments

A. G. thanks A. Šurda for an interesting discussion about the incommensurate phases in the ANNNI model. This work has been partially supported by the Grant-in-Aid for Scientific Research from Ministry of Education, Science, Sports and Culture (Grant No. 09640462 and No. 11640376) and by the Slovak Grant Agency, VEGA No. 2/7201/21 and 2/3118/23. A.G. is also supported by Japan Society for the Promotion of Science (P01192).

## Appendix

### Optimizing process for the local variational weights

We briefly describe the optimizing process of finding out the local variational weights in order to maximize Eq. (5). This optimization is based on the self-consistent equation in the TPVA to achieve the minimum of the free energy. Numerical details in the TPVA has been reported in Refs. TPVA .

In order to maximize the variational partition function in Eq. (5), by a proper tuning of the local variational weights $`V`$, we define two objects. One is the matrix object $``$ that represents a punctured classical system Delgado  defined on the 2-layer spin system which corresponds to the numerator $`\mathrm{\Psi }|𝒯|\mathrm{\Psi }`$ of the variational partition function. It is defined as

$`_{i,0}\{\sigma |\overline{\sigma }\}`$ $`=`$ $`W_{i,0}^\mathrm{B}\{\sigma |\overline{\sigma }\}{\displaystyle \underset{\stackrel{~}{[\sigma ]},\stackrel{~}{[\overline{\sigma }]}}{}}{\displaystyle \underset{ki}{}}{\displaystyle \underset{\mathrm{}0}{}}V_{k,\mathrm{}}\{\sigma \}`$ (9)
$`\times `$ $`W_{k,\mathrm{}}^\mathrm{B}\{\sigma |\overline{\sigma }\}V_{k,\mathrm{}}\{\overline{\sigma }\}.`$
Analogously, the vector object $`𝒜`$ is the punctured system defined on the 1-layer spin system,

$$𝒜_{i,0}\{\sigma \}=\underset{\stackrel{~}{[\sigma ]}}{}\underset{ki}{}\underset{\mathrm{}0}{}V_{k,\mathrm{}}\{\sigma \}V_{k,\mathrm{}}\{\sigma \}.$$
(10)
The configuration sums in Eqs. (9) and (10) are taken over all the spin variables $`\sigma `$ except for the 6 ones at the center of the system. In particular, except for

$$\{\sigma \}=(\sigma _{i,0}\sigma _{i+1,0}\sigma _{i+2,0}\sigma _{i,1}\sigma _{i+1,1}\sigma _{i+2,1})$$
(11)
and analogously for $`\{\overline{\sigma }`$} in Eq. (9). The notations $`_{ki}_\mathrm{}0`$ exclude $`V_{i,0}\{\sigma \}`$ and $`V_{i,0}\{\overline{\sigma }\}`$ from the product. Having defined these two objects, the variational partition function can be transformed into the expression,

$$\lambda _{\mathrm{var}}=\frac{\underset{\{\sigma \},\{\overline{\sigma }\}}{}V_{i,0}\{\sigma \}_{i,0}\{\sigma |\overline{\sigma }\}V_{i,0}\{\overline{\sigma }\}}{\underset{\{\sigma \}}{}V_{i,0}\{\sigma \}𝒜_{i,0}\{\sigma \}V_{i,0}\{\sigma \}}.$$
(12)
Now, consider a variation of $`\lambda _{\mathrm{var}}`$ with respect to variations of the local variational weights

$$\frac{\delta \lambda _{\mathrm{var}}}{\delta \mathrm{\Psi }}\underset{i,j}{}\frac{\delta \lambda _{\mathrm{var}}}{\delta V_{i,j}}.$$
(13)
Carrying out the extremal condition, $`\delta \lambda /\delta V_{i,j}=0`$, the self-consistent equation for the local variational weights $`V_{i,j}`$ is then obtained

$$V_{i,0}^{\mathrm{new}}\{\sigma \}=\underset{\{\overline{\sigma }\}}{}\frac{_{i,0}\{\sigma |\overline{\sigma }\}}{𝒜_{i,0}\{\sigma \}}V_{i,0}\{\overline{\sigma }\}.$$
(14)
The improvement of $`V`$ is performed as

$$V_{i,0}\{\sigma \}=V_{i,0}\{\sigma \}+\epsilon V_{i,0}^{\mathrm{new}}\{\sigma \}$$
(15)
through 64-parameter local parameter adjust. The self-consistent relation, Eq. (14), is a non-linear equation since $`_{i,0}`$ and $`𝒜_{i,0}`$ themselves depend on $`V`$. The convergence parameter $`\epsilon `$ controls the rate at which the improvement process of $`V`$ is performed.

Consequently, we compute the free energy per strip

$$_{\mathrm{new}}=k_\mathrm{B}T\mathrm{ln}\lambda _{\mathrm{var}}$$
(16)
and compare with the free energy $`_{\mathrm{old}}`$ calculated with the previous $`V`$.

### Efficiency of the algorithm

After the trial state $`|\mathrm{\Psi }`$ is optimized, we calculate the spontaneous magnetization at a site

$$\sigma _{i,j}=\frac{\mathrm{\Psi }|\sigma _{i,j}|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}.$$
(17)
Since the competing interactions exist only along the $`x`$ direction, the system is translation invariant with respect to the $`y`$ and $`z`$ directions. Therefore, the spontaneous magnetization $`\sigma _{i,j}`$ is independent on $`j`$ and we used $`\sigma _i`$ instead.

In order to estimate the numerical accuracy of the improved TPVA, we compare the calculation of the critical temperature $`T_\mathrm{c}`$ in the pure Ising model, i. e., when $`J_2=0`$, with other numerical methods. Table II summarizes the obtained $`T_\mathrm{c}`$. It is obvious that the mean-field approximation overestimates $`T_\mathrm{c}`$ and does not give reliable results near the phase boundaries. The improved TPVA with the 64 variational parameters results better $`T_\mathrm{c}`$ than the original TPVA with 16 parameters TPVA ; rem2 .

We set up the convergence parameter $`|\epsilon |=10^2`$. Assuming any $`|\epsilon |10^2`$ is sufficient for the most cases.

In Fig. 12 we illustrate an example which demonstrates the systematic decay of the free energy during the DMRG sweeping process until it finally converges. After the DMRG infinite system method (ISM) is finished, the first left-right sweep (1LR) proceeds followed by the first right-left sweep (1RL) and so on. Each step of the finite system method decreases the free energy until its minimum is reached.