text:
# Decoupling of the 𝑆=1/2 antiferromagnetic zig-zag ladder with anisotropy

## I Introduction

Antiferromagnetic ladder systems have attracted much interest recently.  On the theoretical side they interpolate between the well studied antiferromagnetic chain  and two dimensional antiferromagnets.  The evolution between 1D and 2D spin-1/2 antiferromagnetism is not necessarily smooth, however. In particular, the $`n`$-leg ladder shows a remarkable alternating property in the spectrum as the number of legs is even or odd.  The spectrum has a gap for an even number of legs while it is gapless for an odd number of legs. This is similar to the difference between integer (spin-gap) and half-odd-integer (gapless) spin chains . In the limit of strong-coupling between the two chains the two-leg ladder is essentially composed of weakly interacting singlets that form across the rungs. The lowest excitation is the promotion of a rung singlet to a triplet with an excitation energy of the order of the inter-chain coupling. This spin gap remains nonvanishing even for small inter-chain coupling due to the fact that a single antiferromagnetic chain is critical.  In the case of purely Ising coupled chains the gap appears for all values of $`n`$

The antiferromagnetic zigzag ladder has also attracted interest recently, particularly in the context of experimental systems with low-dimensional magnetic structures like that of $`Cs_2CuCl_4`$ . It is also interesting from a theoretical point of view because it is a frustrated system (see Fig. 1). Indeed, the zig-zag ladder is equivalent to a single antiferromagnetic chain with next-nearest-neighbor interactions. In this paper we consider the spin-1/2 antiferromagnetic zigzag ladder with anisotropy. The isotropic case has been studied before  as a function of the coupling parameter, $`j=J_2/J_1`$, which is the ratio of the next-nearest-neighbor interaction, $`J_2`$, to the nearest-neighbor interaction, $`J_1`$. As $`j`$ increases, the system goes from gapless (single chain) to a dimer phase and then to a spiral phase, where the structure factor has a maximum at a momentum $`\pi /2<q<\pi `$. The system has a spin gap in these last two phases, and it therefore only displays short-range order. In the limit that the intra-chain interaction is much larger than the inter-chain interaction ($`j\mathrm{}`$) the two chains decouple and a gapless single chain behavior is recovered. It has been argued that this only happens, strictly speaking, at $`j=\mathrm{}`$: the spin gap becomes exponentially small as $`j`$ grows, but it remains non-vanishing.  Recently, on the other hand, it has been proposed that incommensurate quasi-long-range spin correlations should be observed if easy-plane ($`XY`$) anisotropy is included in the zigzag ladder.  This is argued to be due to the presence of a “twist” term that results from the inter-chain interaction. It has been proposed that there is one gapless mode and one mode with a gap in the regime of strong $`XY`$ anisotropy in the inter-chain coupling. Another prediction of this work is the existence of spontaneous local spin currents. This, however, has been refuted in ref. 16. Also, other recent numerical work  has failed to confirm the gapless nature of the groundstate in the anisotropic $`XY`$ case at weak interchain coupling. Recent Density Matrix Renormalization Group (DMRG) results  suggest, however, that the zig-zag ladder does indeed show a gapless chiral phase as predicted in ref. 15.

In this paper we use exact diagonalizations, the modified Lanczos method  and the Davidson method  to address the possibility of a transition from a spin-gap regime to a gapless regime as a function of $`j`$ when anisotropy is present in small $`S=1/2`$ antiferromagnetic zig-zag ladders. We compute various probes to identify the different phases and study their behavior close to the transition points. We study in particular the spin stiffness, the dimer correlation function, the structure factor, and analyze in detail the spectrum in the various parameter regimes. Since the zig-zag ladder effectively has both nearest-neighbor and second neighbor interactions, a stiffness tensor is required to account for these two types of interactions. The eigenvalues of this $`2\times 2`$ matrix then become the natural spin rigidities that we use to clarify the behavior of the system in the various regimes. The stiffness of a system is a particularly good measure of the long-range nature of the groundstate. Introducing twisted boundary conditions leads to a response in the energy if the quantum states are extended (gapless case). On the other hand, the energy is insensitive to a change in the boundary conditions if the quantum states are localized (spin-gap case). Therefore, the stiffness with respect to such a twist is positive if the system is gapless and it is zero if the system has a gap. Also, the dimer correlation function naturally signals the dimer phase while the structure factor is a natural way to detect and study the spiral phase.

Our results are consistent with a gapless excitation spectrum in the case of $`XY`$ anisotropy at weak interchain coupling. We obtain qualitatively similar results for the isotropic case, however. This is most likely a finite size effect due to the exponentially small spin gap that persists in the isotropic zig-zag at weak coupling in the thermodynamic limit.  In fact, we show that the phase diagram that is obtained from the analysis of the spectrum for finite systems may be consistent with the field theory prediction  after performing extrapolations to the thermodynamic limit.

The paper is organized as follows: In section II we present the model and the quantities to be calculated. In section III we present our results and in section IV we summarize the work. Technical details concerning exact diagonalization are given in the Appendix. We also briefly review the extrapolation technique to the thermodynamic limit here.

## II Model and Probes

The anisotropic zigzag ladder is defined by the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_1^{XY}{\displaystyle \underset{i}{}}\left(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+\right)+J_1^z{\displaystyle \underset{i}{}}S_i^zS_{i+1}^z`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}J_2^{XY}{\displaystyle \underset{i}{}}\left(S_i^+S_{i+2}^{}+S_i^{}S_{i+2}^+\right)+J_2^z{\displaystyle \underset{i}{}}S_i^zS_{i+2}^z.`$ (2)
The spin operators refer to spin $`S=1/2`$ states, while the summation $`i=1,\mathrm{},N`$ runs along the “rib” of the zig-zag ladder. We shall parameterize the interactions by the coupling parameter $`j=J_2^{XY}/J_1^{XY}`$ and by the anisotropy parameter $`J_1^z/J_1^{XY}=\mathrm{\Delta }=J_2^z/J_2^{XY}`$. (The isotropic case reduces to $`j=J_2/J_1`$ and $`\mathrm{\Delta }=1`$.) We will set $`J_1^{XY}=1`$ henceforth. Consider first the nearest-neighbor Heisenberg chain with anisotropy, which corresponds to both the weak-coupling ($`J_1=0`$) and to the strong-coupling ($`J_2=0`$) limits of the zig-zag ladder (see Fig. 1). The spectrum is gapless for the case of $`XY`$ anisotropy, $`|\mathrm{\Delta }|1`$, as shown by the Bethe ansatz.  The excitation spectrum consists of spin-$`1/2`$ particles dubbed spinons. Since flipping one spin represents a spin-$`1`$ excitation, the spinons can only be created in pairs. Therefore the conventional spin $`1`$ magnons are deconfined into spin-$`1/2`$ spinons that propagate incoherently. In the regime where $`\mathrm{\Delta }1`$, the groundstate is ferromagnetic. When $`1\mathrm{\Delta }1`$ the leading spin configuration is the Néel state with the staggered magnetization lying within the $`XY`$ plane. At $`\mathrm{\Delta }=1`$ the groundstate is again in a Néel state, but with a staggered magnetization that can point in any direction. Last, the spectrum shows a gap in the Ising regime at $`\mathrm{\Delta }>1`$. The groundstate, on the other hand, displays strict long-range Néel order, with the staggered magnetization directed along the $`z`$ axis.

We shall begin our study of the antiferromagnetic zig-zag ladder by analyzing the classical limit of the isotropic Heisenberg case first: $`J_1^{XY}=J_{}=J_1^z`$ and $`J_2^{XY}=J_{}=J_2^z`$ as $`S\mathrm{}`$. A spiral state $`S_i^+=Se^{i\theta _i}`$ yields an energy per site of $`E(\theta )=S^2J_{}\mathrm{cos}(2\theta )+S^2J_{}\mathrm{cos}\theta `$. This magnetic energy is minimized at a pitch angle $`\theta _0`$ that satisfies $`\mathrm{cos}\theta _0=\frac{1}{4}J_{}/J_{}`$ for inter-chain exchange couplings that are below a critical value $`J_{}^c=4J_{}`$. A ferromagnetic state on each chain occurs, on the other hand, at strong coupling $`J_{}>J_{}^c`$, with a pitch angle of $`\theta _0=\pi `$. The spins are thus arranged antiparallel in between chains. To summarize, the system is in a spiral phase for $`J_{}<4J_{}`$, while it is in a collinear phase for $`J_{}>4J_{}`$. The same holds true when only $`XY`$-coupling exists. In the case of Ising coupling only, on the other hand, we have the effective model $`H=J_{}^z_iS_i^zS_{i+2}^z+J_{}^z_iS_i^zS_{i+1}^z`$. There are two possible groundstates. The first is the collinear state defined by $`S_i^z=\frac{1}{2}`$, for $`i`$ even and $`S_i^z=\frac{1}{2}`$ for $`i`$ odd (this has a degeneracy $`2`$) with an energy per site of $`E=S^2(J_{}^zJ_{}^z)`$. The other state is the antiferromagnetic one defined by $`\pm S_i^z=\frac{1}{2}(1)^{i/2}`$ for $`i`$ even and $`\pm S_i^z=\frac{1}{2}(1)^{(i\pm 1)/2}`$, for $`i`$ odd (this has a degeneracy $`2\times 2`$), with an energy per site of $`E=S^2J_{}^z`$. We have an Ising antiferromagnet for $`J_{}^z<2J_{}^z`$ and a collinear Ising ferromagnet for $`J_{}^z>2J_{}^z`$, with a first order transition separating the two phases.

Consider now the stiffness tensor in the classical limit. Imposing a spiral spin configuration on the zig-zag ladder with a pitch angle $`\theta `$, the energy per site of the classical $`J_1J_2`$ model is then given as above by

$$e=\frac{E}{N}=S^2J_1\mathrm{cos}(\theta +\theta _1)+S^2J_2\mathrm{cos}(2\theta +\theta _2),$$
(3)
where we have added small twists $`\theta _1`$ and $`\theta _2`$ to the nearest-neighbor and the next-nearest-neighbor terms, respectively. For $`\theta _1=0=\theta _2`$, we have that $`\mathrm{cos}\theta _0=1`$ for $`J_1>4J_2`$ and that $`\mathrm{cos}\theta _0=\frac{1}{4}(J_1/J_2)`$ for $`J_1<4J_2`$, as stated above. The spin currents are then given by

$`j_1`$ $`=`$ $`{\displaystyle \frac{e}{\theta _1}}|_0=S^2J_1\mathrm{sin}\theta _0`$ (4)
$`j_2`$ $`=`$ $`{\displaystyle \frac{e}{\theta _2}}|_0=S^2J_2\mathrm{sin}2\theta _0`$ (5)
and the rigidity components by

$`\rho _{11}`$ $`=`$ $`{\displaystyle \frac{j_1}{\theta _1}}|_0=S^2J_1\mathrm{cos}\theta _0`$ (6)
$`\rho _{12}`$ $`=`$ $`{\displaystyle \frac{j_1}{\theta _2}}|_0=0`$ (7)
$`\rho _{22}`$ $`=`$ $`{\displaystyle \frac{j_2}{\theta _2}}|_0=S^2J_2\mathrm{cos}2\theta _0.`$ (8)
Note that both the spin currents and the stiffnesses are independent of the anisotropy parameter $`\mathrm{\Delta }`$ in the classical limit: the isotropic and the $`XY`$ anisotropic cases give the same results. In the collinear phase at $`J_1>4J_2`$, the spin currents vanish ($`j_1=0=j_2`$) and $`\rho _{11}=S^2J_1`$, $`\rho _{12}=0`$ and $`\rho _{22}=S^2J_2`$. On the other hand, in the spiral phase at $`J_1<4J_2`$ the local spin currents are non-vanishing: $`j_1=\pm S^2J_1[1(J_1/4J_2)^2]^{1/2}`$ and $`j_2=\frac{1}{2}j_1`$. However, the total spin current $`j_s=j_1+2j_2`$ is null. The stiffness tensor of this spiral phase is given by $`\rho _{11}=\frac{1}{4}S^2(J_1^2/J_2)`$, $`\rho _{12}=0`$ and $`\rho _{22}=S^2[J_2\frac{1}{8}(J_1^2/J_2)]`$. The natural stiffness associated with the total spin $`j_s`$ is the response to an external twist that satisfies $`\theta _2=2\theta _1`$, and is given by $`\rho _s=\rho _{11}+4\rho _{12}+4\rho _{22}`$. It reduces to $`\rho _s=S^2(J_14J_2)`$ in the collinear phase and to $`\rho _s=S^2[4J_2\frac{1}{4}(J_1^2/J_2)]`$ in the spiral phase. These results are displayed in Fig. 2. Here it is shown that $`\rho _s`$ is always positive and that vanishes at the classical transition point between the collinear ferromagnet to the spiral phase. Recall that the nature of the groundstate changes across this transition. Also, we remark that in the spiral phase only the stiffness $`\rho _s`$ and the total current $`j_s`$ are “well behaved”. The other components show spontaneous spin currents, while $`\rho _{22}`$ is not always positive (we will return to this point later while discussing the quantum case).

In the general quantum case we calculate the stiffness in the standard way.  We consider the Hamiltonian (1) with periodic boundary conditions imposing uniform twists around the $`z`$ axis:

$`H(\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_1^{XY}{\displaystyle \underset{i}{}}\left(S_i^+S_{i+1}^{}e^{i\theta _1}+S_i^{}S_{i+1}^+e^{i\theta _1}\right)`$ (9)
$`+`$ $`{\displaystyle \frac{1}{2}}J_2^{XY}{\displaystyle \underset{i}{}}\left(S_i^+S_{i+2}^{}e^{i\theta _2}+S_i^{}S_{i+2}^+e^{i\theta _2}\right)`$ (10)
$`+`$ $`J_i^z{\displaystyle \underset{i}{}}S_i^zS_{i+1}^z+J_2^z{\displaystyle \underset{i}{}}S_i^zS_{i+2}^z.`$ (11)
Here $`\theta _1`$ and $`\theta _2`$ are two independent twists that act separately along the interchain and intrachain directions, respectively. Expanding the exponentials to second order we obtain the form

$`H(\theta _1,\theta _2)`$ $`=`$ $`H(0,0)+\theta _1{\displaystyle \underset{i}{}}J_i^1+\theta _2{\displaystyle \underset{i}{}}J_i^2`$ (12)
$``$ $`{\displaystyle \frac{1}{2}}\theta _1^2{\displaystyle \underset{i}{}}T_i^1{\displaystyle \frac{1}{2}}\theta _2^2{\displaystyle \underset{i}{}}T_i^2`$ (13)
where

$`J_i^1`$ $`=`$ $`{\displaystyle \frac{i}{2}}J_1^{XY}\left(S_i^+S_{i+1}^{}S_i^{}S_{i+1}^+\right)`$ (14)
$`J_i^2`$ $`=`$ $`{\displaystyle \frac{i}{2}}J_2^{XY}\left(S_i^+S_{i+2}^{}S_i^{}S_{i+2}^+\right)`$ (15)
are the spin currents along the interchain and intrachain directions, respectively and where

$`T_i^1`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_1^{XY}\left(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+\right)`$ (16)
$`T_i^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_2^{XY}\left(S_i^+S_{i+2}^{}+S_i^{}S_{i+2}^+\right)`$ (17)
are the kinetic energy operators. Using second-order perturbation theory, we then obtain that the groundstate stiffness tensor, $`\rho _{\alpha \beta }=(^2E/\theta _\alpha \theta _\beta )|_0`$, is given by

$`\rho _{\alpha \beta }=`$ $``$ $`{\displaystyle \underset{i}{}}<0|T_i^\alpha |0>\delta _{\alpha ,\beta }`$ (18)
$``$ $`2{\displaystyle \underset{\nu 0}{}}{\displaystyle \frac{1}{E_\nu E_0}}{\displaystyle \underset{i,j}{}}<0|j_i^\alpha |\nu ><\nu |j_j^\beta |0>,`$ (19)
where $`|0>`$ is the groundstate and $`|\nu >`$ are the excited states ($`\alpha ,\beta =1,2`$). The groundstate is assumed to be non-degenerate. Both the spin current operator $`J^\alpha =_iJ_i^\alpha `$ and the kinetic operator $`T^\alpha =_iT_i^\alpha `$ commute with the translation operator, $`T`$, and conserve total spin $`S_z`$ (where $`\alpha =1,2`$). Therefore, the states $`|\nu >`$ in the stiffness formula (9) are the excited states within the subspace of a given magnetization and momentum that contains the groundstate.

As in the classical case discussed earlier, the stiffness can also be calculated directly  taking numerical derivatives of the energy with respect to small twists, $`\theta _1`$ and $`\theta _2`$, that are imposed on the system. This procedure requires care that $`\theta _1`$ and $`\theta _2`$ are small enough so that there be no level crossings. (We prefer to use the correlation function method, even though the results using both methods agree very well.) The change in energy due to small twists $`\theta _1`$ and $`\theta _2`$ takes the form

$$\delta E=\frac{1}{2}\rho _{11}\theta _1^2+\rho _{12}\theta _1\theta _2+\frac{1}{2}\rho _{22}\theta _2^2$$
(20)
in the absence of spontaneous spin currents. It is then natural to consider the eigenvalues of the stiffness tensor

$$\rho _\pm =\frac{1}{2}\left(\rho _{11}+\rho _{22}\right)\pm \sqrt{\left(\frac{\rho _{11}\rho _{22}}{2}\right)^2+\rho _{12}^2}$$
(21)
and the determinant $`D=\rho _+\rho _{}=\rho _{11}\rho _{22}\rho _{12}^2`$. These eigenvalues will be computed using expression (9) for the stiffness tensor in the next section.

We shall also calculate the correlation function associated with the dimerization via linear response theory. Imposing a small explicit dimerization, we consider the Hamiltonian

$`H(\delta _1,\delta _2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}J_1^{XY}{\displaystyle \underset{i}{}}\left(1+(i)^i\delta _1\right)\left(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+\right)`$ (22)
$`+`$ $`{\displaystyle \frac{1}{2}}J_2^{XY}{\displaystyle \underset{i}{}}\left(1+(i)^i\delta _2\right)\left(S_i^+S_{i+2}^{}+S_i^{}S_{i+2}^+\right)`$ (23)
$`+`$ $`J_1^z{\displaystyle \underset{i}{}}\left(1+(i)^i\delta _1\right)S_i^zS_{i+1}^z`$ (24)
$`+`$ $`J_2^z{\displaystyle \underset{i}{}}\left(1+(i)^i\delta _2\right)S_i^zS_{i+2}^z.`$ (25)
Once again using second order perturbation theory we obtain that the susceptibility $`\chi _{\alpha \beta }=(^2E/\delta _\alpha \delta _\beta )|_0`$ is given by

$`\chi _{\alpha \beta }`$ $`=2\mathrm{Re}{\displaystyle \underset{\nu 0}{}}{\displaystyle \frac{1}{E_\nu E_0}}{\displaystyle \underset{i,j}{}}(1)^{i+j}`$ (27)
$`<0|(T_i^\alpha +M_i^\alpha )|\nu ><\nu |(T_j^\beta +M_j^\beta )|0>`$
where

$`M_i^1`$ $`=`$ $`J_1^zS_i^zS_{i+1}^z`$ (29)
$`M_i^2`$ $`=`$ $`J_2^zS_i^zS_{i+2}^z.`$ (30)
In the dimer correlation function, the factor $`(1)^{i+j}`$ implies that the states $`|\nu >`$ are contained in the same magnetization subspace, but in the $`k=k_0\pm \pi `$ momentum subspace, where $`k_0`$ is the groundstate momentum. Last, the structure factor is defined in the usual way:

$$S(q)=\underset{r}{}e^{iqr}C_r$$
(31)
where the correlation function $`C_r`$ is defined by

$$C_r=\frac{1}{NS(S+1)}\underset{i}{}<0|\stackrel{}{S}_i.\stackrel{}{S}_{i+r}|0>.$$
(32)
and is normalized such that the local correlation function ($`r=0`$) is unity.

## III Results

We now proceed to study the $`XXZ`$ zig-zag model (1) using exact diagonalization of finite systems with sizes up to $`N=28`$. The full energy spectrum is obtained for the smaller system sizes, $`N16`$, while only the groundstate and the first few excited states can be determined for the larger system sizes. The eigenvectors and eigenvalues of Hamiltonian eq. (1) are then substituted into Eqs. (9,13,15) to compute the various correlation functions. Our main aim is to study the transitions between the various phases.

Spectrum. Let us first survey the energy spectrum that is displayed by these small systems. We shall keep track of important quantum numbers associated with each energy level, such as the momentum along the rib of the zig-zag, the spin and the parity. We shall also identify points in parameter space where low-lying levels cross, and use this to identify phase transitions in the system. This procedure is known to yield accurate transition points when applied to even relatively small systems.

Let us begin by determining the quantum numbers of the groundstate as a function of the size $`N`$ for the $`S=1/2`$ zig-zag antiferromagnet. Periodic boundary conditions are imposed throughout. The groundstate is a spin singlet in general due to the antiferromagnetic interactions. For strong enough coupling between chains, $`j=J_2/J_1<1/2`$, it has either momentum $`\pi `$ for $`N=4n+2`$ or momentum $`0`$ for $`N=4n`$. For weak enough coupling between chains $`j>1/2`$, on the other hand, the momentum oscillates between $`0`$ and $`\pi `$ as a function of the coupling parameter $`j`$ and of the system size $`N`$ There are several points along $`j`$ in this regime where the corresponding energy levels for these two momentum values cross. The groundstate is degenerate at these points, and this is reflected by peaks in the dimer correlation function (see Fig. 8). Such level crossings grow in number as the system size grows, and this indicates that the two singlet states in question are in fact degenerate in the thermodynamic limit. By the Lieb-Schultz-Mattis theorem,  this is consistent with a spin gap in the excitation spectrum that survives the thermodynamic limit in the weak-coupling regime $`j>1/2`$.

Consider now the specific case of an $`XXZ`$ zig-zag chain (1) with $`N=16`$ sites under periodic boundary conditions (see Table I). In the isotropic case, $`\mathrm{\Delta }=1`$, the states are organized into spin multiplets due to the $`SU(2)`$ spin invariance. Again, the antiferromagnetic interactions imply that the groundstate is a spin singlet ($`S_z=0`$) in general. The system has three well defined regimes: (a) strong-coupling, (b) intermediate coupling and (c) weak-coupling. In the strong-coupling limit (a), $`j=J_2/J_10`$, the groundstate has momentum $`k=0`$. The first excited state forms a spin triplet in such case, with $`k=\pi `$, while the second excited state is another spin singlet with $`S_z=0`$ and $`k=\pi `$. At $`j=j_{c1}^N0.24`$ there is a level crossing where the first excited states and the second excited state interchange (here $`j_{c1}^N`$ is the value of $`j`$ at which the level crossing occurs for the system with size $`N`$). In the thermodynamic limit the two lowest singlet states ($`S_z=0`$, $`k=0`$ and $`S_z=0`$, $`k=\pi `$) become degenerate and there is a finite gap to the next excited state (the triplet state $`S_z=0,\pm 1`$, $`k=\pi `$). Although the system begins to dimerize at this stage (b), antiferromagnetic correlations along the rib of the zig-zag remain dominant up to the Majumdar-Ghosh point  at $`j=1/2`$. The groundstate is doubly degenerate for any system size at this point, where the two states are perfectly dimerized. Antiferromagnetic correlations within chains of the zig-zag become dominant beyond this point at weaker couplings $`j>1/2`$ (see Fig. 9). Another level crossing occurs as $`j`$ increases to about $`j_{c2}^N1.6`$ such that the first excited states become two triplets with $`S_z=0,\pm 1`$, $`k=\pm \pi /2`$. The groundstate displays quantum numbers $`S_z=0`$ and $`k=0`$ at this stage (c), and it is no longer degenerate. This remains so as $`j\mathrm{}`$. It should be mentioned, however, that White and Affleck predict that a non-zero spin gap persists in the thermodynamic limit at weak-coupling $`j1`$ between chains, and that it becomes exponentially small as $`j`$ grows  (we will return to this point later).

The spectrum of the anisotropic $`S=1/2`$ $`XXZ`$ zig-zag ladder has also been studied previously in the strong-coupling regime up to the Majumdar-Ghosh line ($`0<j<1/2`$).  A gapless regime occurs for $`XY`$ anisotropy $`\mathrm{\Delta }1`$ and strong coupling $`j<j_{c1}(\mathrm{\Delta })`$, an Ising antiferromagnet along the rib of the zig-zag that shows a spin gap in the excitation spectrum occurs for $`\mathrm{\Delta }>1`$ and $`j<j_{c1}(\mathrm{\Delta })`$, and a dimer phase regime that also has a spin gap exists at $`j>j_{c1}(\mathrm{\Delta })`$ and any $`\mathrm{\Delta }`$ (here $`j_{c1}(\mathrm{\Delta })`$ is the transition line obtained from extrapolation to the thermodynamic limit). The line $`j=j_{c1}(\mathrm{\Delta })`$ separates the gapless phase from the dimer phase for $`\mathrm{\Delta }1`$, while it separates the dimer phase from the (Ising) Néel phase for $`\mathrm{\Delta }>1`$. The line at $`\mathrm{\Delta }=1`$ and $`j<j_{c1}`$ separates the $`XY`$ gapless phase from the Ising phase.

Consider again the specific case of $`N=16`$ sites in this instance, with anisotropy parameters $`\mathrm{\Delta }=0`$ or $`\mathrm{\Delta }=0.5`$ (see Table I). The groundstate is a singlet with $`S_z=0`$ and $`k=0`$ as $`j`$ increases from the strong coupling limit at $`j=0`$ up to $`j=1/2`$. The first excited state has degeneracy $`2`$, with spin $`S_z=\pm 1`$ and momentum $`k=\pi `$ inside this regime. The next excited state is again a singlet with $`S_z=0`$ and $`k=\pi `$. As $`j`$ increases there exists a line of points, $`j=j_{c1}^N(\mathrm{\Delta })`$, such that the energy level of the excited state with $`S_z=\pm 1`$, $`k=\pi `$ crosses the energy level of the excited state $`S_z=0`$, $`k=\pi `$ For $`j>j_{c1}^N(\mathrm{\Delta })`$ the two lowest states are the two singlets with $`S_z=0`$, $`k=0`$ and $`S_z=0`$, $`k=\pi `$. Again, these two states become degenerate (with a gap to the first excited state) in the thermodynamic limit. The groundstate must be a spin singlet with $`S_z=0`$ due to the antiferromagnetic couplings. This excludes the possibility that any $`S_z=\pm 1`$ state be degenerate with it. As a result, the level crossing between the $`S_z=\pm 1`$ and singlet states can be used as an accurate criterion to determine the phase transition between the gapless and the spin-gap regimes.  If we increase $`j`$ up to $`j=1/2`$, then the system is exactly degenerate for all system sizes. This is a feature of the Majumdar-Ghosh point, which has a perfectly dimerized groundstate. The behavior of the system does not change much for intermediate coupling $`j>1/2`$ beyond this point, with the exception that the momentum of the two lowest states interchanges between $`k=0`$ and $`k=\pi `$ as $`j`$ grows. At a larger value of $`j=j_{c2}^N`$ between $`1.21.6`$, however, a level crossing occurs between the $`S_z=0`$, $`k=0\mathrm{or}\pi `$ state (first excited state) and a state with $`S_z=\pm 1`$, $`k=\pi /2`$ (second excited state). Notice that the momentum of the first excited state is now $`k=\pi /2`$. This is to be expected in this regime since the two chains are weakly coupled, and the periodicity doubles. For $`j>2`$, in particular, the first excited state is now four-fold degenerate ($`S_z=\pm 1`$, $`k=\pm \pi /2`$), and we might expect to fall back into a gapless regime since the first excited state is not a spin singlet. This level structure is in fact the same as for $`j\mathrm{}`$.

In Fig. 3 we present the crossing points between the states $`S_z=0`$, $`k=\pi `$ and $`S_z=\pm 1`$, $`k=\pm \pi /2`$ that signal this transition as function of the number of sites for several values of the anisotropy parameter $`\mathrm{\Delta }`$. We have attempted to extrapolate the crossing points, $`j_{c2}^N`$, to the thermodynamic limit by employing a standard (BST) algorithm due to Bulirsch and Stoer.  The method is explained in the Appendix. The results of this extrapolation procedure are shown in Fig. 4, where $`j_{c2}`$ is plotted as a function of $`\mathrm{\Delta }`$. The curve $`j_{c2}`$ separates a spin-gap (dimer) phase from a gapless phase at small interchain couplings. As expected the value of $`j_{c2}`$ grows near the isotropic point. (It should tend to $`j=\mathrm{}`$ at $`\mathrm{\Delta }=1`$ according to White and Affleck. )

Physical Probes. Correlation functions can also be used to determine the nature of different thermodynamic groundstates. The spin rigidity, in particular, can discriminate between phases that do and do not have spin gaps. The stiffness of the nearest-neighbor spin-$`1/2`$ Heisenberg chain has been calculated exactly via the Bethe ansatz.  In the thermodynamic limit this solution yields

$$\frac{\rho _{11}}{J_1^{XY}}=\frac{\pi }{4}\mathrm{sin}\frac{\pi }{n}\frac{1}{\frac{\pi }{n}\left(\pi \frac{\pi }{n}\right)}$$
(33)
for anisotropies

$$\mathrm{cos}\frac{\pi }{n}=\frac{J_1^z}{J_1^{XY}}=\mathrm{\Delta },$$
(34)
where $`n`$ is a positive integer. This yields $`\rho _{11}/J_1^{XY}=1/4`$ in the isotropic case and $`\rho _{11}/J_1^{XY}=1/\pi `$ in the $`XY`$-case. In Fig. 5 we plot the exact Bethe ansatz stiffness  for the single chain and compare it with the numerical diagonalization results for $`N=8,12,16,20`$ sites. The spectrum is gapless from the $`XY`$-limit up to the isotropic Heisenberg case. Beyond this critical point the spectrum acquires a gap due to the Ising anisotropy and the stiffness goes to zero. The stiffness remains positive beyond the transition point for the small system sizes that we studied, however. The anisotropy at which $`\rho _{11}=J_1^{XY}/4`$ extrapolates nicely to the Heisenberg point, $`\mathrm{\Delta }=1`$, in the thermodynamic limit, $`N\mathrm{}`$, however.

The zigzag ladder (or the chain with next-nearest neighbor interactions), on the other hand, is not solvable by the Bethe ansatz, and so other methods become necessary to compute the spin rigidity. With this purpose in mind, exact diagonalization calculations of the isotropic spin-1/2 Heisenberg chain with next-nearest neighbors interactions where carried out first by Bonca et al.  on small systems, where the diagonal component $`\rho _{11}`$ of the stiffness tensor was computed (see also ref. 31). We have completed this study by calculating the remaining components of the stiffness tensor (9), including cases with anisotropy. In Fig. 6 we show the results for the various stiffnesses for the zigzag ladder in the isotropic case as a function of $`j`$ for $`N=16`$, while the same set of results are shown for the $`XY`$-case ($`\mathrm{\Delta }=0`$) in Fig. 7. In Figs. 6a and 7a we show the stiffnesses $`\rho _{11}`$, $`\rho _{22}`$ and $`\rho _{12}`$. For small $`j`$ we are in the limit of strong interchain coupling (small next-nearest neighbor interaction) and we recover the previous results:  $`\rho _{11}>0`$ and $`\rho _{22}<0`$. This is also found in the classical case. In the opposite weak-coupling limit of very high $`j`$, $`\rho _{11}`$ remains negative, but $`\rho _{22}`$ becomes positive (as in the classical case). The latter is consistent with the extreme limit of two decoupled antiferromagnetic chains ($`J_1=0`$). Also, $`\rho _{12}`$ becomes non-zero in the intermediate region where quantum fluctuations are stronger and where the transition between weak-coupling and strong-coupling occurs. Note that the classical analysis reveals that the fact that $`\rho _{12}`$ becomes nonzero and the fact that $`\rho _{11}`$ becomes negative are purely quantum effects.

In the initial exact diagonalization study for the isotropic case, Bonca et al. chose to identify the point at which the stiffness component $`\rho _{11}`$ turns negative with a transition point to a quantum disordered phase with a spin gap.  Extrapolating to the thermodynamic limit they estimated the transition point to occur at $`j_c=0.270\pm 0.005`$. We believe that it is a better idea to look at the eigenvalues (11) of the full stiffness tensor to determine possible phase transitions.  In Figs. 6b and 7b we plot $`\rho _+`$ for different system sizes. This stiffness is always positive for all $`j`$, while $`\rho _{}`$ is always negative. For small $`j`$, $`\rho _+`$ is positive and of order unity. As $`j`$ grows $`\rho _+`$ becomes close to zero at a value of $`j`$ that is close to the value where $`\rho _{11}`$ crosses zero, $`j_{\rho 1}^N`$. As $`j`$ increases further $`\rho _+`$ again begins to grow appreciable near the point at which $`\rho _{22}`$ turns positive, $`j_{\rho 2}^N`$. A finite size analysis reveals that as the size $`N`$ grows the first “zero” of $`\rho _+`$ (near $`j=j_{\rho 1}^N`$) occurs at smaller values of $`j`$ and extrapolates to a value close to the dimer transition $`j_{c1}`$ determined by the level-crossing method. Similarly, we expect that the second transition to a nonzero stiffness $`\rho _+`$ (near $`j=j_{\rho 2}^N`$) signals the transition to a gapless regime that extends up to $`j=\mathrm{}`$, and that it might therefore also signal a decoupling transition. This happens, however, both for $`\mathrm{\Delta }=0`$ and for $`\mathrm{\Delta }=1`$. We believe that this is due to a finite-size effect in the latter case as discussed above. The fact that finite-size effects in the spin stiffness of a single chain become larger as the $`XY`$ anisotropy decreases (see Fig. 5) supports this claim.

The level crossing at weak-coupling, $`j_{c2}^N`$, does not, however, correlate well with $`j_{\rho 2}^N`$. The crossings defined by $`j_{\rho 2}^N`$ appear at smaller values of $`j`$ as compared with $`j_{c2}^N`$. Although both have the same trend, apparently finite size effects are stronger in the calculation of the stiffness than in the determination of the level crossings. We have therefore limited ourselves to the extrapolation of the level crossing points $`j_{c2}^N`$.

In Fig. 8 we show the behavior of the dimer correlation function $`\chi _{11}`$. For $`j>j_{c1}^N`$ the dimer correlation function increases signaling the spontaneous dimerization observed in the thermodynamic limit. The various peaks signal level crossings of the groundstate between near-degenerate states with momentum $`k=0`$ and $`k=\pi `$ as $`j`$ varies. The Majumdar-Ghosh point at $`j=1/2`$ is special since the groundstate is doubly degenerate for all system sizes in such case. Beyond $`j1`$ the dimer correlation function $`\chi _{11}`$ becomes small. (The susceptibilities $`\chi _{12}`$ and $`\chi _{22}`$ are always quite small).

Finally, remnants of the spiral phase that exists in the classical limit, $`S\mathrm{}`$, are clearly apparent in the momentum dependence of the structure factor. In Fig. 9 we show the structure factor as a function of momentum for several values of $`j`$ in the Heisenberg case at $`\mathrm{\Delta }=1`$. We see that the location of the maximum shifts from $`k_{max}=\pi `$ (for $`j0.5`$) to a value $`\pi /2<k_{max}<\pi `$ when $`j>0.5`$, thereby signaling the spiral phase. The results are similar for $`\mathrm{\Delta }=0`$.

## IV Conclusions

The $`S=1/2`$ antiferromagnetic zig-zag ladder is a difficult problem to solve due to the intrinsic frustration and to the criticality displayed by both the strong and the weak-coupling limits. The weak-coupling limit is particularly difficult in the absence of anisotropy, in which case antiferromagnetic correlations diverge exponentially.  This renders any numerical study of finite systems hard.

In this paper, we have performed an analysis of the exact properties of such finite systems, looking at various correlation functions and the structure of the spectrum both in the isotropic and the anisotropic cases. We have looked at the spectrum and have computed the spin stiffness of the zig-zag ladder, and have thereby found evidence for a gapless regime at weak coupling that survives the thermodynamic limit in the case of $`XY`$ anisotropy. However, the isotropic and the anisotropic cases look qualitatively similar. We believe that this is due to a strong finite-size effect in the former case. This claim is supported by the increase of finite-size effects in the stiffness of a single $`S=1/2`$ antiferromagnetic chain with decreasing $`XY`$ anisotropy, as shown in Fig. 5.

It was previously shown  that the dimer transition can be accurately determined in relatively small systems by studying the level crossing of the first and second excited states after extrapolating to the thermodynamic limit. We have used a similar criterion to detect a possible transition from the dimer phase to a gapless phase at weak-inter chain coupling. Using standard extrapolation techniques (see Appendix) we have constructed a phase diagram (see Fig. 4) in this regime that is in agreement with a recent proposal for a gapless spiral phase in the presence of $`XY`$ anisotropy.

Also, we expect the spin stiffness $`\rho _+`$ to be a good measure of the nature of the spectrum. In particular, it can be used as an order parameter to distinguish gapless from spin-gap phases. A positive stiffness indicates a gapless excitation spectrum and a null stiffness indicates a net spin-gap. The antiferromagnetic zig-zag for spin $`S=1/2`$ showed an appreciable positive stiffness in the limit of strong interchain coupling (similar to the single chain case), a very small yet positive stiffness in the intermediate (spin-gap) regime, and then again an appreciable positive stiffness in the limit of weak interchain coupling (similar to two decoupled chains). This was true for all values of the ($`XY`$) anisotropy, $`\mathrm{\Delta }1`$. It is known from previous work,  however, that a spin gap is always expected to be present in the isotropic case at weak coupling. The stiffness should therefore remain zero in the thermodynamic limit at $`\mathrm{\Delta }=1`$ in the weak-coupling regime. We believe that the discrepancy between this expectation and our results is a strong finite size effect in the isotropic case that is due to the exponentially diverging antiferromagnetic correlations.  Clearly, larger systems need to be studied.

Acknowledgments

We have benefited from helpful discussions and comments from Didier Poilblanc, Alexander Nersesyan, Fabian Essler and Alexander Gogolin. We are also indebted to Ed Rezayi for suggesting the spectral analysis.

##

Exact Diagonalization. The size of the Hilbert space under consideration can be considerably reduced using the symmetries of the problem. The Hamiltonian commutes with the total spin operator $`S_T^z`$, with the translation operator $`T`$, the spin flip operator $`R`$, and the space reflection operator $`L`$ $`(iN+1i,i=1,\mathrm{},N)`$. In the absence of anisotropies the Hamiltonian also commutes with the total spin operator $`(\stackrel{}{S}_T)^2`$ and the energy levels come in spin multiplets.

The operator $`R`$ commutes with $`T`$ but, although $`L`$ commutes with $`R`$, it does not commute with $`T`$, since $`LT=T^1L`$.

The action of the local operators $`S_i^\pm `$ and $`S_i^z`$ is simply given in the direct product basis $`|m_1^z>\mathrm{}|m_N^z>`$, which are eigenvectors of $`S_T^z`$. In general, these states are not eigenvectors of the additional symmetries. We consider first the translations, obtaining the classes of states which are closed under them. One starts with a state $`|a>`$ and applies $`N_a`$ times the translation operator $`T`$ until one finds $`T^{N_a}|a>=|a>`$, where $`N_a`$ is necessarily a divisor of $`N`$. The state $`|a>`$ is the representative of this class. In combinatorial theory this is called a necklace and its periodic part, of length $`N_a`$, a Lyndon word. The other classes are formed proceeding in the same manner starting with other states, not already used, until all the states have been exhausted. In each class, since $`T^{N_a}=\widehat{1}`$, the possible eigenvectors of $`T`$ are $`t_{N_a}(k)=e^{\frac{i2\pi k}{N_a}}`$, with $`k=0,\mathrm{},N_a1`$, corresponding to the momentum $`p_k=\frac{2\pi k}{N_a}`$. The operators $`P_N(k)=\frac{1}{N}_{i=0}^{i=N1}(\frac{T}{t_N(k)})^i`$ are projectors, i.e, they satisfy $`P_N(k)P_N(k^{})=\delta _{k,k^{}}P_N(k)`$. They also have the property $`TP_N(k)=t_N(k)P_N(k)`$. The projector $`P_N(k)`$ acting on the states of a class with $`N_a`$ elements reduces to $`\frac{N}{N_a}P_{N_a}(k)`$, if $`t_N(k)`$ is one of the $`t_{N_a}(k)`$ and gives zero otherwise. The normalized eigenvector of momentum $`p_k`$ formed from the class of the state $`|a>`$ is given by $`|k,a>=\sqrt{N_a}P_{N_a}(k)|a>`$. Since the Hamiltonian commutes with $`S_T^z`$ and $`T`$, the subspaces of fixed magnetization and momentum are invariant subspaces and it is important to consider each of these subspaces separately. A general state is written as a linear combination of the states $`|k,a>`$ corresponding to that magnetization. The state $`|k,a>`$ is represented by $`|a>`$. This allows us to reduce the size of the basis to the number of representatives.

The spin reflection changes the sign of the magnetization. Therefore it is relevant only for the classification of the states when the magnetization is zero. IF $`R|a>`$ is not in the class of $`|a>`$ both eigenvectors with $`r=1`$ and $`r=1`$ can be formed. IF $`R|a>`$ is in the class of $`|a>`$, it must be obtained from this state by $`\frac{N_a}{2}`$ translations (implying that $`N_a`$ must be even), since for half-integer spin it is linearly independent of $`|a>`$. The state $`|k,a>`$ is an eigenvalue of $`R`$ with $`r=1`$ for $`k`$ even and $`r=1`$ for $`k`$ odd.

The space reflection operator reverses the momentum, since one has $`LP_N(k)=P_N(k)L`$. The state $`L|a>`$ can be in a different class or in the same class. The states $`|k,a>`$ and $`|k,a>`$ are linearly independent, except for $`k=0`$, with eigenvalue $`l=1`$, or $`k=\frac{N_a}{2}`$, for $`N_a`$ even, with eigenvalue $`l=1`$. As a result, one sees that, for zero magnetization and $`t=\pm 1`$, it is possible to use all the operators $`T,R`$ and $`L`$ in the construction of the states. This allows a further reduction of the size of the subspace under consideration. Since some classes have a definite eigenvalue of $`R`$ or $`L`$ they are simply excluded if their values are not those of the state which we are studying.

Since the complete Hamiltonian commutes with the operator $`T`$ only transitions to states with the same momentum are allowed, even if separate terms in the Hamiltonian allow them. If $`H|a>=_b\alpha T^{p(a,b)}|b>`$, where $`p(a,b)`$ is an integer and $`|b>`$ is the representative of a class, one finally finds $`<k,b|H|k,a>=_b^{}\alpha \sqrt{\frac{N_a}{N_b}}\left(t_{N_a}(k)\right)^{p(a,b)}`$, where the summation is restricted to the allowed transitions.

The modified Lanczos method  is very useful to obtain the ground state of an Hamiltonian. Restricting our attention to a subspace of fixed magnetization and momentum we will find the ground state of the block Hamiltonian in that subspace. To obtain an approximate ground state we choose a trial state $`|\psi _0>`$ that can not be orthogonal to the true ground-state. We define a state $`|\psi _1>`$ as

$$|\psi _1>=\frac{\widehat{H}|\psi _0><H>|\psi _0>}{\left(<H^2><H>^2\right)^{1/2}}$$
(35)
where $`<\psi _1|\psi _1>=1`$, $`<\psi _1|\psi _0>=0`$ and $`<H^n>=<\psi _0|\widehat{H}^n|\psi _0>`$. Defining a matrix of the Hamiltonian in the basis $`|\psi _0>,|\psi _1>`$ we can diagonalize it  obtaining a next order approximation for the energy $`\stackrel{~}{ϵ}_0`$, and ground state $`|\stackrel{~}{\psi }_0>`$, with

$`\stackrel{~}{ϵ}_0`$ $`=`$ $`<H>+b\alpha `$ (36)
$`|\stackrel{~}{\psi }_0>`$ $`=`$ $`{\displaystyle \frac{|\psi _0>+\alpha |\psi _1>}{\left(1+\alpha ^2\right)^{1/2}}}`$ (37)
where

$`b`$ $`=`$ $`\left(<H^2><H>^2\right)^{1/2}`$ (38)
$`\alpha `$ $`=`$ $`f\left(f^2+1\right)^{1/2}`$ (39)
and

$$f=\frac{<H^3>3<H><H^2>+2<H>^3}{2\left(<H^2><H>^2\right)^{3/2}}.$$
(40)
Taking $`|\stackrel{~}{\psi }_0>`$ as the new $`|\psi _0>`$ we can iterate the method to obtain a better estimate for the energy and the ground state.

Since the ground state and the first excited of the system have different momenta they are both ground states of different subspaces, and we apply the modified Lanczos method to these subspaces. If the system does not have anisotropies one can alternatively look for the component of the spin multiplet of the first excited state with magnetization $`1`$.

Extrapolation. The results for the infinite system can be estimated using standard extrapolation methods like the BST method.  In the BST algorithm we look for the limit of a sequence of the type $`T(h_j)=T+a_1h_j^\omega +a_2h_j^{2\omega }+\mathrm{}`$, $`j=0,\mathrm{},N_p1`$ ($`N_p`$ being the number of data points), where $`h_j`$ is a sequence converging to zero as $`j\mathrm{}`$, corresponding to different system sizes $`N_j`$. Typically $`h_j=\frac{1}{N_j}`$, with $`N_j=a+bj`$, for some values of $`a,b`$. The value of the $`m^{th}`$ iteration for the sequence is obtained from

$$T_m^{(j_m)}=T_{m1}^{(j_m+1)}+\frac{T_{m1}^{(j_m+1)}T_{m1}^{(j_m)}}{\left(\frac{h_{j_m}}{h_{j_m+m}}\right)^\omega \left(1\frac{T_{m1}^{(j_m+1)}T_{m1}^{(j_m)}}{T_{m1}^{(j_m+1)}T_{m2}^{(j_m+1)}}\right)1}$$
(41)
with $`j_m=0,\mathrm{},M(m)`$ and where $`M(m)`$ is the number of values of $`T_m^{(j_m)}`$ at each iteration. It decreases by one, at each iteration, from $`M(0)=N_p`$ to $`M(N_p1)=1`$, when the iteration process is fulfilled. As initial values one defines $`T_1^{(j)}=0`$ and $`T_0^{(j)}=T(h_j)`$. The extrapolated value is $`T_{N_p1}^{(0)}`$ and the estimated error is $`ϵ=|T_{N_p2}^{(1)}T_{N_p2}^{(0)}|`$. Finally, $`\omega `$ is a free parameter which is adjusted such that the estimate of the error is a minimum.

Table I- Lowest energy levels for $`N=16`$. The states are represented by their $`S_z`$ values and momenta $`(S_z;k)`$. When more than one state is represented this means they are degenerate. When the momenta are not $`0`$, $`\pi `$ or $`\pi /2`$ the momentum of the state is represented by an integer, $`n`$, such that $`k=\frac{2\pi }{N}n`$.

$`\begin{array}{cccc}j\hfill & \mathrm{\Delta }=0\hfill & \mathrm{\Delta }=0.5\hfill & \mathrm{\Delta }=1\hfill \\ & & & \\ & (2;\pi )\hfill & (0;0)\hfill & (0,1,2;0,\pi )\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0;\pi /2)\hfill \\ & (1;\pi /2)\hfill & (1;\pi /2)\hfill & (0,1;\pi /2)\hfill \\ 10\hfill & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \\ & & & \\ & (2;\pi )\hfill & (0;0)\hfill & (0;\pi /2)\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0;\pi )\hfill \\ & (1;\pi /2)\hfill & (1;\pi /2)\hfill & (0,1;\pi /2)\hfill \\ 2\hfill & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \\ & & & \\ & (2;\pi )\hfill & (0;\pi /2)\hfill & (0;\pi /2)\hfill \\ & (0;\pi )\hfill & (1;\pi /2)\hfill & (0,1;\pi /2)\hfill \\ & (1;\pi /2)\hfill & (0;\pi )\hfill & (0;\pi )\hfill \\ 1.5\hfill & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \\ & & & \\ & (0;\pi )\hfill & (0;\pi )\hfill & \\ & (1;\pi /2)\hfill & (1;\pi /2)\hfill & (0,1;\pi /2)\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0;\pi )\hfill \\ 1\hfill & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \\ & & & \\ & (1;3)\hfill & (1;3)\hfill & (0,1;\pi )\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0,1;3)\hfill \\ & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \\ 0.75\hfill & (0;\pi )\hfill & (0;\pi )\hfill & (0;\pi )\hfill \\ & & & \\ & (0;0)\hfill & (1;1)\hfill & (0,1;1)\hfill \\ & (1;\pi )\hfill & (1;\pi )\hfill & (0;0)\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0,1;\pi )\hfill \\ 0.5\hfill & (0;0,\pi )\hfill & (0;0,\pi )\hfill & (0;0,\pi )\hfill \\ & & & \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0,1;1)\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0,1;\pi )\hfill \\ & (1;\pi )\hfill & (1;\pi )\hfill & (0;\pi )\hfill \\ 0.25\hfill & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \\ & & & \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0,1;1)\hfill \\ & (0;\pi )\hfill & (0;\pi )\hfill & (0;\pi )\hfill \\ & (1;\pi )\hfill & (1;\pi )\hfill & (0,1;\pi )\hfill \\ 0.1\hfill & (0;0)\hfill & (0;0)\hfill & (0;0)\hfill \end{array}`$

Figure Captions

Fig. 1- Diagram of the zig-zag ladder. The nearest-neighbor interaction is parameterized by $`J_1`$ and the next-nearest-neighbor interaction is parameterized by $`J_2`$.

Fig. 2- Classical stiffnesses $`\rho _{11}/(J_1S^2)`$, $`\rho _{22}/(J_1S^2)`$ and $`\rho _s/(J_1S^2)`$ as a function of $`j=J_2/J_1`$.

Fig. 3- Crossing points of the first excited state as a function of $`N`$ for different values of $`\mathrm{\Delta }`$ at weak coupling.

Fig. 4- Phase diagram in the regime of weak interchain coupling. The critical coupling parameter $`j_{c2}`$ is obtained from the crossing of the first excited state with the second excited state after extrapolating to the thermodynamic limit (see the Appendix).

Fig. 5- Comparison of the exact Bethe ansatz result for the spin stiffness of a single chain to the numerical results obtained from diagonalizing small systems as a function of $`\mathrm{\Delta }`$.

Fig. 6- (a) Stiffnesses $`\rho _{11}`$, $`\rho _{12}`$ and $`\rho _{22}`$ as a function of $`j`$ for $`\mathrm{\Delta }=1`$ and $`N=16`$. (b) $`\rho _+`$ as a function of $`j`$ for various system sizes for $`\mathrm{\Delta }=1`$.

Fig. 7- (a) Stiffnesses $`\rho _{11}`$, $`\rho _{12}`$ and $`\rho _{22}`$ as a function of $`j`$ for $`\mathrm{\Delta }=0`$ and $`N=16`$. (b) $`\rho _+`$ as a function of $`j`$ for various system sizes for $`\mathrm{\Delta }=0`$.

Fig. 8- Dimer correlation function for $`N=16`$ and $`\mathrm{\Delta }=1`$ and $`\mathrm{\Delta }=0`$.

Fig. 9- Structure factor, $`S(q)`$, for $`N=16`$ and $`\mathrm{\Delta }=1`$ for various coupling strengths. The peak is located at $`q=\pi `$ in the single-chain limit, $`j=0`$ but this momentum tends to $`q=\pi /2`$ as $`j`$ increases. This indicates the presence of important spiral correlations in the weakly-coupled $`S=1/2`$ zig-zag antiferromagnet.


format: text

text:
# Multiscale Reference Function Analysis of the 𝒫𝒯 Symmetry Breaking Solutions for the 𝑃²+𝑖𝑋³+𝑖𝛼𝑋 Hamiltonian

## I Introduction

There has been much interest, recently, in understanding the symmetry breaking mechanism for $`𝒫𝒯`$ invariant Hamiltonians of the type $`P^2+_{j=0}^JC_J(iX)^j`$. The mathematical interest in these systems originated from a conjecture by D. Bessis, and WKB analysis confirmation by Bender and Boettcher (1998), that the class of potentials of the form $`V(x)=(iX)^n`$ only allow for $`𝒫𝒯`$-invariant solutions, and thus can only have real discrete spectra. The recent literature testifies to the great interest in these problems, as can be found in the cited references by Bender et al (2001), and Mezincescu (2000, 2001).

An important work establishing that $`𝒫𝒯`$-symmetry breaking systems do exist was the recent study by Delabaere and Trinh (2000) which used asymptotic methods to analyze the Hamiltonian $`H_\alpha P^2+iX^3+i\alpha X`$. Their analysis showed the existence of symmetry breaking solutions for large $`\alpha `$ values; however, for moderate values, their results were near the limits of their analytical validity. Despite this, they made certain questionable predictions for moderate $`\alpha `$ values.

Our immediate objective is to check the validity of the Delabaere and Trinh study, by explicitly computing the low-lying complex (and real) eigenenergies for moderate $`\alpha `$ values. We do this in two ways. The first makes use of the very efficient Multiscale Reference Function (MRF) formalism of Tymczak et al (1998a,b). The results of this eigenenergy estimation analysis are then confirmed through application of a recently developed eigenenergy bounding theory proposed by Handy (2001a,b), which can generate converging bounds to the complex eigenenergies. This bounding approach, referred to as the Eigenvalue Moment Method (EMM), is exact, although numerically slower in its implementation than the MRF procedure.

As a footnote to the above, we emphasize that the EMM theory generates an infinite hierarchy of closed, finite dimensional, algebraic, eigenenergy constraints. These are then solved, numerically. In this regard, the EMM procedure is very different from other numerical schemes, such as numerical integration, which are intrinsically of an approximating nature, and cannot provide any fundamental theoretical insight into the underlying physical processes. The algebraic constraints generated through the EMM formalism can provide such insight; although such analysis have not been attempted, heretofore. However, for the immediate purposes of this study, we solely defer to EMM in order to check the validity of the MRF results.

In this work, we provide the essentials of the MRF theoretical structure, as applied to the $`H_\alpha `$ Hamiltonian. The EMM theory is not discussed. Only the numerical bounds are quoted in the Tables. This work validates the relevancy of MRF theory in the computation of complex eigenenergies, for the class of problems referenced above.

Both the MRF and EMM methods are dependent on a moments’ representation for the given system. This in turn is readily realizable for any (multidimensional) rational fraction potential.

Any moment based analysis is inherently multiscale in nature. That is, as the number of moments used increases, one is probing the system at successively smaller scales. Consistent with this, the MRF basis representation, particularly within configuration space, has important ties with (complex) turning point quantization (Handy et al (2000)), and wavelet analysis (Handy and Brooks (2001)).

We outline the basic MRF theory and its implementation, in the next section. The last section contains a detailed enumeration, and illustration, of the MRF results, which make precise the qualitative spectral structure conjectured by Delabaere and Trinh.

## II  The MRF Representation

### A The Moment Equation

The starting point for the MRF analysis is the transformation of the Schrodinger equation into the Fourier representation, assuming that one is working with the physical, $`L^2`$, solutions. Thus, for the configuration space Schrodinger equation studied by Delabaere and Trinh (2000)

$$_x^2\mathrm{\Psi }(x)+(ix^3+i\alpha x)\mathrm{\Psi }(x)=E\mathrm{\Psi }(x),$$
(1)
its Fourier transform counterpart is (i.e. $`xi_k`$, $`_xik`$)

$$k^2\widehat{\mathrm{\Psi }}(k)+(_k^3\alpha _k)\widehat{\mathrm{\Psi }}(k)=E\widehat{\mathrm{\Psi }}(k),$$
(2)
where

$$\widehat{\mathrm{\Psi }}(k)=\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^+\mathrm{}𝑑xe^{ikx}\mathrm{\Psi }(x).$$
(3)
It is important to note that, for the physical solutions, a simple application of WKB analysis (Bender and Boettcher (1998)) tells us that the asymptotic behavior of the configuration space representation yields an entire Fourier transform. Because of this, the $`k`$-power series expansion is absolutely convergent, and defined in terms of the power moments:

$$\widehat{\mathrm{\Psi }}(k)=\frac{1}{\sqrt{2\pi }}\underset{p=0}{\overset{\mathrm{}}{}}(ik)^p\frac{\mu _p}{p!},$$
(4)
where

$$\mu _p_{\mathrm{}}^+\mathrm{}𝑑xx^p\mathrm{\Psi }(x),$$
(5)
define the Hamburger power moments.

We can generate the recursion relation for the $`\mu _p`$’s from the standard power series expansion methods for linear differential equations (Bender and Orszag (1978)). Alternatively, we can apply $`_{\mathrm{}}^+\mathrm{}𝑑xx^p`$ to both sides of Eq.(1), combined with integration by parts, and obtain the necessary Moment Equation:

$$\mu _{p+3}=\alpha \mu _{p+1}iE\mu _pip(p1)\mu _{p2},$$
(6)
for $`p0`$. This corresponds to a homogeneous, linear, finite difference equation, of effective order $`1+m_s=3`$, since specification of the independent moments $`\{\mu _0,\mu _1,\mu _2\}`$, plus the (complex) energy parameter, $`E`$, generates all of the remaining moments. The independent moments are referred to as the missing moments.

The linear dependence of the moments, on the missing moments, can be expressed through the relation

$$\mu _p=\underset{\mathrm{}=0}{\overset{m_s}{}}M_{p,\mathrm{}}(E)\mu _{\mathrm{}},$$
(7)
where the energy dependent coefficients satisfy the moment equation, with respect to the $`p`$-index, as well as the initial conditions $`M_{\mathrm{}_1,\mathrm{}_2}=\delta _{\mathrm{}_1,\mathrm{}_2}`$, for $`0\mathrm{}_1,\mathrm{}_2m_s`$.

### B  Defining an Analytic Basis in the Fourier Space

The physical solutions in the Fourier representation must also be $`L^2`$. One would like to find an appropriate basis into which to transform the Fourier power series expansion:

$$\frac{1}{\sqrt{2\pi }}\underset{p=0}{\overset{\mathrm{}}{}}(ik)^p\frac{\mu _p}{p!}=\frac{1}{\sqrt{2\pi }}\underset{j=0}{\overset{\mathrm{}}{}}a_j_j(k).$$
(8)
The easiest choice, leading to a rapid, analytic generation of the $`a_j`$ coefficients, is to take

$$_j(k)=(ik)^j\widehat{R}(k),$$
(9)
where $`\widehat{R}`$ is some arbitrary “reference” function yielding a complete (if not orthogonal) basis. Thus, if $`\frac{1}{\widehat{R}(k)}`$ is analytic, then one can generate the $`a_j`$’s by expanding $`\frac{\widehat{\mathrm{\Psi }}(k)}{\widehat{R}(k)}`$. In particular, we can take $`\widehat{R}(k)=e^{\beta k^2}`$, where $`\beta >0`$, and otherwise arbitrary, yielding

$$A(k)\underset{j=0}{\overset{\mathrm{}}{}}a_j(ik)^j=e^{\beta k^2}\underset{p=0}{\overset{\mathrm{}}{}}(ik)^p\frac{\mu _p}{p!}.$$
(10)
It is then clear that the $`a_j`$’s will become linear in the missing moments, $`\{\mu _{\mathrm{}}|0\mathrm{}m_s\}`$. Specifically,

$$a_j(E;\mu _0,\mu _1,\mu _2)=\underset{p+2q=j}{}\frac{(\beta )^q\mu _p}{q!p!},$$
(11)
or

$$a_j(E;\mu _0,\mu _1,\mu _2)=\underset{\mathrm{}=0}{\overset{m_s=2}{}}\left(\underset{p+2q=j}{}\frac{(\beta )^qM_{p,\mathrm{}}(E)}{q!p!}\right)\times \mu _{\mathrm{}},$$
(12)
for $`j0`$. Clearly, the maximum (Hamburger) moment order generated, $`P_{max}`$, determines the maximum order of $`a_j`$’s generated, $`0jP_{max}`$.

### C  The MRF Quantization Prescription

It has been argued by Tymczak et al (1998a,b) that the convergent zeroes of the coefficient functions

$$a_j(E_n^{(j)})=0,$$
(13)
converge to the exact discrete state energies:

$$\underset{j\mathrm{}}{lim}E_n^{(j)}=E_n^{physical}.$$
(14)
The above root equation must be adapted to the, $`1+m_s`$, linear, missing moment structure of the Hamiltonian in question. Thus, to any expansion order $`J`$, we impose that the last $`1+m_s`$, $`a_j`$-coefficients be zero (i.e. $`a_J\mathrm{}=0`$, for $`0\mathrm{}m_s`$). This results in a $`1+m_s=3`$ dimensional, determinantal equation for the energy:

$$\mathrm{\Delta }_J(E)=0,$$
(15)
where

$$\mathrm{\Delta }_J(E)Det\left(\begin{array}{c}_{0,0}^{(J)}(E,\beta ),_{0,1}^{(J)}(E,\beta ),_{0,2}^{(J)}(E,\beta )\\ _{1,0}^{(J)}(E,\beta ),_{1,1}^{(J)}(E,\beta ),_{1,2}^{(J)}(E,\beta )\\ _{2,0}^{(J)}(E,\beta ),_{2,1}^{(J)}(E,\beta ),_{2,2}^{(J)}(E,\beta )\end{array}\right),$$
(16)
and

$$_{\mathrm{}_1,\mathrm{}_2}^{(J)}(E,\beta )\underset{p+2q=J\mathrm{}_1}{}\frac{(\beta )^qM_{p,\mathrm{}_2}(E)}{q!p!},$$
(17)
for $`0\mathrm{}_1,\mathrm{}_2m_s=2`$. In the Tables, the $`P_{max}`$ parameter corresponds to $`JP_{max}`$.

## III Numerical Implementation of MRF

The MRF analysis is implemented for $`\beta =.5`$. This value of the arbitrary parameter generates the fastest converging results.

In Figure I we plot the $`𝒫𝒯`$-symmetry breaking solutions ($`Im(E)0`$), and the $`𝒫𝒯`$-symmetry invariant solutions ($`Im(E)=0`$), for moderate $`\alpha `$ values at the limits of Delabaere and Trinh’s asymptotic analysis. There are four branches depicted. In Tables I-IV we specify some of the points plotted, and in addition, compare the MRF eigenenergy estimates with the EMM bounds. The MRF results in these tables generally correspond to $`J50`$. The EMM results were calculated to the same moment order, $`P_{max}50`$, for the complex energies. In those cases where the MRF method predicts a purely real eigenenergy, we used the faster EMM formalism corresponding to explicitly assuming that the underlying solutions are $`𝒫𝒯`$ invariant (Handy (2001a)). For this case, the underlying problem is of Stieltjes character, and a Stieltjes moment order of $`P_{max}^{(S)}=30`$, corresponds to a Hamburger moment order of 60 (i.e. $`P_{max}^{(S)}=\frac{P_{max}}{2}`$). For this reason, in some cases the MRF estimates (computed at $`J50`$) lie outside of the EMM bounds. However, for some of the EMM bounds, numerical instability concerns required that $`P_{max}^{(S)}<30`$ (these were computed on a standard, double precision, IBM platform). Such cases are relatively few in number, and are explicitly identified. By working at a larger precision order, the corresponding bounds can be improved.

It will also be noted that for $`\alpha =0`$, we quote the EMM bounds generated at higher Stieltjes moment order ($`P_{max}^{(S)}60`$), given in the work by Handy (2001a). These bounds are very consistent with the corresponding entries in Tables V-VI.

In Tables V and VI we quote the MRF results, for $`J=100`$ (using CRAY-double precision). Only the stable digits are given (that is, Tables V and VI correspond to the (empirically determined) stable digits within the MRF generated sequence, for $`0J100`$).

The results are consistent, that is, the MRF results in Tables I-IV, for the most part, lie within the bounds. This is definitely the case for Tables V - VI. Note that in the Tables we quote the imaginary parts of the energy as $`\pm `$, this is because we cannot tell which $`𝒫𝒯`$-invariant branch continues into the $`𝒫𝒯`$-breaking branch. For $`𝒫𝒯`$ invariant Hamiltonians, complex energies come in conjugate pairs. Thus, if $`E`$ is a solution, so too is $`E^{}`$.

In Figures II - III we narrow in on the smaller of the two critical $`\alpha `$ values, $`\alpha _{cr_{1,2}}`$, as shown in Fig. I. At these critical points, the energy goes from being real ($`𝒫𝒯`$ \- invariant solutions) to complex ($`𝒫𝒯`$ \- breaking solutions). They are:

$$\alpha _{cr_1}=2.611809356,$$
(18)
corresponding to $`E_{cr_1}=1.28277353562`$; and

$$\alpha _{cr_2}=5.375879629,$$
(19)
corresponding to $`E_{cr_2}=4.181388093.`$

## IV Conclusion

We have confirmed the asymptotic analysis prediction of Delabaerre and Trinh (DT) on the existence of symmetry breaking solutions for the $`H_\alpha `$ Hamiltonian. Our methods enable the precise analysis of the complex-real spectra, particularly for moderate $`\alpha `$ values at the limits of their (DT) asymptotic validity. The results of both an eigenenergy estimation method (MRF) and an eigenenergy bounding method (EMM) were presented. The algebraic simplicity, and ease of computational implementability, of the MRF method recommend it highly for application to similar problems. Through the use of readily available algebraic programming software, the MRF approach can be extended to arbitrary precision (indeed, Tymczak et al (1998b) were able to generate the quartic anharmonic oscillator ground state energy to more than 250 decimal places), making it a very powerful tool in these types of investigations.

Acknowledgments

This work was supported through a grant from the National Science Foundation (HRD 9632844) through the Center for Theoretical Studies of Physical Systems (CTSPS). The author extends his appreciation to Professors D. Bessis, G. Japaridze, and G. A. Mezincescu for useful discussions.

References

Bender C M and Boettcher S 1998 Phys. Rev. Lett. 80 5243

Bender C M, Berry M, Meisinger P N, Savage V M, and Simsek M 2001 J. Phys. A: 34 L31

Bender C M and Orszag S A, Advanced Mathematical Methods for Scientists and Engineers (New York: McGraw Hill 1978).

Delabaere E and Trinh D T 2000 J. Phys. A: Math. Gen. 33 8771

Handy C R and Brooks H 2001 J. Phys. A J. Phys. A , to appear.

Handy C R, Murenzi R, Bouyoucef K, and Brooks H 2000 J. Phys. A 33 2151

Handy C R 2001a CAU preprint (Submitted to J. Phys. A)

Handy C R 2001b CAU preprint (Submitted to J. Phys. A)

Mezincescu G A 2000 J. Phys. A: Math. Gen. 33 4911

Mezincescu G A 2001 J. Phys. A: Math. Gen. 34 3329

Shohat J A and Tamarkin J D, The Problem of Moments (American Mathematical Society, Providence, RI, 1963).

Tymczak C J, Japaridze G S, Handy C R, and Wang Xiao-Qian 1998 Phys. Rev. Lett. 80, 3678

Tymczak C J, Japaridze G S, Handy C R, and Wang Xiao-Qian 1998 Phys. Rev. A58, 2708


text:
# Horizontal-Branch Morphology as an Age Indicator

## 1. Introduction

Because the horizontal-branch (HB) stars in globular clusters are much brighter than main-sequence (MS) stars, the interpretation of HB morphology in terms of relative age differences would be of great value in the study of distant stellar populations where the MS turnoff is fainter than the detection limit. In this paper, we report our progress in the use of HB as a reliable age indicator.

## 2. New HB Population Models

Some seven years ago, Lee et al. (1994, hereafter LDZ) have concluded that age is the most natural candidate for the global second parameter, because other candidates can be ruled out from the observational evidence, while supporting evidence do exist for the age hypothesis. Although this conclusion is generally accepted in the community, there is still some debate about this issue (see Table 1 of Lee et al. 1999). Among others, critics are arguing that the relative age differences inferred from the MS turnoffs are often too small compared to the amounts predicted from the HB models for the several second parameter pairs, including the well-known M3 and M13 pair.

We found, however, several recent developments can affect the relative age dating technique from HB morphology. First of all, there is now a reason to believe that absolute age of the oldest Galactic globular clusters is reduced to about 12 Gyr, as suggested by new Hipparcos distance calibration and other improvements in stellar models (Reid 1997; Gratton et al. 1997; Chaboyer et al. 1998). As LDZ already demonstrated in their paper, this has a strong impact in the relative age estimation from HB morphology. Also, it is now well established that $`\alpha `$ elements are enhanced in halo populations (\[$`\alpha `$/Fe\] = 0.4). Finally, Reimers’(1975) empirical mass-loss law suggests more mass-loss at larger ages. The result of this effect was also presented in LDZ, but unfortunately most widely used diagram (their Fig. 7) is the one based on fixed mass-loss. We found all of the above effects make the HB morphology to be more sensitive to age (see Lee & Yoon 2001 for details). Therefore, as illustrated in Figure 1 , now the required age difference is much reduced compared to Figure 7 of LDZ. Now, only 1.1 Gyr of age difference, rather than 2 Gyr, is enough to explain the systematic shift of the HB morphology between the inner and outer halo clusters. Also, to within the observational uncertainty, age difference of about 1.5-2 Gyr is now enough to explain the observed difference in HB morphology between the remote halo clusters (Pal 3, Pal 4, Pal 14, and Eridanus) and M3. These values are consistent with the recent relative age datings both from the HST and high-quality ground-based data.

## 3. Case of M3 and M13

We obtained high quality color-magnitude data for the classic second parameter globular clusters M3 and M13. The clusters were observed during the same nights with the same instruments (MDM 2.4m), allowing us to determine the accurate relative ages. From the color difference method between the turnoff and the base of the red giant-branch (RGB), we now confirm M13 is indeed 1.7$`\pm `$0.5 Gyr older than M3 (Fig. 2; see Rey et al. 2001 for details). This is consistent with the age difference predicted from our new HB models (see Fig. 1), which provides a further support that the HB morphology is a reliable age indicator in most population II stellar systems.

## References

Chaboyer, B., Demarque, P., Kernan, P. J., & Krauss, L. M. 1998, ApJ, 494, 96

Gratton, R. G., Fusi Pecci, F., Carretta, E., Clementini, G., Corsi, C. E. & Lattanzi, M. 1997, ApJ, 491, 749

Lee, Y.-W., & Yoon, S.-J. 2001, in preparation

Lee, Y.-W., Demarque, P., & Zinn, R. 1994, ApJ, 423, 248 (LDZ)

Lee, Y.-W., Yoon, S.-J., Lee, H.-c., & Woo, J.-H. 1999, in ASP Conf. Ser. Vol. 192, Spectrophotometric Dating of Stars and Galaxies, ed. I. Hubeny, S. Heap, & Cornett (San Francisco: ASP), 185

Reid, I. N. 1998, AJ, 115, 204

Reimers, D. 1975, Mem. Soc. Roy. Sci. Leige, 6th Ser., 8, 369

Rey, S.-C., Yoon, S.-J., Lee, Y.-W., Chaboyer, B., & Sarajedini, A. 2001, submitted

VandenBerg, D. A., Bolte, M., & Stetson, P. B. 1990, AJ, 100, 445

### Acknowledgments.


text:
# A reciprocity theorem for domino tilings

## 1 Introduction

It has undoubtedly been observed many times that the Fibonacci sequence, when run backwards as well as forwards from its initial conditions, yields a doubly-infinite sequence that is symmetrical modulo some minus-signs:

$$\mathrm{},8,5,3,2,1,1,0,1,1,2,3,5,8,\mathrm{}$$

We will see that this is just one instance of a more general phenomenon involving domino tilings of rectangles. For positive integers $`m`$ and $`n`$, let $`T(m,n)`$ denote the number of ways to cover an $`m`$-by-$`n`$ rectangle with 1-by-2 rectangles (“dominos”) with pairwise disjoint interiors. When $`m=2`$, the values $`T(m,n)`$ form the Fibonacci sequence; more generally, for each fixed $`m`$, the values of $`T(m,n)`$ form a sequence satisfying a higher-order linear recurrence relation (see section 4 for an explanation of why such a recurrence relation must exist).

These recurrences allow one to extrapolate $`T(m,n)`$ to negative values of $`n`$ in a natural way, for each fixed positive integer $`m`$. There are in fact infinitely many recurrences one could employ; but, as section 4 explains, the extrapolated values of $`T(m,n)`$ do not depend on which of the recurrence relations is used.

One finds that as long as $`m`$ is not congruent to 2 (mod 4), the sequence of values $`T(m,n)`$ (with $`n`$ varying over the positive and negative integers) exhibits palindromicity:

$$T(m,2n)=T(m,n)\text{for all }n.$$

When $`m`$ is congruent to 2 (mod 4), a modified version of this holds:

$$T(m,2n)=(1)^nT(m,n)\text{for all }n.$$

These properties of the sequence of values $`T(m,n)`$ can also be expressed in terms of the generating function $`F_m(x)=_{n=0}^{\mathrm{}}T(m,n)x^n`$, with $`T(m,0)`$ taken to be 1. This generating function is a rational function of $`x`$, as a consequence of the recurrence relation governing the values $`T(m,n)`$. Denote this rational function by $`P_m(x)/Q_m(x)`$. Note that the polynomial $`Q_m(x)`$ encodes the recurrence relation that governs the sequence of coefficients $`T(m,n)`$, at least eventually; if the degree of $`P_m(x)`$ is less than the degree of $`Q_m(x)`$, then the recurrence relation applies to all the coefficients of the generating function, not just ones with $`n`$ sufficiently large.

Stanley showed St  (see parts (d) and (e) of his main Theorem) that for $`m>1`$, the degree of $`Q_m`$ is 2 more than the degree of $`P_m`$, $`Q_m(1/x)=Q_m(x)`$ when $`m2`$ (mod 4) and $`Q_m(1/x)=Q_m(x)`$ when $`m2`$ (mod 4). These reciprocity relations, first observed empirically by Klarner and Pollack KP , can be shown to be equivalent to the claims made earlier about the extrapolated values of $`T(m,n)`$.

Stanley’s proof makes use of algebraic methods; in particular, it relies upon formulas proved by Kasteleyn Ka  and Temperley and Fisher TF  that express $`T(m,n)`$ in terms of complex roots of unity (or, if one prefers, trigonometric functions).

I will present a combinatorial way of thinking about the two-sided infinite sequence $`(T(m,n))_{n=\mathrm{}}^{\mathrm{}}`$. All features of the reciprocity relation — including the factor $`ϵ_{m,n}`$ — will have direct combinatorial interpretations. Moreover, the method used here can also be used to demonstrate similar reciprocity phenomena for three- and higher-dimensional analogues of domino tiling. For instance, if $`T(k,m,n)`$ denotes the number of ways to tile a $`k`$-by-$`m`$-by-$`n`$ box with 1-by-1-by-2 blocks (in any orientation), then $`T(k,m,2n)=\pm T(k,m,n)`$, where $`T(k,m,2n)`$ is obtained by taking the recurrence relation for $`T(k,m,n)`$ (with $`k,m`$ fixed and $`n`$ varying) and running it backward. Such higher-dimensional tiling problems are not amenable to Stanley’s method, since the formulas of Kasteleyn and Temperley and Fisher do not extend beyond the two-dimensional setting.

## 2 Signed matchings

Instead of using domino-tilings, we will use an equivalent mathematical model: perfect matchings of grid-graphs, also known as dimer covers or dimer configurations111The difference between the titles of St  and KP , compared with the difference in the names of the journals in which the two articles appear, clarifies the difference between dimers and dominos: Studying dimers configurations is applied mathematics. Studying tilings by dominos is not.. A perfect matching of a graph $`𝒢`$ (hereafter called just a matching for short) is a subset of the edge-set of $`𝒢`$ such that each vertex of $`𝒢`$ belongs to one and only one edge of the subset. The domino tilings of an $`m`$-by-$`n`$ rectangle are in obvious bijection with the perfect matchings of a graph with $`mn`$ vertices arranged in $`m`$ rows of $`n`$ vertices each, with edges joining horizontal and vertical neighbors. See Figure 1.

Fig. 1: The graph $`𝒢(2,5)`$.

Given that some of the extrapolated numbers $`T(m,n)`$ are negative, it is natural to try to provide them with a combinatorial interpretation by way of a signed version of graph-matching. A signed graph is a graph whose edges are labeled $`+1`$ and $`1`$, and a matching counts as either positive or negative according to whether the number of $`1`$ edges is even or odd. (This is a special case of weighted matching, in which each edge is assigned a weight, and the weight of a matching is the product of the weights of its edges.) When we count the matchings of a signed graph, we count them according to sign. For example, the Fibonacci number $`T(2,5)=3`$ will arise as the “number” of signed matchings of the signed graph shown in Figure 2.

Fig. 2: The signed graph $`𝒢(2,5)`$.

Here the horizontal edges are all $`+1`$ edges and the vertical edges are $`1`$ edges (where the presence of a negative edge is indicated by use of dashed lines). Each of the three matchings of the graph has an odd number of vertical edges and so counts as a negative matching. Hereafter, when I refer to the number of matchings of a signed graph, I always mean the number of positive matchings minus the number of negative matchings.

By a “signed graph of width $`m`$”, I will mean a subgraph of the $`m`$-by-$`n`$ grid graph for some $`n`$, with some edges having sign 1 and others having sign $`1`$. We will be pay special attention to some specific graphs $`𝒢(m,n)`$, defined for $`m1`$ and $`n`$ an arbitrary integer. When $`n`$ is positive, $`𝒢(m,n)`$ is just the $`m`$-by-$`n`$ grid graph described above, with all edges having sign 1. When $`n`$ is less than or equal to zero, $`𝒢(m,n)`$ is a modified version of the graph $`𝒢(m,2n)`$ in which the rightmost and leftmost flanks of vertical edges have been removed, and all remaining vertical edges are decreed to have sign $`1`$. Thus Figure 1 shows the graph $`𝒢(2,5)`$ while Figure 2 shows the graph $`𝒢(2,5)`$.

It will emerge that the (signed) number of matchings of $`𝒢(m,n)`$ is $`T(m,n)`$. Note that for all $`m`$, $`𝒢(m,0)`$ and $`𝒢(m,2)`$ each have a unique matching, whose sign is $`+1`$. This verifies the reciprocity theorem in the case $`n=0`$ (and the case $`n=1`$ is trivial). To prove the claim for other values of $`n`$ (it suffices to consider $`n>0`$), we must do more work.

## 3 Adjunction of graphs

The two definitions of $`𝒢(m,n)`$ given above (one for $`n>0`$, the other for $`n0`$) are not just two distinct combinatorial notions, patched together at 0; they fit “seamlessly”. To make this clearer, we will situate the graphs $`𝒢(m,n)`$ in a broader context, in which they function as building blocks. Suppose $`𝒢_1`$ and $`𝒢_2`$ are signed graphs of width $`m`$. Define the “adjoined graph” $`𝒢_1𝒢_2`$ as the signed graph we get when we place $`𝒢_1`$ to the left of $`𝒢_2`$ and join the rightmost vertices of $`𝒢_1`$ to the respective leftmost vertices of $`𝒢_2`$ using $`n`$ edges of type $`+1`$. Figure 3 shows the graph $`𝒢(2,3)𝒢(2,3)`$.

Fig. 3: The signed graph $`𝒢(2,3)𝒢(2,3)`$.

Clearly adjunction is an associative operation, so a product of three or more signed graphs of width $`m`$ is well-defined.

Let $`M(𝒢)`$ denote the (signed) number of matchings of a signed graph $`𝒢`$. I can now state the precise sense in which the two parts of the definition of $`𝒢(m,n)`$ fit together:

$$M(𝒢(m,n_1)𝒢(m,n_2)\mathrm{}𝒢(m,n_k))$$

$$=M(𝒢(m,n_1+n_2+\mathrm{}+n_k))$$

for all $`n_1,\mathrm{},n_k`$, regardless of sign.

To indicate why the preceding equation is true, we will prove the case $`k=2`$, which contains all the essential ideas.

When $`n_1`$ and $`n_2`$ are positive, we have $`𝒢(m,n_1)𝒢(m,n_2)`$ literally equal to the graph $`𝒢(m,n_1+n_2)`$ by virtue of the definition of adjunction, so this case requires no proof.

When $`n_1`$ and $`n_2`$ are both non-positive, we get situations like the one shown in Figure 4, depicting $`𝒢(2,2)𝒢(2,3)`$.

Fig. 4: The signed graph $`𝒢(2,2)𝒢(2,3)`$.

It is a well-known lemma in the theory of matchings that the number of matchings of a graph $`𝒢`$ is unaffected if a chain of vertices $`u,v,w`$ (where $`v`$ is connected to $`u`$ and $`w`$ but to no other vertices in the graph) is shrunk down to a single vertex (i.e., $`v`$ is removed and $`u`$ and $`w`$ are replaced by a single vertex adjacent to all the vertices other than $`v`$ to which $`u`$ and $`w`$ were originally connected). The same observation applies to signed graphs as well if the edges $`uv`$ and $`vw`$ are $`+1`$ edges. Applying this lemma to $`𝒢(2,2)𝒢(2,3)`$ twice (once in the top row and once in the bottom row), we obtain a copy of the graph $`𝒢(2,5)`$. More generally, $`M(𝒢(m,a)𝒢(m,b))=M(𝒢(m,(a+b)))`$ whenever $`a,b0`$.

When $`n_1`$ and $`n_2`$ are of opposite type (one positive, the other non-positive), we get situations like the one shown in Figure 3, depicting $`𝒢(2,3)𝒢(2,3)`$. If we apply the shrinking lemma in this case, once in each row, we get the signed graph shown in Figure 5.

Fig. 5: $`𝒢(2,3)𝒢(2,3)`$, after shrinking.

Here the edge marked $`+`$ is actually two edges, one a $`+1`$ edge and the other a $`1`$ edge. Any signed matching that uses the $`+1`$ edge can be paired with a signed matching that uses the $`1`$ edge; these matchings have opposite sign, and so together contribute 0 to the signed number of matchings. Hence the number of matchings of the graph is unaffected if we remove both of these edges, obtaining the signed graph shown in Figure 6.

Fig. 6: $`𝒢(2,3)𝒢(2,3)`$, reduced to $`𝒢(2,2)𝒢(2,2)`$.

But this is just $`𝒢(2,2)𝒢(2,2)`$. More generally, applying the “shrinking lemma” to $`𝒢(m,a)𝒢(m,b)`$ yields $`𝒢(m,a1)𝒢(m,(b1))`$ as long as $`a,b`$ are both positive. This sets the stage for a proof by induction, and all that remains is to verify the base case where $`a=0`$ or $`b=0`$. This is an easy verification that I leave to the reader. The case of $`𝒢(m,a)𝒢(m,b)`$ is identical.

## 4 Recurrence relations

We start by showing that for every width-$`m`$ signed graph $``$, the sequence $`M(𝒢(m,n))`$ (with $`n=1,2,3,\mathrm{}`$) satisfies a linear recurrence relation of degree at most $`2^m`$ that depends only on $`m`$, not on $``$.

Let $``$ be a width-$`m`$ signed graph. For any sequence $`c_1,\mathrm{},c_m`$ of $`m`$ bits, and for $`n1`$, let $`S(n;c)`$ be the sum of the signs of the matchings of the signed graph obtained from $`𝒢(m,n)`$ by deleting the subset of the rightmost $`m`$ vertices specified by the bit-pattern $`c`$; specifically, we remove the rightmost vertex in row $`i`$ (and all edges incident with it) for precisely those $`i`$ with $`c_i=0`$. Let $`S(n)`$ be the row-vector of length $`2^m`$ whose components are the numbers $`S(n;c)`$. Then (cf. KP  and Re ) there is a $`2^m`$-by-$`2^m`$ matrix $`M`$ such that $`S(n+1)=S(n)M`$. More specifically, the $`c,c^{}`$ entry of $`M`$ gives the number of ways in which each matching counted by $`S(n;c)`$ can be extended to a matching counted by $`S(n+1;c^{})`$, adding only horizontal edges between the $`n`$th and $`n+1`$st columns and vertical edges in the $`n+1`$st column. Hence we have $`S(n)=S(1)M^{n1}`$. By the Cayley-Hamilton Theorem, the sequence of powers of $`M`$ satisfies the linear recurrence associated with the characteristic polynomial of $`M`$, and it follows that the sequence of vectors $`S(1),S(1)M,S(1)M^2,\mathrm{}`$ also satisfies this recurrence. In particular, the entry associated with $`c=(1,1,\mathrm{},1)`$ satisfies this recurrence.

For instance, with $`m=2`$, $`M`$ is the 4-by-4 matrix

$$\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 1\end{array}\right),$$

where rows and columns are indexed in the order $`(0,0),(0,1),(1,0),(1,1)`$. The characteristic polynomial of this matrix is $`x^4x^32x^2+x+1`$, so the one-sided sequence of numbers $`a_n=M(𝒢(m,n))`$ (with $`n`$ going from 1 to $`\mathrm{}`$) must satisfy the recurrence $`a_{n+4}a_{n+3}2a_{n+2}+a_{n+1}+a_n=0`$, regardless of the nature of the width-$`m`$ signed graph $``$.

I now claim that the two-sided sequence of numbers $`a_n`$ (with $`n`$ going from $`\mathrm{}`$ to $`\mathrm{}`$) satisfies this same recurrence relation. If we take $`=𝒢(m,N)`$ with $`N`$ large, and use the fact that $`M(𝒢(m,N)𝒢(m,n))=M(𝒢(m,nN))=T(m,nN)`$, we see that if we take the doubly-infinite sequence $`\mathrm{},T(m,1),T(m,0),T(m,1),\mathrm{}`$ and start it at $`T(m,N)`$, we get a singly-infinite sequence that (as we have shown) must satisfy the Cayley-Hamilton recurrence. Since $`N`$ was arbitrary, the two-sided infinite sequence of $`a_n`$’s satisfies the Cayley-Hamilton recurrence as well.

Having constructed one recurrence relation satisfied by the doubly-infinite sequence $`𝐯=(\mathrm{},T(m,1),T(m,0),T(m,1),\mathrm{})`$, or equivalently, one linear operator $`A`$ that annihilates $`𝐯`$, we must now consider others. Could there be a linear operator $`B`$ that annihilates the singly-infinite sequence $`T(m,1),T(m,2),T(m,3),\mathrm{}`$ but not the doubly-infinite sequence $`𝐯`$? If so, then $`B`$ would yield a different extrapolation of $`T(m,n)`$ to negative values of $`n`$ than $`A`$ would.

To show that no such operator $`B`$ exists, note that $`BA`$ annihilates $`𝐯`$. Since $`BA=AB`$, $`B`$ must send $`𝐯`$ into something annihilated by $`A`$. If $`B`$ does not itself annihilate $`𝐯`$, then $`B`$ sends $`𝐯`$ into some other sequence that vanishes for all sufficiently large indices but does not vanish for all indices. It is easy to show that such a sequence cannot be annihilated by $`A`$ or indeed by any linear operator with constant coefficients.

The above arguments show that the numbers $`T(m,n)`$, defined for non-positive values of $`n`$ as above, are the unique way to extrapolate so that a linear recurrence is satisfied.

## 5 Reciprocity

With all the ingredients in place, we can now give a very simple explanation of the reciprocity relation for domino tilings of rectangles of fixed width. For concreteness, we start with the special case $`T(2,5)`$. Referring to Figure 2, we see that every matching of $`𝒢(2,5)`$ must contain the two leftmost horizontal edges and the two rightmost horizontal edges. If we remove these edges from the graph (and the vertices incident with those edges, and the edges incident with those vertices), we get a copy of $`𝒢(2,3)`$ in which the vertical edges have sign 1. Since every matching of $`𝒢(2,3)`$ has an odd number of vertical edges, we see that $`T(2,5)=T(2,3)`$.

More generally, when $`n`$ is positive, every matching of $`𝒢(m,2n)`$ must contain the $`m`$ leftmost horizontal edges and the $`m`$ rightmost horizontal edges. Removing these edges (and concomitant vertices and edges) leaves a copy of $`𝒢(m,n)`$ in which the vertical edges have sign 1. It is well-known (see e.g. Th ) that every domino-tiling of a rectangle (or indeed more general regions) can be obtained from every other by means of moves in which a 2-by-2 block of horizontal dominos is rotated to give a 2-by-2 block of vertical dominos, or vice versa. Since moves of this kind do not affect the parity of the number of vertical dominos (or, in our terms, the parity of the number of vertical edges in a matching), we know that that all the signed matchings of our modified version of $`𝒢(m,n)`$ carry the same sign; it remains to determine what this sign is.

If one is willing to appeal to the invariance-of-parity result mentioned in the preceding paragraph, it is simple to evaluate the common sign of all the matchings by considering the all-horizontal matching or the all-vertical matching. But in the interests of making the article self-contained, I give a direct argument (of a fairly standard kind). Say we have a matching of $`𝒢(m,n)`$ that involves $`k_i`$ vertical edges joining row $`i`$ and row $`i+1`$, for $`i`$ ranging from 1 to $`m1`$. When these edges and their vertices are removed from the graph, along with all other edges joining row $`i`$ and row $`i+1`$, our graph splits into two subgraphs, each of which must have an even number of vertices (since we are assuming that we have a matching). Hence $`nik_i`$ and $`n(mi)k_i`$ must be even. Using the congruence $`k_ini`$ (mod 2), we find that $`k_1+\mathrm{}+k_{m1}`$, the number of vertical edges, must be congruent to $`n(1+2+\mathrm{}+(m1))=n(m1)m/2`$ (mod 2). If $`m`$ is congruent to 0 or 1 (mod 4), $`(m1)m/2`$ is even, so the number of vertical edges is even. If $`m`$ is congruent to 3 (mod 4), then (in order for a matching to exist) $`n`$ must be even, so that $`n(m1)m/2`$ is even. However, if $`m`$ is congruent to 2 (mod 4), then $`(m1)m/2`$ is odd, so that $`n(m1)m/2`$ is even or odd according to the parity of $`n`$.

Hence, if we define $`ϵ_{m,n}`$ to be $`1`$ when $`m`$ is congruent to 2 (mod 4) and $`n`$ is odd, and $`+1`$ otherwise, we have

$$T(m,2n)=ϵ_{m,n}T(m,n)$$

as claimed.

## 6 Motivation

The definition of $`𝒢(m,n)`$ for $`n<0`$ arose not by ad hoc insight but by a particular vision of combinatorics as a specialized form of algebra. Under this vision, a collection of combinatorial objects ought to be represented by a multivariate polynomial whose coefficients are all equal to 1, where the individual terms represent the combinatorial objects themselves in some fashion.

In the case of matchings of $`𝒢(2,n)`$ with $`n>0`$, we assign a formal variable (or “weight”) to each edge, define the weight of a matching as the product of the weights of its constituent edges, and define the polynomial $`𝐌(𝒢(m,n))`$ as the sum of the weights of all the matchings. This is a polynomial in which every coefficient equals 1, and in which the constituent monomials encode the respective matchings of the graph. If we think of $`𝒢(2,n)`$ as being embedded in $`𝒢(2,n+1)`$ as the induced subgraph on the leftmost $`n`$ columns, we find that the edge-variables occurring in $`𝐌(𝒢(2,n))`$ form a subset of the edge-variables occurring in $`𝐌(𝒢(2,n+1))`$, and that indeed, there is a recurrence expressing $`𝐌(𝒢(2,n+1))`$ in terms of $`𝐌(𝒢(2,n))`$ and $`𝐌(𝒢(2,n1))`$:

$$𝐌(𝒢(2,n+1))=y_{n+1}𝐌(𝒢(2,n))+w_nx_n𝐌(𝒢(2,n1)),$$

where $`y_{n+1}`$ is the weight of the rightmost vertical edge of $`𝒢(2,n+1)`$ and $`w_n`$ and $`x_n`$ are the weights of the rightmost horizontal edges of $`𝒢(2,n+1)`$.

We can run this recurrence backward, obtaining in succession some rational functions of the variables $`w_n,x_n,y_n`$ with $`n`$ negative. These rational functions are in fact Laurent polynomials; moreover, all the coefficients of these Laurent polynomials are $`\pm 1`$. Just as the monomials of $`𝐌(𝒢(2,n))`$ with $`n>0`$ encode the matchings of $`𝒢(2,n)`$, one might assume that the monomials in the extrapolated Laurent polynomials “$`𝐌(𝒢(2,n))`$” with $`n0`$ encode combinatorial objects that will yield a combinatorial interpretation of the extrapolated values $`T(2,n)`$ with $`n0`$. The monomials lead one to matchings of the graphs $`𝒢(2,n)`$ with $`n<0`$, and once one has this signed-matching model, the idea is readily extended to $`T(m,n)`$ for all $`m`$.

It is worth remarking that the reciprocity theorem satisfied by the Laurent polynomials $`𝐌(𝒢(m,n))`$ is a bit more complicated than the reciprocity theorem satisfied by the numbers $`T(m,n)`$; in addition to the sign-factor $`ϵ_{m,n}`$, there is also a monomial factor involving the variables. A comprehensive reciprocity theorem for domino tilings of rectangles would take this into account. I intend to write a longer article that treats these issues, and also applies them to other models, such as matchings of honeycomb hexagons (see CLP ). I expect that algebraic recurrences for sequences of Laurent polynomials will also shed light on the issues raised in Problem 32 (the last problem) in Pr .