Bulk and Surface Critical Behaviors of quantum Heisenberg antiferromagnet on a two-dimensional coupled diagonal ladders

We study several physical quantities to investigate the bulk phases and related phase transitions. The two transitions are associated with the spontaneously breaking of the spin rotational symmetry. The staggered magnetization is used to describe the Néel order,

where the staggered phase factor $\phi_{i}=\pm 1$ according to the sublattice to which site $i$ belongs. The Binder cumulant $U_{2}$ [26, 27] is defined using $m_{s}^{z}$

which is dimensionless at the critical point. $U_{2}$ converging to 1 with increasing system size indicates the existence of magnetic order, while tending to zero with increasing system size implies that the system is in the magnetic disordered phase.

The mean spin stiffness $\rho_{s}$ over the $x$ and $y$ directions is also calculated. It is related to the fluctuations of the winding number[28, 29]

where

is the winding number along the $a=x,y$ direction. Here, $N^{+}_{a}$ and $N^{-}_{a}$ denotes the total number of operators transporting spin in the positive and negative $a$ direction, respectively.

$\rho_{s}$ is none zero if the state is magnetically ordered and goes to zero when the system is in the magnetically disordered phase. The size dependence of the spin stiffness exactly at a QCP is expected as follows [30]

with $z$ the dynamic exponent and $d=2$ the dimensions of the model. In the case that $z=1$, the $\rho_{s}L$ is expected to be dimensionless at critical point.

Figure 3 plots $\rho_{s}L$ and $U_{2}$ as functions of $J_{\perp}$ for different system sizes. Clearly, the model is in the antiferromagnetic ordered state when $J_{\perp}$ is between 0.17 and 3. Since $U_{2}$ and $\rho_{s}L$ are dimensionless at a critical point, the crossings of curves for different sizes roughly indicate two transition points.

We adopt the standard $(L,2L)$ crossing analysis for $U_{2}$ and $\rho_{s}L$ to estimate the critical point and critical properties.[31] For $Q=U_{2}$ or $\rho_{s}L$, we define the finite-size estimator of the critical points $J_{c}^{(Q)}(L)$ as the crossing point of $Q(J_{\perp})$ curves for $L$ and $2L$, which drifts toward the critical point $J_{c}$ in the thermodynamic limit in the following way

where $\nu$ is the correlation length exponent, $\omega>0$ is the effective irrelevant exponent, and $a$ is an unknown constant. At the crossing point $J_{c}^{(Q)}(L)$, we define the finite-size estimator of exponent $\nu$ as

where $s^{(Q)}(L)$ is the slope of the curve $Q(J_{\perp})$ for size $L$ at $J_{c}^{(Q)}(L)$. $\nu^{(Q)}(L)$ converges to the exponent $\nu$ in the following way

with $b$ an unknown constant.

For $U_{2}$ and $\rho_{s}L$, the analyses yield consistent estimates of $J_{c}$ and $\nu$ within error bars. The results with higher accuracy are selected as the final results. All the results are listed in Table 1. In particular, the final estimates of the two critical points are $J_{c1}=0.17425(3)$ and $J_{c2}=2.99046(5)$.

To further determine the universal properties of the two critical points, we calculate the static spin structure factor and the spin correlation at the longest distance in a finite system at the two estimated critical points $J_{c1}$ and $J_{c2}$. The two quantities are defined based on the spin correlation function

where $\mathbf{r}_{ij}$ is the vector from site $i$ to $j$. The static spin structure factor at wave vector $(\pi,\pi)$ is defined as follows

where $\epsilon_{ij}=\pm 1$, depending on whether $i$ and $j$ belong to the same sublattice. The spin correlation function of half lattice size $C(L/2,L/2)$ averages $C(\mathbf{r}_{ij})$ between two spins $i$ and $j$ at the longest distance $\mathbf{r}_{ij}=(L/2,L/2)$.

$S(\pi,\pi)$ and $C(L/2,L/2)$ are used to extract the scaling dimension $y_{h}$ of the staggered magnetic field $h$ and the anomalous dimension $\eta$. At QCP, $S(\pi,\pi)$ and $C(L/2,L/2)$ satisfy the following finite size scaling forms

and

respectively, in which $d=2$ is the spatial dimension, $z=1$ is the dynamical critical exponent, and $\omega$ the effective correction to scaling exponent. The two exponents $y_{h}$ and $\eta$ are not independent and are expected to obey the following scaling relation

The numerical results of $S(\pi,\pi)/L^{2}$ and $C(L/2,L/2)$ as functions of system size $L$ at two critical points are shown in Fig. 4. We fit the data of $S(\pi,\pi)/L^{2}$ and $C(L/2,L/2)$ according to Eqs. (13) and (14), respectively, and find the critical exponents $y_{h}$ and $\eta$, as listed in Table 1. The obtained $y_{h}$ and $\eta$ satisfy the scaling relations Eq. (15).

Comparing with the best known exponents of the 3D $O(3)$ universality class[32, 33], we conclude that both critical points belong to the 3D $O(3)$ universality class. This also shows that the topological order does not change the universality class of the bulk phase transition.

The spin-1/2 Heisenberg diagonal ladder is shown as the composite spin representation of a spin-1 chain[23] by representing the spin-1 operator $\boldsymbol{\sigma}_{i}$ of the chain as the sum of two spin-1/2 operators $\boldsymbol{\sigma}_{i}=\mathbf{S}_{i,0}+\mathbf{S}_{i,1}$ on two legs, as illustrated in Fig. 5. The low-energy spectrum of the ladder is identical to that of the spin-1 chain. When periodic boundary conditions are applied, all spins are bound to form valence bonds and the ground state is unique. The ground state of this ladder can be well described by the AKLT state[34], which is a short-ranged valence-bond (VB) state. A typical configuration is shown in Fig. 2. The Haldane gap is related to the energy needed to break a valence bond. More importantly, with open boundaries, the ground states have two spin-1/2 spins localized at the ends of the ladder. This is evident in the diagonal ladders as shown in Fig. 5. The ladder is in a Haldane phase with symmetry protected topological order.

The Haldane phase with symmetry protected topological order is characterized by a nonlocal string order, which is evident when the $z$ component of the spins on the same rung are summed. The total $S_{i}^{z}=S^{z}_{i,0}+S^{z}_{i,1}$ can take the values of 1,0, -1. When the sites with $S_{i}^{z}=0$ are removed, the remaining sites have a Néel order, which means a string order, as shown in Fig. 5(a). This order can be described by the following string order parameter: [35, 23]

The name $\mathcal{S}_{\rm odd}$ comes from the topology of the VB’s in the diagonal ladder [23], which is determined by the parity of the number of VB’s crossing an arbitrary vertical line. Note that the VB state shown in Fig. 5 is odd.

Other ladders in the Haldane phase may show a topologically distinct string order, which can be defined as

which is nonzero when the VB configuration is even,i.e., an even number of VB’s crossing an arbitrary vertical line[23]. In the case of diagonal ladder, as shown in Fig. 5(b), when the $z$ component of the spins along the plaquette diagonals are summed, there is no such a string order. Apparently, the parity of the VB ground state is intimately related to the type of string order.

The finite value of a string order parameter at the limit $|i-j|\to\infty$ characterizes a stable topological order in the thermodynamic limit. In the simulations of a system of size $L$ with periodic boundaries, we calculate $\mathcal{S}_{\rm odd}(L/2)$ ($\mathcal{S}_{\rm even}(L/2)$) by averaging $\mathcal{S}_{\rm odd}(i,j)$ ($\mathcal{S}_{\rm even}(i,j)$ ) at the maximum available distance $|i-j|=L/2$ along an individual ladder. As shown in Fig. 6(a) and (b)), in the case $J_{\perp}=0$, $S_{\rm odd}(L/2)$ is finite, $S_{\rm even}(L/2)$ vanishes when $L$ goes infinity, as expected for a diagonal ladder. We find $\mathcal{S}_{\rm odd}$ converges to 0.374325(7).

Now consider that the diagonal ladders are coupled to form a 2D lattice by interladder couplings $J_{\perp}>0$. Due to the Haldane gap, the properties of the ground state are robust against a weak higher-dimensional coupling between ladders. There is no phase transition when $J_{\perp}$ is turned on to finite values less than $J_{c1}$. This means that the model has a gapped Haldane phase that is adiabatically connected to the AKLT states of diagonal ladders. However, it was predicted theoretically that the string order of the coupled spin-1 Haldane chains is not stable and decays exponentially for arbitrarily weak interchain coupling [36]. This prediction has been verified numerically recently in the 2D CHC model[13] for sufficient large system sizes.

We obtain similar results in the SPT Haldane phase of the current model. We find that the string order parameter $\mathcal{S}_{\rm odd}(L/2)$ decays exponentially with $L$, but much slower than the decay of string order parameter in the CHC model at the same inter-ladder/inter-chain couplings. The numerical results are shown in Fig. 6(a). Fitting according to the following formula [36]

we find $\alpha=0.0007437(8)$ for $J_{\perp}=0.04$, and $\alpha=0.003089(1)$ for $J_{\perp}=0.08$.

We also calculated $\mathcal{S}_{\rm even}(L/2)$ inside the Haldane phase. As expected, the even string order parameter values are much smaller than the odd one. However, interestingly, we find that $\mathcal{S}_{\rm even}(L/2)$ decays algebraically with system size, as shown in Fig. 6(b). We have tried to fit the data using the following scaling form

and find $\beta=4.2(1)$ for a single ladder, $\beta=3.97(5)$ at $J_{\perp}=0.04$, and $\beta=3.74(3)$ at $J_{\perp}=0.08$.

However, the hallmark of the SPT phase is not the string order, but the presence of nontrivial surface states that are gapless or degenerate. Our model is spatial anisotropy, we here consider two different surfaces, the Y surface and the X surface, exposed by cutting the lattice, see Fig. 1. To study the surface states, we calculate the surface parallel correlation $C_{\parallel}(L/2)$ which averages $C(\mathbf{r}_{ij})$ between two surface spins $i$ and $j$ at the longest distance $L/2$.

The results for the Y surface at two couplings $J_{\perp}$ in the SPT phase are plotted in Fig. 7(a). We see that $C_{\parallel}(L/2)$ decays with system size $L$ in a power law,

with $p=1.30(4)$ at $J_{\perp}=0.04$ and $p=1.12(2)$ at $J_{\perp}=0.08$, meaning that the surface states are gapless. This is easy to understand from the AKLT state shown in Fig. 5(a). With open boundaries, each ladder carries a spin-1/2 excitation at each end. An effective spin-1/2 antiferromagnetic Heisenberg (AFH) chain is formed at the Y surface by coupling these spins, which is gapless according to the Lieb-Schultz-Mattis theorem[17].

Meanwhile, the results of $C_{\parallel}(L/2)$ for the X surface at the same two couplings in the SPT phase have completely different finite-size behavior, as shown in Fig. 7(b). The data can be fitted using straight lines on a linear-log scale, meaning the correlation decay exponentially. Fitting the curves with

we obtain $a=10.56(2)$ at $J_{\perp}=0.04$ and $a=11.49(3)$ at $J_{\perp}=0.08$. The surface states on the X surface are gapped, because the X surface is a gapped diagonal ladder, in its AKLT state, weakly coupled to the bulk.

Now we move to the RS phase with $J_{\perp}>J_{c2}$. The nature of this trivial disordered phase can be understood by examining the limit $J_{\perp}\to\infty$, at which the lattice reduces to disjoint bonds, see Fig. 2. The ground state is the direct product of these rung singlets.

We have calculated $C_{\parallel}(L/2)$ along the Y surface. The results are plotted in Fig. 7(b). We see the correlations at $J_{\perp}=3.2$ and $J_{\perp}=4.0$ decaying exponentially. Fitting data according to Eq. (21), we find $a=7.2(2)$ and $3.46(6)$, respectively, indicating the Y surface states are gapped. Apparently, the surface can be understood as sitting in a state adiabatically connected to a product state of dimers.

However, the X surface can be considered as a chain of dangling spins, weakly coupled to the bulk, in the RS phase. Thus, we expect the surface states are gapless, forming a spin-1/2 Heisenberg chain. This is proven by our numerical results. As shown in Fig. 7(a), the correlation $C_{\parallel}(L/2)$ along the X surface at $J_{\perp}=3.2$ and $4.0$ decay in a power law. Fitting according to Eq. (20), we obtain the power $p=0.827(2)$ and $0.886(2)$, respectively.

We first study the surface critical behaviors associated with the bulk critical point $J_{c2}$ separating the Néel ordered phase from the RS phase.

We start with checking the SCB on the Y surface, referred to as “Y-c2”. The numerical result of $\chi_{s1}$ as a function of size $L$ is graphed in Fig. 8(a), and the results of $C_{\parallel}(L/2)$ and $C_{\perp}(L/2)$ as functions of $L$ are plotted in Fig. 8(b).

We fit the data of $C_{\parallel}(L/2)$ and $C_{\perp}(L/2)$ according to Eqs. (24) and (25) and find statistically sound estimates of $\eta_{\parallel}=1.318(31)$ and $\eta_{\perp}=0.682(6)$.

The finite-size scaling form Eq. (26) supplemented with a constant $c$ as non-singular contribution, i.e.,

is used to fit the data of $\chi_{s1}$. The estimate of $y_{h1}$ is 0.852(34), with $\omega$ setting to 1.

These surface exponents are listed in Tab. 2. They obey the scaling relations in Eqs. (27) and (28), and agree well with the universal class of the ordinary transition associated with the 3D O(3) universality class found in various classical and quantum phase transitions (see Table 3). This behavior is expected since the surface state on Y surface in the RS phase is gapped, as shown in Sec.III.3.

We then check the SCBs at critical point $J_{c2}$ on the X surface, referred to as “X-c2”. This is the case that the surface is made up of dangling spins.

The numerical result of $\chi_{s1}$ as a function of size $L$ is graphed in Fig. 8(a), and the results of $C_{\parallel}(L/2)$ and $C_{\perp}(L/2)$ as functions of $L$ are shown in Fig. 8(b). Data fitting according to Eqs. (24), (25), and (26) finds statistically sound estimations $\eta_{\parallel}=-0.560(8),\eta_{\perp}=-0.259(4)$, and $y_{h1}=1.780(2)$, satisfying the scaling relations Eqs. (27) and (28).

The three exponents are also listed in Tab.2. They are consistent or very close to the nonordinary SCBs found in the quantum critical points of the 3D O(3) universality class [4, 9, 10, 13]. This result supports the scenario that nonordinary SCB can be induced by the gapless surface mode on the dangling spin-1/2 surface, as explored in Sec.III.3.

We then study the SCBs associated with the bulk critical point $J_{c1}$ at which the SPT Haldane phase transfers to the O(3) symmetry broken Néel phase.

We first study the Y surface with associated SCB denoted by “Y-c1”. The surface does not consist of dangling spins. However, we have shown in Sec.III.2 that the state of the Y surface in the gapped SPT Haldane phase is gapless due to the spin-1/2 excitations located at the ends of each diagonal ladder.

The numerical results of $C_{\parallel}(L/2)$, $C_{\perp}(L/2)$, and $\chi_{s1}(L)$ as functions of $L$ are plotted in Fig. 8. Apparently, the SCBs on the Y surface (Y-c1) are similar to those of X-c2 at $J_{c2}$.

We fit the data of $\chi_{s1}$, $C_{\parallel}(L/2)$, and $C_{\perp}(L/2)$ according to Eqs. (26), (24), and (25), respectively, and obtain statistically sound results $y_{h1}=1.756(3)$, $\eta_{\parallel}=-0.511(2)$ and $\eta_{\perp}=-0.237(2)$, as listed in Tab. 2. These exponents satisfy the scaling relations Eqs. (27) and (28), and are consistent with or very close to the nonordinary SCBs found on the X surface at $J_{c2}$, as well as those nonordinary SCBs found at other quantum critical points of the 3D O(3) universality class [4, 9, 10, 13].

Finally, we check the SCBs on the X surface, referred to as “X-c1”. The surface is gapped in the Haldane phase, as shown in Sec.III.2, as a result, the surface transition must belong to the ordinary class. Our numerical results verified that the SCBs are of the ordinary type.

The numerical results of $C_{\parallel}(L/2)$, $C_{\perp}(L/2)$, and $\chi_{s1}(L)$ as functions of $L$ are plotted in Fig. 8. The curves share similar slopes of the corresponding Y-c2 curves. Fitting these results according to Eqs. (24), (25), and (26), we obtain $y_{h1}=0.82(1)$, $\eta_{\parallel}=1.36(5)$ and $\eta_{\perp}=0.69(2)$, as listed in Tab. 2. Again, they satisfy the scaling relations Eqs. (27) and (28), and agree well with the exponents of the ordinary class in the 3D O(3) universality class.