text:
# Thermodynamics of spin 𝑺=𝟏/𝟐 antiferromagnetic uniform and alternating-exchange Heisenberg chains

## I Introduction

An antiferromagnetic alternating-exchange Heisenberg chain is one in which nearest-neighbor spins in the chain interact via a Heisenberg interaction, but with two antiferromagnetic (AF) exchange constants $`J_2J_1,J_1,J_20`$ which alternate from bond to bond along the chain; the alternation parameter is $`\alpha J_2/J_1`$. Here we will be concerned with the magnetic susceptibility $`\chi `$ and specific heat $`C`$ versus temperature $`T`$ of alternating-exchange chains consisting of spins $`S=1/2`$. The uniform AF Heisenberg chain is one limit of the alternating chain in which the two exchange constants are equal ($`\alpha =1,J_1=J_2J`$). At the other limit is the isolated dimer in which one of the exchange constants is zero ($`\alpha =0`$). The present work is a combined theoretical and experimental study of $`\chi (T)`$ and $`C(T)`$ of the $`S=1/2`$ AF alternating-exchange chain over the entire range $`0\alpha 1`$ of the alternation parameter, with the emphasis on the regime $`\alpha 1`$ at and close to the uniform chain limit. This latter regime is relevant for compounds showing second order spin dimerization transitions with decreasing $`T`$. The present work was originally motivated by our desire to accurately extract the temperature dependent energy gap $`\mathrm{\Delta }(T)`$ for magnetic excitations, the “spin gap”, from experimental $`\chi (T)`$ data for the $`S=1/2`$ chain/two-leg ladder compound $`\mathrm{NaV}_2\mathrm{O}_5`$ below its spin dimerization temperature $`T_\mathrm{c}34`$ K. We found that existing theory for the alternating-exchange chain was insufficient to accomplish this goal. In the present work we critically examine the predictions of previous theory, perform the required additional theoretical calculations, and then apply the results to extract $`\mathrm{\Delta }(T)`$ at $`TT_\mathrm{c}`$ from our $`\chi (T)`$ data for $`\mathrm{NaV}_2\mathrm{O}_5`$ single crystals. We have extended the original goal so that we also include theoretical and experimental studies of $`C(T)`$ and how this quantity relates to $`\chi (T)`$. In the remainder of this introduction we briefly review the prior theoretical results pertaining to $`\chi (T)`$ and $`C(T)`$ of the uniform and alternating-exchange chain to place our work in the proper context. We then review the experimental and theoretical background on $`\mathrm{NaV}_2\mathrm{O}_5`$ and describe the plan for the rest of the paper.

### A Theory

The $`\chi (T)`$ and $`C(T)`$ of both limits of the $`S=1/2`$ AF alternating-exchange Heisenberg chain are known exactly. For the dimer, the $`\chi (T)`$ is given by the exact Eq. (14) below and the exact $`C(T)`$ is also easily calculated. The $`\chi (T)`$ and $`C(T)`$ of the uniform chain for $`T0.4J/k_\mathrm{B}`$ ($`k_\mathrm{B}`$ is Boltzmann’s constant) were estimated from calculations for chains with $`11`$ spins by Bonner and Fisher in 1964; they extended their results by extrapolating to $`T=0`$, and in the case of $`\chi (T)`$ to the exact $`T=0`$ value. The exact solution for $`\chi (T)`$ of the uniform chain was obtained using the Bethe ansatz in 1994 by Eggert, Affleck and Takahashi, and compared with their low-$`T`$ results from conformal field theory. They found, remarkably, that $`\chi (T0)`$ has infinite slope. Their numerical $`\chi (T)`$ values are up to $`10`$% larger than the Bonner-Fisher extrapolation for $`T0.25J/k_\mathrm{B}`$ (for a comparison of the two predictions, see Fig. 8.1 in Ref. ). Their conformal field theory calculations showed that the leading order correction to the zero temperature limit is of the form $`\chi (T)=\chi (0)\{1+1/[2\mathrm{ln}(T_0/T)]\}`$, where the value of the scaling temperature $`T_0`$ is not predicted by the field theory. Such log terms are called “logarithmic corrections” in the literature. One of us recently presented numerical Bethe ansatz calculations of $`\chi (T)`$ with a higher absolute accuracy for $`\chi (T)`$ estimated to be $`1\times 10^7`$, and showed that the data are consistent with the above field theory prediction, with an additional higher order logarithmic correction, over the temperature range $`5\times 10^{25}k_\mathrm{B}T/J10^3`$. Corresponding $`C(T)`$ calculations were also carried out, and logarithmic corrections were studied for this quantity as well. Lukyanov has recently presented an exact theory for $`\chi (T)`$ and $`C(T)`$ at low $`T`$, including the exact value of $`T_0`$. In the present work, we compare the very recent numerical Bethe ansatz results of Klümper and Johnston with the predictions of Lukyanov’s theory and find agreement for $`\chi (T)`$ to high accuracy ($`1\times 10^6`$) over a temperature range spanning 18 orders of magnitude, $`5\times 10^{25}k_\mathrm{B}T/J5\times 10^7`$; the agreement in the lower part of this temperatures range is much better, $`𝒪(10^7)`$. For $`C(T)`$, the logarithmic correction in Lukyanov’s theory is insufficient to describe the Bethe ansatz data sufficiently accurately even at very low temperatures, so we derive the next two logarithmic corrections from the Bethe ansatz $`C(T)`$ data. For various applications, it would be desirable to have fits to the $`\chi (T)`$ and $`C(T)`$ Bethe ansatz data which extend to higher temperatures. We describe the formulation and implementation of fit functions, incorporating the influence of the logarithmic correction terms, which yield extremely precise fits to the data for both quantities over the entire 25 decades in temperature of the calculations, $`5\times 10^{25}k_\mathrm{B}T/J5`$.

The $`\chi (T)`$ in the intermediate regime $`0<\alpha <1`$ has been investigated analytically in the Hartree-Fock approximation and using numerical techniques. Of particular interest here is the regime $`\alpha 1`$, close to the uniform limit, which is the regime relevant to materials exhibiting a dimerization transition with decreasing $`T`$ such as occurs in materials exhibiting a spin-Peierls transition. There are no accurate theoretical predictions available for $`\chi (T)`$ of the alternating-exchange Heisenberg chain in this regime, which is the property usually used to initially characterize the occurrence of such a transition experimentally. To address this deficiency and to also cover a more extended $`\alpha `$ range, we carried out extensive quantum Monte Carlo (QMC) simulations and transfer-matrix density-matrix renormalization group (TMRG) calculations of $`\chi (T)`$ for $`0.05\alpha 1`$ over the temperature range $`0.002k_\mathrm{B}T/J_110`$.

An interesting issue is how the spin gap $`\mathrm{\Delta }`$ evolves with alternation parameter $`\alpha `$ as the uniform limit is approached, $`\alpha 1`$. Because the uniform chain is a gapless quantum-critical system, the introduction of alternating exchange along the chain has been theoretically predicted to yield a nonanalytic $`\mathrm{\Delta }(\alpha )`$ behavior for $`\alpha 1`$. We derive $`\mathrm{\Delta }(\alpha )`$ by fitting our low-$`t`$ TMRG $`\chi (T)`$ data by an expression which we formulated. The $`\mathrm{\Delta }(\alpha )`$ results are compared with those of previous numerical calculations and with the theoretical prediction. We infer from our data that the asymptotic critical regime is only entered for $`\alpha 0.99`$.

In order to be optimally useful for accurately modeling experimental $`\chi (T)`$ data for alternating-exchange chain compounds, our QMC and TMRG $`\chi (\alpha ,T)`$ results must first be accurately fitted by a continuous function of both $`\alpha `$ and $`T`$. We will introduce a general fit function which eventually proves capable of fitting these combined data for the alternating-exchange Heisenberg chain very accurately. We first fit the $`\chi (T)`$ of the uniform chain and isolated dimer using this function and then use the obtained fitting parameters as end-point parameters in the fit to our combined QMC and TMRG data for intermediate values of $`\alpha `$. The final fit function is a single two-dimensional function of $`\alpha `$ and $`T`$ for $`0\alpha 1`$ which can be used to extract the (possibly temperature-dependent) alternation parameter, exchange constants and spin gap from experimental $`\chi (T)`$ data for compounds for which the $`S=1/2`$ AF alternating-exchange Heisenberg chain Hamiltonian is appropriate. Our fit function will also be useful as a reference for $`\chi (T)`$ calculated from other related $`S=1/2`$ Hamiltonians such as that incorporating the spin-phonon interaction for spin-Peierls systems.

### B NaV2O5

Vanadium oxides show a remarkable variety of electronic behaviors. For example, the metallic fcc normal-spinel structure compound LiV2O4 shows local momentlike behaviors above $`50`$ K, crossing over to heavy fermion behaviors below $`10`$ K. On the other hand, the $`d^1`$ compound CaV2O5 has a two-leg trellis-ladder-layer structure in which all of the V atoms are crystallographically equivalent and is a Mott-Hubbard insulator. The $`\chi (T)`$ shows a spin-gap $`\mathrm{\Delta }/k_\mathrm{B}660`$ K arising from strong coupling of the V $`S=1/2`$ spins across a rung. Modeling of $`\chi (T)`$ by QMC simulations confirmed that this compound consists magnetically of V2 dimers, with an intradimer AF exchange constant $`J/k_\mathrm{B}665`$ K and with very weak interdimer interactions.

The compound NaV2O5 can also be formed. The crystal structure was initially reported in 1975 to consist of two-leg ladders as in CaV2O5, but in a non-centrosymmetric (acentric) structure (space group P21mn) in which charge segregation occurs such that one leg of each ladder consists of V+4 and the other of crystallographically inequivalent V+5 ions. However, recently five different crystal structure investigations showed that the structure is actually centrosymmetric (space group Pmmn), with all V atoms crystallographically equivalent at room temperature, so that (static) charge segregation between the V atoms does not, in fact, occur. This result is consistent with 51V NMR investigations which showed the presence of only one type of V atom at room temperature. This compound is thus formally a mixed-valent $`d^{0.5}`$ system, which has been considered in a one-electon-band picture to be a quarter-filled ladder compound. We note that from modeling optical excitations in the energy range 4 meV–4 eV, Damascelli and coworkers initially concluded that the room-temperature structure of NaV2O5 is acentric; their analysis was consistent with the V atoms on a rung of a ladder having oxidation states of 4.1 and 4.9, respectively. However, this group subsequently explained that length- and/or time-scale-of-measurement issues may be involved in their interpretation, such that charge disproportionation between V atoms may only occur locally and possibly dynamically, which could then be consistent with the (average long-range) crystal structure refinements and NMR measurements. Theoretical support for this scenario was provided by Nishimoto and Ohta. Factor group analyses of the possible IR- and Raman-active phonon modes and comparisons with experimental observations at room temperature are consistent with the centrosymmetric space group for the compound. A first-principles electronic structure study based on the density functional method within the generalized gradient approximation showed that the total energy of the centric structure is about 1.0 eV/(formula unit) lower than that of the acentric structure, consistent with the recent structural studies.

One might expect that the hole-doping which occurs upon replacing Ca in CaV2O5 by Na would result in metallic properties for NaV2O5, because of the nonintegral oxidation state of the V cations and of the crystallographic equivalence of these atoms. However, NaV2O5 is a semiconductor. This has been explained by the formation of $`d^1`$ V-O-V molecular clusters on the rungs of the two-leg ladders, again resulting in a Mott-Hubbard insulator due to the on-site Coulomb repulsion, where in this case a “site” is a V-O-V molecular cluster. Thus a nonintegral oxidation state and crystallographic equivalence of transition metal atoms in a compound are not sufficient to guarantee metallic character simply by symmetry; all nearest-neighbor pairs, triplets, …, of transition metal atoms must also be crystallographically equivalent, which is not the case in NaV2O5, since a V2 pair on a rung is not crystallographically equivalent to one on a leg in the two-leg ladders. In contrast, all V atoms and pairs of V atoms in mixed-valent fcc LiV2O4 are respectively crystallographically equivalent, resulting in metallic character as demanded by symmetry.

The V-O-V rung molecular clusters which are coupled along the ladder direction in NaV2O5 may be considered to form an effective $`S=1/2`$ one-dimensional (1D) chain. Experimental support for this picture, often quoted in the literature, is that the magnetic susceptibility (above $`T_\mathrm{c}`$, see below) is in agreement with the Bonner-Fisher prediction for the $`S=1/2`$ Heisenberg chain, as reported by Isobe and Ueda. Angle-resolved photoemission spectroscopy (ARPES) measurements on NaV2O5 by Kobayashi et al. showed that the electronic structure is essentially 1D, despite the ostensibly 2D nature of the trellis layer, with dispersion in the oxygen and copper bands (below the Fermi energy) occurring only in the ladder direction ($`b`$-axis). Interestingly, the dispersion in the lowest binding energy part of the occupied Cu Hubbard band showed a lattice periodicity of 2$`b`$, which may reflect dynamical short-range AF and/or crystallographic ordering in the ladder direction. Temperature-dependent ARPES measurements on Na0.96V2O5 by the same group from 120 to 300 K showed evidence for the predicted spin-charge separation in 1D magnetic systems.

A phase transition occurs in NaV2O5 at a critical temperature $`T_\mathrm{c}33`$–36 K, below which the spin susceptibility $`\chi ^{\mathrm{spin}}0`$ as $`T0`$ and a lattice distortion occurs. The lattice distortion results in a supercell with lattice parameters $`2a\times 2b\times 4c`$. Therefore the transition was initially characterized as a possible spin-Peierls transition, which by definition is driven by magnetoelastic (spin-phonon) coupling, and in which there is no change in the charge/spin distribution within the rungs/V-O-V molecular clusters. The superstructure in the $`a`$ and $`c`$ directions, perpendicular to the V chains which run in the $`b`$ direction, would be a result of the phasing of the distortions in adjacent chains/ladders. In this interpretation, and within the adiabatic approximation (discussed later), one would expect that the magnetic properties above $`T_\mathrm{c}`$ should be close to those of the $`S=1/2`$ Heisenberg uniform chain, and of an $`S=1/2`$ alternating-exchange Heisenberg chain below $`T_\mathrm{c}`$.

It has become clear, however, that the phase transition occurring at $`T_\mathrm{c}`$ in NaV2O5 is accompanied by charge ordering, in contrast to a classic spin-Peierls transition. Therefore, the magnetoelastic coupling may only play a secondary role, and the spin gap may be a secondary order parameter. In particular, 51V NMR experiments showed the presence of (inequivalent) V+4 and V+5 below $`T_\mathrm{c}`$, whereas only one V species was present above $`T_\mathrm{c}`$. This result is consistent with the solution of the superstructure below $`T_\mathrm{c}`$ by Lüdecke and co-workers using synchrotron x-ray diffraction. Lüdecke et al. found that there are modulated and unmodulated chains of V atoms below $`T_\mathrm{c}`$, tentatively assigned to magnetic and nonmagnetic chains. One interpretation of the results is that the $`d^1`$ V+4 cations segregate into alternate two-leg ladders which are isolated from each other within the $`\mathrm{V}_2\mathrm{O}_3`$ trellis layer by intervening two-leg ladders containing only nonmagnetic V+5. The anomalous strong increase in the thermal conductivity below $`T_\mathrm{c}`$ may also be due to charge ordering. From ultrasonic measurements of shear and longitudinal elastic constants, Schwenk and co-workers have suggested that the charge ordering is of the zig-zag type within each ladder. In each of these scenarios for charge ordering, static charge disproportionation occurs such that 1/2 of the V atoms have oxidation state +4 and the other half +5, consistent with the average formal oxidation state of +4.5 in the compound.

Köppen et al. have concluded from thermal expansion measurements that the phase transition at $`T_\mathrm{c}`$ actually consists intrinsically of two closely-spaced phase transitions separated by $`1`$ K, where the upper transition is thermodynamically of second order whereas the lower one is first order. However, a double transition was not found in their specific heat measurements on the same crystal, which they attributed to the 50 mK temperature oscillations required by their ac measurement technique which were thought to broaden the two transitions and render them indistinguishable.

The nature of the possible charge ordering pattern has been studied theoretically by several groups. Seo and Fukuyama predicted (at $`T=0`$) a static zig-zag chain of V+4 ions on each two-leg ladder, with an interpenetrating zig-zag chain of V+5 ions. They proposed that pairs of V+4 spins, one each on adjacent ladders, would form spin singlets, resulting in the observed spin gap. A similar zig-zag charge configuration in each ladder was inferred by Mostovoy and Khomskii, with subsequent experimental support by Smirnov et al., and by Gros and Valenti. Motivated in part by the above thermal expansion measurement results of Köppen et al., Thalmeier and Fulde proposed that the charge ordering transition would result in one linear chain of V+4 and one linear chain of V+5 on each two-leg ladder, thereby then allowing a conventional spin-Peierls transition to occur at a slightly lower temperature, resulting in a double transition as reported by Köppen et al. A similar picture was put forward by Nishimoto and Ohta. Thalmeier and Yaresko have extensively discussed the linear-chain and zig-zag scenarios for charge ordering, and in addition have considered the alternating two-leg ladder charge ordering pattern of the type suggested by Lüdecke et al. They point out that in both the linear and zig-zag patterns, a secondary spin-Peierls dimerization or spin exchange anisotropy (in spin space) may be necessary to give a spin gap, whereas the two-leg ladder ordering has a spin gap even with no lattice distortion. Thalmeier and Yaresko describe the characteristic signatures of each of the charge-ordered models to be compared with experimental inelastic neutron scattering measurements. Riera and Poilblanc have discussed the influence of electron-phonon coupling on the derived charge- and spin-order phase diagrams.

We have carried out $`\chi (T)`$ measurements from 2 to 750 K on single crystals of NaV2O5 along the ladder ($`b`$ axis) direction to further characterize and clarify the nature of the magnetic interactions and ordering below and above $`T_\mathrm{c}`$. We find that the magnetic properties above $`T_\mathrm{c}`$ are not accurately described by the $`S=1/2`$ Heisenberg uniform chain prediction with a $`T`$-independent $`J`$, although a mean-field ferromagnetic interchain coupling can explain these data. Using our theoretical $`\chi (\alpha ,T)`$ fit function for the AF alternating-exchange chain below $`T_\mathrm{c}`$, we find that $`\delta (0)(1\alpha )/(1+\alpha )=0.034(6)`$ and that the zero-temperature spin-gap of NaV2O5 is $`\mathrm{\Delta }(0)/k_\mathrm{B}=103(2)`$ K. The $`\delta (T)`$ and $`\mathrm{\Delta }(T)`$ below $`T_\mathrm{c}`$ are extracted. A spin pseudogap is found to occur above $`T_\mathrm{c}`$ with a rather large magnitude. From our specific heat measurements on two crystals, we find that the magnetic specific heat at low temperatures $`T15`$ K is too small to be resolved experimentally, and that the spin entropy at $`T_\mathrm{c}`$ is too small to account for the entropy of the transition. A quantitative analysis shows that at least 77 % of the entropy change at $`T_\mathrm{c}`$ due to the transition(s) and associated order parameter fluctuations must arise from the lattice and/or charge degrees of freedom and less than 23 % from the spin degrees of freedom.

### C Plan of the paper

The rest of the paper is organized as follows. Our notation for the Heisenberg spin Hamiltonian and for the reduced susceptibility, temperature and spin gap are given immediately in Sec. II. Some general features of the high-temperature series expansion (HTSE) for $`\chi (T)`$ and $`C(T)`$ of $`S=1/2`$ Heisenberg spin lattices and the low-temperature limits of these quantities for one-dimensional (1D) systems with a spin gap are then given. We then specialize to the $`S=1/2`$ AF alternating-exchange Heisenberg chain in Sec. II C, where we discuss the HTSEs, the spin gap and the one-magnon dispersion relations $`E(\mathrm{\Delta },k)`$. In the latter subsection, we derive a one-parameter approximation for $`E(\mathrm{\Delta },k)`$ which correctly extrapolates to the $`\alpha 0`$ limit and which we will need in order to later fit the TMRG $`\chi (T)`$ data to extract $`\mathrm{\Delta }(\alpha )`$. We also show that the expressions for the low-$`T`$ limits of both $`\chi (T)`$ and $`C(T)`$ depend only on the spin gap (in addition to $`T`$). In Sec. III, we discuss overall features of the $`\chi (T)`$ and $`C(T)`$ for the uniform chain and then focus on the low-$`T`$ behavior. The explicit forms of the logarithmic corrections previously found for $`\chi (T)`$ are first discussed. We show that a low-$`T`$ expansion of the theory of Lukyanov gives the same first two corrections, and in addition gives the next higher-order term. We then compare the Bethe ansatz $`\chi (T)`$ results directly with the theory with no adjustable parameters or approximations. Logarithmic corrections are also found to be important to accurately describe the Bethe ansatz data for $`C(T)`$. We show that the lowest order correction is not sufficient to fit the data, and we derive the next two higher-order corrections by fitting the data at very low temperatures.

General features of our scheme to fit numerical $`\chi (T)`$ data are described in Sec. IV A, followed by a fit to the exact $`\chi (T)`$ for the antiferromagnetic Heisenberg dimer and two fits to the numerical $`\chi (T)`$ data for the uniform chain. Due to the special requirements of, and constraints on, the two-dimensional fit function necessary to accurately fit $`\chi (\alpha ,T)`$ data for the alternating-exchange chain over large ranges of both $`\alpha `$ and $`T`$, a separate section, Sec. IV E, is devoted to formulating and discussing this fit function. Using a fit function similar to that used to fit numerical $`\chi (T)`$ data, in the next section an extremely accurate and precise fit is obtained over 25 decades in temperature to the Bethe ansatz $`C(T)`$ data. Our QMC and TMRG $`\chi (T)`$ data for the alternating-exchange chain are presented and fitted in Sec. V, using as end-point parameters those determined for the uniform chain and the dimer, respectively. The spin gap $`\mathrm{\Delta }(\alpha )`$ is extracted for $`0.8\alpha 0.995`$ by fitting the TMRG $`\chi (\alpha ,T)`$ data at low temperatures in Sec. VI. Section VII contains a comparison of our numerical results with previous work. The $`\mathrm{\Delta }(\alpha )`$ values are compared with previous numerical results and with the theoretical prediction for the asymptotic critical behavior in Sec. VII A. Our $`\chi (T)`$ calculations are shown in Sec. VII B to be in good agreement with the previous numerical results of Barnes and Riera for $`0.2\alpha 0.8`$. The numerical calculations of Bulaevskii have been extensively used in the past by experimentalists to fit the $`\chi (T)`$ of spin-Peierls compounds, but up to now a detailed analysis of the predictions of this theory has not been given. We present such an analysis in Sec. VII C and compare our results with these predictions.

We begin the experimental part of the paper by studying the anisotropic magnetic susceptibility of a high quality NaV2O5 single crystal in Sec. VIII A, where literature data on the anisotropy of the $`g`$ factor and Van Vleck susceptibility are compared with our results. In the following sections we illustrate the utility and application of many of the theoretical results derived and presented previously in the paper. In Sec. VIII B we present experimental $`\chi (T)`$ data for single crystals of NaV2O5 and model these data in detail in Sec. VIII C using our QMC and TMRG $`\chi (T)`$ data fit function for the AF alternating-exchange Heisenberg chain. We show that qualitatively and quantitatively new information about the temperature dependences of the spin dimerization parameter and spin gap below $`T_\mathrm{c}`$ can be obtained from our modeling. This analysis also shows that spin dimerization fluctuations and a spin pseudogap are present above $`T_\mathrm{c}`$, and we quantitatively determine their magnitudes. Our specific heat measurements of NaV2O5 single crystals and our extensive modeling of these data are presented in Sec. VIII D, where we obtain quantitative limits on the relative contributions of the lattice, spin and charge degrees of freedom to the change in the entropy due to the transition at $`T_\mathrm{c}`$ and to associated order parameter fluctuations. A summary and concluding discussion are given in Sec. IX.

## II Theory

In this paper we will only be concerned with the spin $`S=1/2`$ antiferromagnetic (AF) Heisenberg Hamiltonian

$$=\underset{<ij>}{}J_{ij}𝑺_i𝑺_j,$$
(1)
where $`J_{ij}`$ is the Heisenberg exchange interaction between spins $`𝑺_i`$ and $`𝑺_j`$ and the sum is over unique exchange bonds. A $`J_{ij}>0`$ corresponds to AF coupling, whereas $`J_{ij}<0`$ refers to ferromagnetic coupling. Note that magnetic nearest neighbors $`𝑺_j`$ of a given spin $`𝑺_i`$ in Eq. (1) need not be crystallographic nearest neighbors. A magnetic nearest neighbor of a given spin is any other spin with which the given spin has an exchange interaction.

For notational convenience, we define the reduced spin susceptibilities $`\chi ^{}`$ and $`\overline{\chi ^{}}`$, reduced temperatures $`t`$ and $`\overline{t}`$ and reduced spin gaps $`\mathrm{\Delta }^{}`$ and $`\overline{\mathrm{\Delta }^{}}`$ as

$$\chi ^{}\frac{\chi J^{\mathrm{max}}}{Ng^2\mu _\mathrm{B}^2},\overline{\chi ^{}}\frac{\chi J}{Ng^2\mu _\mathrm{B}^2},$$
(2)
$$t\frac{k_\mathrm{B}T}{J^{\mathrm{max}}},\overline{t}\frac{k_\mathrm{B}T}{J},$$
(3)
$$\mathrm{\Delta }^{}\frac{\mathrm{\Delta }}{J^{\mathrm{max}}},\overline{\mathrm{\Delta }^{}}\frac{\mathrm{\Delta }}{J},$$
(4)
where $`J^{\mathrm{max}}`$ and $`J`$ are, respectively, the largest and average exchange constants in the system, $`N`$ is the number of spins, $`g`$ is the spectroscopic splitting factor appropriate to the direction of the applied magnetic field relative to the crystallographic axes, and $`\mu _\mathrm{B}`$ is the Bohr magneton.

### A High-temperature series expansions for the spin susceptibility and magnetic specific heat

For any Heisenberg spin lattice (in any dimension) in which the spins are magnetically equivalent, i.e. where each spin has identical magnetic coordination spheres, the first three to four terms of the exact quantum mechanical high temperature series expansion of $`\chi ^{}(t)`$ have the same form, with a particularly simple form if the series is inverted. For $`S=1/2`$, one obtains

$$\frac{1}{4\chi ^{}t}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{d_n}{t^n},$$
(6)
$$d_0=1,d_1=\frac{1}{4J^{\mathrm{max}}}\underset{j}{}J_{ij},d_2=\frac{1}{8J_{}^{\mathrm{max}}{}_{}{}^{2}}\underset{j}{}J_{ij}^2,$$
(7)
$$d_3=\frac{1}{24J_{}^{\mathrm{max}}{}_{}{}^{3}}\underset{j}{}J_{ij}^3.$$
(8)
Equation (7) is universal, but Eq. (8) holds only for spin lattices which are not geometrically frustrated for AF ordering and in which the magnetic and crystallographic nearest neighbors of a given spin are the same. Geometrically frustrated lattices typically contain closed triangular exchange paths within the spin lattice structure, such as in the 2D triangular lattice or in the 3D $`B`$ sublattices of the fcc $`AB_2\mathrm{O}_4`$ oxide normal-spinel and $`A_2B_2\mathrm{O}_7`$ oxide pyrochlore structures. The uniform and alternating-exchange chains considered in this paper are not geometrically frustrated, and the magnetic and crystallographic nearest neighbors of a given spin are the same. It has been found that the terms to $`𝒪(1/t^3)`$ on the right-hand-side of Eq. (6) are sufficient to quite accurately describe the susceptibilities of a variety of nonfrustrated zero-, one-, and two-dimensional $`S=1/2`$ AF Heisenberg spin lattices to surprisingly low temperatures $`t1`$. Higher order $`d_n/t^n`$ terms with $`n4`$ are dependent on the structure and dimensionality of the spin lattice. The Weiss temperature $`\theta `$ in the Curie-Weiss law $`\chi (T)=C/(T\theta )`$ is given by the universal expression $`\theta =d_1J^{\mathrm{max}}/k_\mathrm{B}`$.

Because the spin susceptibility and the magnetic contribution $`C(T)`$ to the specific heat can both be expressed, via the fluctuation-dissipation theorem and the Heisenberg Hamiltonian, respectively, in terms of the spin-spin correlation functions, there is a close relationship between these two quantities. In particular, just as there is a universal expression for the first three to four HTSE terms for $`\chi (T)`$ of a Heisenberg spin lattice as discussed above, a universal expression for the first one to two HTSE terms for $`C(T)`$ of such a spin lattice exists and is given for $`S=1/2`$ by

$$\frac{C(t)}{Nk_\mathrm{B}}=\frac{3}{32}\left[\frac{_jJ_{ij}^2}{t^2J_{}^{\mathrm{max}}{}_{}{}^{2}}+\frac{_jJ_{ij}^3}{2t^3J_{}^{\mathrm{max}}{}_{}{}^{3}}+𝒪\left(\frac{1}{t^4}\right)\right].$$
(9)
The sums are over all magnetic nearest-neighbor bonds of any given spin $`𝑺_i`$. The first term is universal but the second term holds only for geometrically nonfrustrated spin lattices in which the crystallographic and magnetic nearest-neighbors of any given spin are the same. Higher order terms all depend on the structure and dimensionality of the spin lattice.

A common misconception is that $`C=0`$ if the magnetic susceptibility of a local-moment system obeys the Curie-Weiss law. This is only true classically. For Heisenberg spin lattices, one can easily show that the Weiss temperature $`\theta `$ in the Curie-Weiss law arises from the first HTSE term \[$`𝒪(1/t)`$\] of the magnetic nearest-neighbor spin-spin correlation function, which is the same quantity that the first HTSE term of $`C(t)`$ arises from. Thus, e.g., for $`S=1/2`$ Heisenberg spin lattices at temperatures $`t1`$ at which the Curie-Weiss law holds, the magnetic specific heat is given by the universal first term of Eq. (9).

### B Low-temperature limit of the spin susceptibility and specific heat of 1D systems with a spin gap

##### Magnetic susceptibility.

For one-dimensional (1D) $`S=1/2`$ Heisenberg spin systems with a spin gap such as the $`S=1/2`$ two-leg ladder (and the alternating-exchange chain), Troyer, Tsunetsugu, and Würtz derived a general expression for $`\chi ^{}(t)`$ which approximately takes into account kinematic magnon interactions, given by

$$\chi ^{}(t)=\frac{1}{t}\frac{z(t)}{1+3z(t)},$$
(11)
$$z(t)=\frac{1}{\pi }_0^\pi \mathrm{e}^{\epsilon _k/t}d(ka),$$
(12)
where $`\epsilon _kE(k)/J^{\mathrm{max}}`$, $`E(k)`$ is the nondegenerate one-magnon (triplet) dispersion relation (the Zeeman degeneracy is already accounted for) and $`a`$ is the (average) distance between spins. This expression is exact in both the low- and high-temperature limits. For the isolated dimer, for which $`\epsilon _k=\mathrm{\Delta }^{}=1`$, Eq. (11) is exact at all temperatures. Inserting $`z(t)=\mathrm{e}^{1/t}`$ for the dimer into Eq. (11) yields the correct result

$$\chi ^{}(t)=\frac{\mathrm{e}^{1/t}}{t}\frac{1}{1+3\mathrm{e}^{1/t}},(\mathrm{dimer})$$
(14)
$$\chi ^{}(t0)=\frac{\mathrm{e}^{1/t}}{t}.$$
(15)
The $`\chi ^{}(t)`$ in Eq. (14) for the antiferromagnetic Heisenberg dimer is plotted in Fig. 1; the fit shown in the figure will be presented and discussed later in Sec. IV B.

At low temperatures $`t\mathrm{\Delta }^{}`$ and $`t`$ one-magnon bandwidth/$`J^{\mathrm{max}}`$, and for a dispersion relation with a parabolic dependence on wave vector $`k`$ near the band minimum

$$\epsilon _k\frac{E(k)}{J^{\mathrm{max}}}=\mathrm{\Delta }^{}+c^{}(ka)^2,$$
(16)
one can replace $`\epsilon _k`$ in Eq. (12) by the approximation (16) and replace the upper limit of the integral in Eq. (12) by $`\mathrm{}`$, yielding $`z(t)=\mathrm{e}^{\mathrm{\Delta }^{}/t}\sqrt{t}/(2\sqrt{\pi c^{}})`$. Substituting this result into Eq. (11) gives the low-$`t`$ limit

$$\chi ^{}(t0)=\frac{A}{t^\gamma }\mathrm{e}^{\mathrm{\Delta }^{}/t},$$
(18)
$$A=\frac{1}{2\sqrt{\pi c^{}}},\gamma =\frac{1}{2}.$$
(19)
This result is correct for any 1D $`S=1/2`$ Heisenberg spin system with a spin gap and with a nondegenerate (excluding Zeeman degeneracy) lowest-lying excited triplet magnon band which is parabolic at the band minimum.

On the other hand, the low-temperature limit of $`\chi ^{}(t)`$ for the isolated dimer in Eq. (15) is of the same form as Eq. (18), but with $`\gamma =1`$. Thus, for 1D systems consisting essentially of dimers which are weakly coupled to each other, a crossover from $`\gamma =1`$ to $`\gamma =1/2`$ is expected with decreasing $`t`$.

The parameters $`A`$ and $`\gamma `$ can be determined if very accurate $`\chi ^{}(t)`$ and $`\mathrm{\Delta }^{}`$ data are available. Taking the logarithm of Eq. (18) yields the low-$`t`$ prediction

$$\mathrm{ln}[\chi ^{}(t)]+\frac{\mathrm{\Delta }^{}}{t}=\mathrm{ln}A\gamma \mathrm{ln}t,$$
(21)
so plotting the left-hand-side vs $`\mathrm{ln}t`$ allows these two parameters to be determined. Alternatively, assuming $`\gamma =1/2`$, one can obtain estimates of $`A`$ and $`\mathrm{\Delta }^{}`$ using Eq. (18), according to

$$\mathrm{ln}(\chi ^{}\sqrt{t})=\mathrm{ln}A+\frac{\mathrm{\Delta }^{}}{t}$$
(22)
and/or

$$\frac{\mathrm{ln}(\chi ^{}\sqrt{t})}{(1/t)}=\mathrm{\Delta }^{}.$$
(23)
##### Specific heat.

The low-$`t`$ limit of the magnetic contribution $`C(T)`$ to the specific heat for the same model is calculated to be

$`{\displaystyle \frac{C(t0)}{Nk_\mathrm{B}}}={\displaystyle \frac{3}{2}}`$ $`\left({\displaystyle \frac{\mathrm{\Delta }^{}}{\pi c^{}}}\right)^{1/2}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{t}}\right)^{3/2}`$ (25)
$`\times \left[1+{\displaystyle \frac{t}{\mathrm{\Delta }^{}}}+{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{t}{\mathrm{\Delta }^{}}}\right)^2\right]\mathrm{e}^{\mathrm{\Delta }^{}/t}.`$
Note that, in addition to the ratio $`t/\mathrm{\Delta }^{}=k_\mathrm{B}T/\mathrm{\Delta }`$ of the thermal energy to the spin gap, the magnitude of $`\chi ^{}`$ in Eqs. (16) is determined by the actual value of the curvature $`c^{}`$ at the triplet one-magnon band minimum, whereas the magnitude of $`C`$ in Eq. (25) depends only on the ratio of $`c^{}`$ to $`\mathrm{\Delta }^{}`$. These formulas have been applied in the literature to model experimental data for alternating-exchange chain and two-leg spin ladder compounds. However, with one exception to our knowledge, these modeling studies have not recognized that the prefactor parameter and the spin gap are not independently adjustable parameters. For a given spin lattice, they are in fact uniquely related to each other. Their relationship for the $`S=1/2`$ two-leg Heisenberg ladder was studied in Ref. . For the alternating-exchange chain, we estimate the relationship between $`c^{}`$ and $`\mathrm{\Delta }^{}`$ below in Sec. II C 3.

### C Alternating-exchange chain

The $`S=1/2`$ AF alternating-exchange Heisenberg chain Hamiltonian is written in three equivalent ways as

$``$ $`=`$ $`{\displaystyle \underset{i}{}}J_1\stackrel{}{S}_{2i1}\stackrel{}{S}_{2i}+J_2\stackrel{}{S}_{2i}\stackrel{}{S}_{2i+1}`$ (27)
$`=`$ $`{\displaystyle \underset{i}{}}J_1\stackrel{}{S}_{2i1}\stackrel{}{S}_{2i}+\alpha J_1\stackrel{}{S}_{2i}\stackrel{}{S}_{2i+1}`$ (28)
$`=`$ $`{\displaystyle \underset{i}{}}J(1+\delta )\stackrel{}{S}_{2i1}\stackrel{}{S}_{2i}+J(1\delta )\stackrel{}{S}_{2i}\stackrel{}{S}_{2i+1},`$ (29)
where

$`J_1`$ $`=`$ $`J(1+\delta )={\displaystyle \frac{2J}{1+\alpha }},`$ (31)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{J_2}{J_1}}={\displaystyle \frac{1\delta }{1+\delta }},`$ (33)
$`\delta `$ $`=`$ $`{\displaystyle \frac{J_1}{J}}1={\displaystyle \frac{J_1J_2}{2J}}={\displaystyle \frac{1\alpha }{1+\alpha }},`$ (35)
$`J`$ $`=`$ $`{\displaystyle \frac{J_1+J_2}{2}}=J_1{\displaystyle \frac{1+\alpha }{2}},`$ (37)
with AF couplings $`J_1J_20`$, $`0(\alpha ,\delta )1`$. The uniform undimerized chain corresponds to $`\alpha =1,\delta =0,J_1=J_2=J`$. The form of the Hamiltonian in Eq. (29) is most appropriate for chains showing a second-order dimerization transition at $`T_\mathrm{c}`$ with decreasing $`T`$. If the exchange modulation $`\delta 1`$ ($`\alpha 1`$), the (average) $`J`$ below $`T_\mathrm{c}`$ can be identified with the exchange coupling in the high-$`T`$ undimerized state.

In spin-Peierls systems, the spin-phonon interaction causes a lattice dimerization to occur below the spin-Peierls transition temperature, together with a spin-gap due to the formation of spin singlets on the dimers. The Hamiltonian can be mapped onto the spin Hamiltonian (II C) (with renormalized exchange constants) only in the adiabatic parameter regime, in which the relevant phonon energy is much smaller than $`J`$. If such a mapping cannot be made, dynamical phonon effects (quantum mechanical fluctuations) become important and the $`\chi (T)`$ can be significantly different from that predicted from Hamiltonian (II C). This issue will be discussed further when modeling the $`\chi (T)`$ data for NaV2O5 in Sec. VIII B.

#### 1 High-temperature series expansions

##### Magnetic Susceptibility.

For the alternating-exchange chain, according to our definition one has $`J^{\mathrm{max}}=J_1`$. Then using the definition for $`\alpha `$ in Eq. (33), the $`d_n`$ HTSE coefficients in Eqs. (7) and (8) become

$$d_0=1,d_1=\frac{1+\alpha }{4},d_2=\frac{1+\alpha ^2}{8},d_3=\frac{1+\alpha ^3}{24}.$$
(38)
One can change variables from $`\alpha `$ and $`J_1`$ in $`\chi ^{}(\alpha ,t)`$ to $`\delta `$ and $`J`$ in $`\overline{\chi ^{}}(\delta ,\overline{t})`$ using Eqs. (II C) which give

$$t=\frac{\overline{t}}{1+\delta },$$
(40)
$$\overline{\chi ^{}}(\delta ,\overline{t})=\frac{1}{1+\delta }\chi ^{}(\frac{1\delta }{1+\delta },\frac{\overline{t}}{1+\delta }).$$
(41)
We write the resulting HTSE for the inverse of $`\overline{\chi ^{}}(\delta ,\overline{t})`$ as

$$\frac{1}{4\overline{\chi ^{}}\overline{t}}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{\overline{d}_n}{\overline{t}^n},$$
(43)
where we find

$$\overline{d}_0=1,\overline{d}_1=\frac{1}{2},\overline{d}_2=\frac{1+\delta ^2}{4},\overline{d}_3=\frac{1+3\delta ^2}{12}.$$
(44)
An important feature of this HTSE of $`\overline{\chi ^{}}(\delta ,\overline{t})`$ is that it is an even (analytic) function of $`\delta `$ for any finite temperature. This constraint must be true in general and not just for the terms listed, because $`\overline{\chi ^{}}(\delta ,\overline{t})`$ cannot depend on the sign of $`\delta =(J_1J_2)/(2J)`$: the Hamiltonian in Eq. (29) is invariant upon such a sign change. Physically, a negative $`\delta `$ would simply correspond to relabeling all $`𝑺_i𝑺_{i+1}`$, which cannot change the physical properties. We will use this constraint that $`\overline{\chi ^{}}(\delta ,\overline{t})`$ must be an even function of $`\delta `$ to help formulate our fitting function (after a change back in variables) for our QMC and TMRG $`\chi ^{}(\alpha ,t)`$ calculations for the alternating-exchange chain. This constraint is important because it allows a fit function for $`\chi ^{}(\alpha ,t)`$ to be formulated which is accurate for $`\alpha 1`$ ($`\delta 1`$), a parameter regime relevant to compounds exhibiting second-order spin-dimerization transitions with decreasing temperature.

##### Magnetic specific heat.

Using $`J^{\mathrm{max}}`$ = $`J_1`$ and $`\alpha =J_2/J_1`$, the general HTSE expression in Eq. (9) yields the two lowest-order HTSE terms for the magnetic specific heat $`C(T)`$ of the $`S=1/2`$ AF alternating-exchange Heisenberg chain as

$$\frac{C(t)}{Nk_\mathrm{B}}=\frac{3}{32}\left[\frac{1+\alpha ^2}{t^2}+\frac{1+\alpha ^3}{2t^3}+𝒪\left(\frac{1}{t^4}\right)\right].$$
(45)
#### 2 Spin gap

The spin gap $`\mathrm{\Delta }^{}(\alpha )`$ of the alternating-exchange chain was determined to high ($`1\%`$) accuracy for $`0\alpha 0.9`$, in $`\alpha `$ increments of 0.1, using multiprecision methods by Barnes, Riera, and Tennant (BRT). They found that their calculations could be parametrized well by

$$\mathrm{\Delta }^{}(\alpha )\frac{\mathrm{\Delta }(\alpha )}{J_1}(1\alpha )^{3/4}(1+\alpha )^{1/4},$$
(47)
$$\overline{\mathrm{\Delta }^{}}(\delta )\frac{\mathrm{\Delta }(\delta )}{J}2\delta ^{3/4}.$$
(48)
The same $`\overline{\mathrm{\Delta }^{}}(\delta )`$ was found in numerical calculations by Ladavac et al. for $`0.01\delta 1`$, whereas calculations for $`0.03\delta 0.06`$ by Augier et al. yielded somewhat smaller values of $`\overline{\mathrm{\Delta }^{}}`$ than predicted by Eq. (48). The asymptotic critical behavior of $`\overline{\mathrm{\Delta }^{}}`$ as the uniform limit is approached ($`\alpha 1,\delta 0`$) has been given as

$$\overline{\mathrm{\Delta }^{}}(\delta )\frac{\delta ^{2/3}}{|\mathrm{ln}\delta |^{1/2}};$$
(49)
thus the parametrization in Eq. (48) evidently indicates that the fitted data do not reside within the asymptotic critical regime. Alternatively, Barnes, Riera, and Tennant suggested that Eq. (49) may not be the correct form for the asymptotic critical behavior. On the other hand, Uhrig et al. fitted their $`T=0`$ density matrix renormalization group (DMRG) calculations of $`\overline{\mathrm{\Delta }^{}}(\delta )`$ for $`0.004\delta 0.1`$ to a power-law behavior without the log correction and obtained $`\overline{\mathrm{\Delta }^{}}1.57\delta ^{0.65}`$. We will further discuss the above spin gap calculation results later in Sec. VII A after deriving our own $`\mathrm{\Delta }^{}(\alpha )`$ values from our TMRG $`\chi ^{}(\alpha ,t)`$ data in Sec. VI.

#### 3 One-magnon dispersion relations

Barnes, Riera, and Tennant have computed the dimer series expansion of the dispersion relation $`\epsilon (\alpha ,k)E(\alpha ,k)/J_1`$ for the one-magnon ($`S=1`$) energy $`E(\alpha ,k)`$ vs wave vector $`k`$ along the chain for the lowest-lying one-magnon band, up to fifth order in $`\alpha `$, which we write as

$$\epsilon (\alpha ,k)=\underset{n=0}{\overset{\mathrm{}}{}}a_n(\alpha )\mathrm{cos}(2nka),$$
(50)
where $`a`$ is the (average) spin-spin distance, which is 1/2 the lattice repeat distance along the chain in the dimerized state. Plots of $`\epsilon (\alpha ,k)`$ for $`\alpha =0`$, 0.2, 0.4, 0.6, and 0.8 up to fifth order in $`\alpha `$, as given in Fig. 4 of Ref. , are shown as the dashed curves in Fig. 2. The curves are symmetric about $`ka=\pi /2`$, so the same spin gap $`\mathrm{\Delta }^{}(\alpha )=_{n=0}^{\mathrm{}}a_n(\alpha )`$ occurs at $`ka=0`$ and $`\pi `$. This fifth-order approximation yields $`\mathrm{\Delta }^{}(\alpha )`$ values for $`\alpha 0.9`$ in rather close agreement with BRTs’ results discussed in the previous section. For a dimer series expansion we expect the average energy of the one-magnon band states to be nearly independent of $`\alpha `$, i.e.,

$$\frac{1}{\pi }_0^\pi \epsilon (\alpha ,ka)d(ka)=1.$$
(51)
Indeed, upon inserting BRTs’ fifth order expansion coefficients into Eq. (50) and the result into Eq. (51), we find that this sum rule is satisfied to within 1% for $`0\alpha 1`$.

Also shown as a solid curve in Fig. 2 is the exact result $`\epsilon (k)=(\pi /2)|\mathrm{sin}(ka)|`$ for the uniform chain ($`\alpha =1`$). This $`\epsilon (k)`$ has a cusp with infinite curvature (at $`ka=0`$

and $`\pi `$) which cannot be accurately approximated by a Fourier series with a small number of terms. This singular behavior is evidently closely related to the critical behavior of $`\mathrm{\Delta }^{}(\alpha 1)`$ discussed above. In order to later model our TMRG $`\chi ^{}(\alpha ,t)`$ data close to, but not in, the low-$`t`$ limit, we will need an expression for $`\epsilon (\alpha ,k)`$ which is correct in the limit $`\alpha 1`$ and which also reproduces reasonably well the $`\epsilon (\alpha ,k)`$ of BRT. We found that the simple one-parameter ($`\mathrm{\Delta }^{}`$) form suggested earlier by one of us in the context of the $`S=1/2`$ two-leg ladder

$$\epsilon (\mathrm{\Delta }^{},k)=[\mathrm{\Delta }_{}^{}{}_{}{}^{2}+f^2(\mathrm{\Delta }^{})\mathrm{sin}^2(ka)]^{1/2},$$
(53)
is satisfactory in these regards for the AF alternating-exchange chain over the entire range $`0\mathrm{\Delta }^{}1`$. The function $`f(\mathrm{\Delta }^{})`$ is determined here by the sum rule (51), which yields the condition

$$\mathrm{E}\left[\frac{f^2(\mathrm{\Delta }^{})}{\mathrm{\Delta }_{}^{}{}_{}{}^{2}}\right]=\frac{\pi }{2\mathrm{\Delta }^{}},$$
(54)
where E($`x`$) is the complete elliptic integral of the second kind. From Eq. (54), $`f`$ varies nonlinearly with $`\mathrm{\Delta }^{}`$ from $`f(\mathrm{\Delta }^{}=0)=\pi /2`$ to $`f(\mathrm{\Delta }^{}=1)=0`$, as shown in Fig. 3. From an independently determined dependence of $`\mathrm{\Delta }^{}`$ on $`\alpha `$ as in Eq. (47), one can then determine $`f(\alpha )`$ as also shown in Fig. 3. Using the fifth-order $`\mathrm{\Delta }^{}(\alpha )`$ values of BRT in Fig. 2, the resulting dispersion relations (3) for $`\alpha =0,`$ 0.2, 0.4, 0.6, and 0.8 were calculated and are shown as the solid curves in Fig. 2, where they are seen to be in close agreement with the respective dashed curves of BRT. An important difference for large $`\alpha `$, however, is that the $`\epsilon (\mathrm{\Delta }^{},k)`$ in Eqs. (3) properly reduces by construction to the exact $`\epsilon (\alpha ,k)`$ for $`\alpha 1`$, whereas the one in Eq. (50) with a finite number of terms does not.

Close to the one-magnon band minimum, the square root and the sine function in the dispersion relation in Eq. (53) can be expanded, yielding

$$\epsilon (\mathrm{\Delta }^{},k0)\mathrm{\Delta }^{}+\frac{1}{2}\frac{f^2(\mathrm{\Delta }^{})}{\mathrm{\Delta }^{}}(ka)^2.$$
(55)
A comparison of Eqs. (55) and (16) shows that the parameter $`c^{}`$ in the formulas for $`\chi ^{}(t0)`$ \[Eqs. (16)\] and $`C(t0)`$ \[Eq. (25)\] is a unique function of $`\mathrm{\Delta }^{}`$ which in our approximation is given by

$$c^{}(\mathrm{\Delta }^{})=\frac{1}{2}\frac{f^2(\mathrm{\Delta }^{})}{\mathrm{\Delta }^{}},$$
(56)
with $`f^2(\mathrm{\Delta }^{})`$ given by Eq. (54). Thus both $`\chi ^{}(t0)`$ and $`C(t0)`$ for the alternating-exchange chain only depend on the single parameter $`\mathrm{\Delta }^{}`$ (in addition to $`t`$). Explicitly, we obtain

$$\chi ^{}(t0)=\frac{1}{\sqrt{2\pi }f(\mathrm{\Delta }^{})}\left(\frac{\mathrm{\Delta }^{}}{t}\right)^{1/2}\mathrm{e}^{\mathrm{\Delta }^{}/t}.$$
(57)
As might have been anticipated, the only thermodynamic variable is the ratio $`t/\mathrm{\Delta }^{}=k_\mathrm{B}T/\mathrm{\Delta }`$ of the thermal energy to the spin gap. The numerical prefactor depends explicitly (only) on the reduced spin gap $`\mathrm{\Delta }^{}\mathrm{\Delta }/J_1`$. Similarly, the magnetic specific heat is obtained as

$`{\displaystyle \frac{C(t0)}{Nk_\mathrm{B}}}={\displaystyle \frac{3}{\sqrt{2\pi }}}`$ $`{\displaystyle \frac{\mathrm{\Delta }^{}}{f(\mathrm{\Delta }^{})}}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{t}}\right)^{3/2}`$ (59)
$`\times \left[1+{\displaystyle \frac{t}{\mathrm{\Delta }^{}}}+{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{t}{\mathrm{\Delta }^{}}}\right)^2\right]\mathrm{e}^{\mathrm{\Delta }^{}/t},`$
where again the same characteristics are present as just discussed for $`\chi ^{}(t)`$. The variations of the prefactors with $`\mathrm{\Delta }^{}`$ for $`\chi ^{}(t)`$ and $`C(t)`$ can both be ascertained from the plot of $`f(\mathrm{\Delta }^{})`$ in Fig. 3. In particular, when $`\alpha 1`$ ($`\delta 1`$) for which $`\mathrm{\Delta }^{}1`$, $`f`$ is nearly a constant. For our and our readers’ convenience when modeling materials showing small spin gaps, we have fitted our numerical $`f(\mathrm{\Delta }^{})`$ calculations for $`0\mathrm{\Delta }^{}0.4`$ by a third-order polynomial to within 2 parts in $`10^4`$, given by

$$f(\mathrm{\Delta }^{})=\frac{\pi }{2}0.034289\mathrm{\Delta }^{}1.18953\mathrm{\Delta }_{}^{}{}_{}{}^{2}+0.40030\mathrm{\Delta }_{}^{}{}_{}{}^{3}.$$
(60)
By a change in variables to $`(J,\delta )`$ and using the $`\mathrm{\Delta }^{}(\alpha )`$ in Eq. (II C 2), we obtain the following forms which are more useful for modeling materials with small spin gaps, especially those showing second-order spin dimerization transitions with decreasing $`T`$:

$$\overline{\chi ^{}}(\overline{t}0)=\frac{1}{(1+\delta )\sqrt{2\pi }f(\overline{\mathrm{\Delta }^{}})}\left(\frac{\overline{\mathrm{\Delta }^{}}}{\overline{t}}\right)^{1/2}\mathrm{e}^{\overline{\mathrm{\Delta }^{}}/\overline{t}}.$$
(62)
$`{\displaystyle \frac{\overline{C}(t0)}{Nk_\mathrm{B}}}=`$ $`{\displaystyle \frac{3}{(1+\delta )\sqrt{2\pi }}}{\displaystyle \frac{\overline{\mathrm{\Delta }^{}}}{\overline{f}(\overline{\mathrm{\Delta }^{}})}}\left({\displaystyle \frac{\overline{\mathrm{\Delta }^{}}}{\overline{t}}}\right)^{3/2}`$ (64)
$`\times \left[1+{\displaystyle \frac{\overline{t}}{\overline{\mathrm{\Delta }^{}}}}+{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{\overline{t}}{\overline{\mathrm{\Delta }^{}}}}\right)^2\right]\mathrm{e}^{\overline{\mathrm{\Delta }^{}}/\overline{t}},`$
$$\overline{f}(\overline{\mathrm{\Delta }^{}})=\frac{\pi }{2}0.033933\overline{\mathrm{\Delta }^{}}1.19607\overline{\mathrm{\Delta }^{}}^2+0.92430\overline{\mathrm{\Delta }^{}}^3.$$
(65)
Note that in these formulas $`\mathrm{\Delta }^{}/t=\overline{\mathrm{\Delta }^{}}/\overline{t}=\mathrm{\Delta }/(k_\mathrm{B}T)`$.

## III Theory: $`𝑺\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$ Uniform Heisenberg Chain

### A Magnetic spin susceptibility

The uniform $`S=1/2`$ chain is one limit of the alternating-exchange chain with $`J_{ij}J,\alpha =1,\delta =0`$, and with no spin gap \[the $`\chi ^{}(t0)`$ and $`C(t0)/t`$ are finite\]. The spin susceptibility was calculated accurately by Eggert, Affleck and Takahashi in 1994, and recently refined by Klümper as shown in Fig. 4(a) where only the calculations up to $`t=2`$ are shown. An expanded plot of the data for $`t0.02`$, including the exact value $`1/\pi ^2`$ at $`t=0`$, is shown in Fig. 4(b), along with a fit (Fit 2) to the data to be derived and discussed in Sec. IV C. The most recent calculations of Ref.  have an absolute accuracy estimated to be $`1\times 10^9`$ and show a broad maximum at a temperature $`T^{\mathrm{max}}`$ with a value $`\chi ^{\mathrm{max}}`$. By fitting data points near the maximum by up to 6th order polynomials, we determined these numerical values to be given by

$$T^{\mathrm{max}}=\mathrm{0.6\hspace{0.17em}408\hspace{0.17em}510}(4)J/k_\mathrm{B},$$
(67)
$$\frac{\chi ^{\mathrm{max}}J}{Ng^2\mu _\mathrm{B}^2}=\mathrm{0.146\hspace{0.17em}926\hspace{0.17em}279}(1),$$
(68)
$$\chi ^{\mathrm{max}}T^{\mathrm{max}}=\mathrm{0.0\hspace{0.17em}941\hspace{0.17em}579}(1)\frac{Ng^2\mu _\mathrm{B}^2}{k_\mathrm{B}}.$$
(69)
These values are consistent within the errors with those found by Eggert et al., but are much more accurate. For one mole of spins, setting $`N=N_\mathrm{A}`$ (Avogadro’s number) in Eq. (69) yields

$$\chi ^{\mathrm{max}}T^{\mathrm{max}}=\mathrm{0.0\hspace{0.17em}353\hspace{0.17em}229}(3)g^2\frac{\mathrm{cm}^3\mathrm{K}}{\mathrm{mol}}.$$
(70)
Note that the product $`\chi ^{\mathrm{max}}T^{\mathrm{max}}`$ in Eqs. (69) and (70) is independent of $`J`$, and hence is a good initial test of whether the $`S=1/2`$ AF uniform Heisenberg chain model might be applicable to a particular compound.

#### 1 High-temperature series expansions

The coefficients $`c_n`$ of the HTSE for $`\chi ^{}(t)`$,

$$4\chi ^{}t=\underset{n=0}{\overset{\mathrm{}}{}}\frac{c_n}{t^n},$$
(72)
are given up to $`𝒪(1/t^7)`$ by

$`c_0`$ $`=`$ $`1,c_1={\displaystyle \frac{1}{2}},c_2=0,c_3={\displaystyle \frac{1}{24}},c_4={\displaystyle \frac{5}{384}},`$ (73)
$`c_5`$ $`=`$ $`{\displaystyle \frac{7}{1280}},c_6={\displaystyle \frac{133}{30720}},c_7={\displaystyle \frac{1}{4032}}.`$ (74)
Inverting the series, we obtain the corresponding $`d_n`$ coefficients in Eq. (6) as

$`d_0`$ $`=`$ $`1,d_1={\displaystyle \frac{1}{2}},d_2={\displaystyle \frac{1}{4}},d_3={\displaystyle \frac{1}{12}},d_4={\displaystyle \frac{1}{128}},`$ (75)
$`d_5`$ $`=`$ $`{\displaystyle \frac{29}{3840}},d_6={\displaystyle \frac{317}{92160}},d_7={\displaystyle \frac{11}{71680}}.`$ (76)
The $`d_n`$ coefficients with $`n=0,`$ 1, 2, and 3 are of course in agreement with Eq. (38) for $`\alpha =1`$.

#### 2 Logarithmic corrections at low temperatures

At low temperatures a simple expansion of thermodynamic properties in for instance the variable $`t`$ is not possible. Such a nonanalyticity in $`t`$ can be viewed as due to the strong correlation of the quasiparticles, i.e., the elementary excitations of the system are not strictly free; they show rather nontrivial scattering processes. Spinons with low energies $`ϵ_1`$ and $`ϵ_2`$ have a scattering phase $`\varphi (ϵ_1,ϵ_2)\varphi _0+\text{const}/|\mathrm{log}(ϵ_1ϵ_2)|`$. From this property it is clear that an expansion in the single variable $`t`$ is not possible, but has to be supplemented by a term $`1/\mathrm{log}(t)`$. Although being very intuitive, this physical picture on the basis of scattering processes of spinons has not played any important role in the investigation of logarithmic corrections until recently.

Logarithmic dependencies of physical quantities have been known for the spin-1/2 Heisenberg chain for a rather long time. Usually, a quantum chain with critical couplings leads to critical correlations only in the thermodynamic limit $`1/L=0`$ and at $`T=0`$, where $`L`$ is the length of the chain. If one of these conditions is not met the physical properties receive nonanalytic contributions in terms of $`1/L`$ or $`T`$. From the renormalization group point of view the existence of logarithmic corrections is reflected by the perturbation of the (critical) fixed point Hamiltonian by some marginal operator. Such operators usually exist only for isotropic systems.

The investigation of the size dependence of energy levels of critical quantum chains was started more than a decade ago. For the isotropic spin-1/2 Heisenberg chain, expansions in $`1/L`$ and additional logarithmic corrections ($`1/L\mathrm{log}L`$ etc.) were found in lattice approaches (Bethe ansatz) as well as in field theory \[RG study of the Wess-Zumino-Witten (WZW) model with topological term $`k=1`$ (Refs. )\].

Many of these earlier results are still relevant for the issues discussed in this section due to an equivalence of many-particle systems at $`T=0`$, $`1/L>0`$ (groundstate properties of finite chains) and those at $`T>0`$, $`1/L=0`$ (thermodynamics of the bulk). This leads to asymptotic series where $`T`$ and $`1/L`$ play very similar roles. To our knowledge the first explicit report on $`\mathrm{log}(T)`$ corrections in the magnetic susceptibility resulting in an infinite slope at $`T=0`$ was given in Ref. . Including higher order terms, the asymptotic expansion $`\chi _{\mathrm{lt}}^{}(t)`$ for $`\chi ^{}(t)`$ is

$`\chi _{\mathrm{lt}}^{}(t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi ^2}}\left[1+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{ln}\left(+\frac{1}{2}\right)}{(2)^2}}+\mathrm{}\right],`$ (78)
$`=`$ $`{\displaystyle \frac{1}{\pi ^2}}\left[1+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{ln}}{(2)^2}}{\displaystyle \frac{1}{(2)^3}}+\mathrm{}\right],`$ (79)
$$\mathrm{ln}(t_0/t),$$
(80)
where $`t_0`$ is a nonuniversal (undetermined) parameter. In Ref.  the field theoretical prediction on the basis of the WZW model was compared with the results of thermodynamic Bethe ansatz calculations and showed convincing agreement in an intermediate temperature regime. Using up to the first logarithmic correction term in Eq. (78), Eggert, Affleck, and Takahashi estimated $`t_07.7`$, so at low temperatures $`t0.01`$ the parameter $`1`$.

A general feature of field theoretical and lattice approaches is their restriction to “low” and “high” temperatures, respectively. Field theoretical studies suffer at high temperatures from the neglect of (infinitely many) irrelevant operators. Lattice studies show convergence problems at low temperatures as increasingly larger systems have to be studied in order to avoid finite-size effects. In addition, the comparison of field theory and lattice results can only verify or falsify the universal aspects of an asymptotic expansion. Non-universal quantities like $`t_0`$ which derive from some coupling constant of a marginal or irrelevant operator (undetermined within the field theory) can at best be fitted as done in Ref. .

The latter problem of determining the coupling constants in an effective field theory was solved by Lukyanov who used a bosonic representation of the Heisenberg chain. The coupling constants were fixed by a comparison of the susceptibility data $`\chi (T=0,h)`$ obtained by him with Bethe ansatz calculations for magnetic field $`h`$ at $`T=0`$. Eventually, the $`\chi (T>0,h=0)`$ data could be calculated without any need of a fit parameter.

Lukyanov obtained the following analytical low-temperature expansion of $`\chi ^{}(t)`$,

$`\chi _{\mathrm{lt},\mathrm{g}}^{}(t)={\displaystyle \frac{1}{\pi ^2}}\{1`$ $`+`$ $`{\displaystyle \frac{g}{2}}+{\displaystyle \frac{3g^3}{32}}+𝒪(g^4)`$ (82)
$`+`$ $`{\displaystyle \frac{\sqrt{3}}{\pi }}t^2[1+𝒪(g)]\},`$ (83)
where $`g`$ obeys the transcendental equation

$$\frac{1}{2}\mathrm{ln}g+\frac{1}{g}=$$
(84)
or equivalently

$$\sqrt{g}\mathrm{e}^{1/g}=\frac{t_0}{t},$$
(85)
with a unique value of $`t_0`$ given by

$$t_0=\sqrt{\frac{\pi }{2}}\mathrm{e}^{\gamma +(1/4)}\mathrm{2.866\hspace{0.17em}257\hspace{0.17em}058},$$
(86)
where $`\gamma \mathrm{0.577\hspace{0.17em}215\hspace{0.17em}665}`$ is Euler’s constant. Lukyanov showed that his $`\chi _{\mathrm{lt},\mathrm{g}}^{}(t)`$ is in agreement with the numerical data of Eggert, Affleck and Takahashi at low temperatures $`t0.003`$.

In the following, we will compare high-accuracy numerical Bethe ansatz calculations carried out to much lower temperatures by Klümper and Johnston with this theory in some detail because this theory is exact at low temperatures with no adjustable parameters. The calculations of Ref.  are based on lattice studies, however without suffering from the usual shortcomings. By means of a lattice path integral representation of the finite temperature Heisenberg chain and the formulation of a suitable quantum transfer matrix (both quite analogous to the numerical TMRG calculations presented later in this paper) a set of numerically well-posed expressions for the free energy was derived. In more physical terms the method can be understood as an application of the familiar though often rather vague concept of quasiparticles to a quantitative description of the many particle system valid for all temperatures $`T`$ and magnetic field values $`h`$. The work can be understood as an evaluation of the full scattering theory of spinons and antispinons.

Our iterative solution of Eq. (84) yields the expansion

$$g=\frac{1}{}\left\{1\frac{\mathrm{ln}}{2}+\frac{(\mathrm{ln})^2\mathrm{ln}}{(2)^2}+𝒪\left[\frac{1}{(2)^3}\right]\right\}.$$
(87)
A log-log plot of $`g`$ vs $`t`$ obtained by numerically solving Eq. (85) is shown in Fig. 5 (solid curve), along with its lowest-order approximation $`g(t)1/=1/\mathrm{ln}(t_0/t)`$ (dashed curve). This approximation is 1.1% larger than the exact result at $`t=10^{30}`$, with the discrepancy increasing steadily to 5.8% at $`t=10^{15}`$ and 8.5% at $`t=10^7`$. Substituting Eq. (87) into (83) gives

$`\chi _{\mathrm{lt},\mathrm{log}}^{}(t)={\displaystyle \frac{1}{\pi ^2}}\{1`$ $`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{ln}}{(2)^2}}`$ (88)
$`+`$ $`{\displaystyle \frac{(\mathrm{ln})^2\mathrm{ln}+(3/4)}{(2)^3}}+𝒪\left[{\displaystyle \frac{1}{(2)^4}}\right]`$ (89)
$`+`$ $`{\displaystyle \frac{\sqrt{3}}{\pi }}t^2[1+𝒪\left({\displaystyle \frac{1}{2}}\right)]\}.`$ (90)
The first three terms are identical with those in Eq. (79), but the constant term in the numerator of the fourth term is not the same as in Eq. (79), indicating that Eq. (79) is not accurate to order higher than $`𝒪[1/(2)^2]`$.

An important issue is the accuracy to which the log expansion approximation $`\chi _{\mathrm{lt},\mathrm{log}}^{}(t)`$ in Eq. (90) reproduces the $`\chi _{\mathrm{lt},\mathrm{g}}^{}(t)`$ prediction of the original Eqs. (III A 2). We have calculated both quantities to high accuracy and plot the difference vs $`t`$, for the range $`10^{30}t0.5`$, in Fig. 6. The $`\chi _{\mathrm{lt},\mathrm{log}}^{}(t)`$ is seen to increasingly diverge from $`\chi _{\mathrm{lt},\mathrm{g}}^{}(t)`$ with increasing $`t`$.

When comparing the predictions of Lukyanov’s theory with numerical results such as obtained from the Bethe ansatz, it is important to know at what temperature the low temperature expansion in Eqs. (III A 2) ceases to be accurate (“accurate” must be defined) with increasing temperature. There are three aspects of this issue that need to be addressed. The first and second aspects concern the temperatures at which the unknown $`𝒪(g^4)`$ and $`𝒪(g)`$ terms in Eq. (83) become significant, respectively; we will return to these two issues shortly. The third aspect is whether the log expansion approximation $`\chi _{\mathrm{lt},\mathrm{log}}^{}(t)`$ in Eq. (90) can be used in this comparison. The absolute accuracy of the most recent Bethe ansatz calculations is estimated to be $`1\times 10^9`$. From Fig. 6, we see that $`\chi _{\mathrm{lt},\mathrm{log}}^{}(t)`$ approximates $`\chi _{\mathrm{lt},\mathrm{g}}^{}(t)`$ to this degree only for temperatures $`t10^{30}`$ \[we infer that the previous Eqs. (III A 2) only apply to this accuracy at similarly very low temperatures\]. Therefore, to avoid this unnecessary approximation as a source of error at higher temperatures, we will henceforth compare the numerical Bethe ansatz calculations with $`\chi _{\mathrm{lt},\mathrm{g}}^{}(t)`$ calculated from Lukyanov’s original Eqs. (III A 2).

A comparison of the low-temperature Bethe ansatz $`\chi ^{}(t)`$ calculations and Lukyanov’s theoretical $`\chi _{\mathrm{lt},\mathrm{g}}^{}`$ prediction is shown in Fig. 7(a). On the scale of this figure, the two results are identical. The (small) quantitative differences between them are shown as the filled circles in Fig. 7(b). The lower error bar on each data point in Fig. 7(b) is $`1\times 10^7`$ to indicate the scale. The upper error bar is the estimated uncertainty in $`\chi _{\mathrm{lt},\mathrm{g}}^{}`$ arising from the presence of the unknown $`𝒪(g^4)`$ and higher-order terms in Eq. (83), which was set to $`g^4(t)/\pi ^2`$; the uncertainty in the $`t^2`$ contribution, $`\sqrt{3}t^2g(t)/\pi ^3`$, is negligible at low $`t`$ compared to this. At the lower temperatures, the data agree extremely well with the prediction of Lukyanov’s theory. At the highest temperatures $`t10^3`$, higher or-

der $`t^n`$ terms also become important, as inferred from our empirical fits (Fits 1 and 2) below to the numerical data.

Irrespective of the uncertainties in the theoretical prediction at high temperatures just discussed, we can safely conclude directly from Fig. 7(b) that the Bethe ansatz $`\chi ^{}(t)`$ data are in agreement with the exact theory of Lukyanov to within an absolute accuracy of $`1\times 10^6`$ (relative accuracy $`10`$ ppm) over a temperature range spanning 18 orders of magnitude from $`t=5\times 10^{25}`$ to $`t=5\times 10^7`$. The agreement is much better than this at the lower temperatures.

### B Magnetic specific heat

The magnetic specific heat $`C`$ of the $`S=1/2`$ AF uniform Heisenberg chain was recently calculated to high accuracy by Klümper and Johnston over the temperature range $`5\times 10^{25}k_\mathrm{B}T/J5`$. The accuracy is estimated to be $`3\times 10^{10}C(t)`$. The results for $`T2J/k_\mathrm{B}`$ are shown in Fig. 8(a). The initial $`T`$ dependence is approximately (see below) linear, and is given exactly in

the $`t=0`$ limit by

$$\frac{C(t0)}{Nk_\mathrm{B}}=\frac{2}{3}t.$$
(91)
The data show a maximum with a value $`C^{\mathrm{max}}`$ at a temperature $`T_C^{\mathrm{max}}`$. By fitting 3–7 data points in the vicinity of the maximum by up to 6th order polynomials, these values were found to be

$`{\displaystyle \frac{k_\mathrm{B}T_C^{\mathrm{max}}}{J}}`$ $`=`$ $`\mathrm{0.48\hspace{0.17em}028\hspace{0.17em}487}(1),`$ (92)
$`{\displaystyle \frac{C^{\mathrm{max}}}{Nk_\mathrm{B}}}`$ $`=`$ $`\mathrm{0.3\hspace{0.17em}497\hspace{0.17em}121\hspace{0.17em}235}(2).`$ (94)
The electronic specific heat coefficient $`C(T)/T`$ is plotted vs temperature in Fig. 8(b). As expected from Eq. (91), the data approach the value $`(2/3)Nk_\mathrm{B}^2/J`$ for $`t0`$. The initial deviation from this constant value is positive and approximately (see below) quadratic in $`t`$. The data exhibit a smooth maximum with a value $`(C/T)^{\mathrm{max}}`$ at a temperature $`T_{\mathrm{C}/\mathrm{T}}^{\mathrm{max}}`$, values which we determined by fitting polynomials to the data in the vicinity of the peak to be

$`{\displaystyle \frac{k_\mathrm{B}T_{C/T}^{\mathrm{max}}}{J}}`$ $`=`$ $`\mathrm{0.30\hspace{0.17em}716\hspace{0.17em}996}(2),`$ (95)
$`{\displaystyle \frac{(C/T)^{\mathrm{max}}J}{Nk_\mathrm{B}^2}}`$ $`=`$ $`\mathrm{0.8\hspace{0.17em}973\hspace{0.17em}651\hspace{0.17em}576}(5).`$ (97)
The magnetic entropy $`S(T)`$ is determined by integrating the $`C(T)/T`$ data in Fig. 8(b) vs $`T`$ and the result, normalized by $`S(T\mathrm{})=Nk_\mathrm{B}\mathrm{ln}2`$, is plotted vs $`T`$ in Fig. 9. This figure allows one to estimate the maximum magnetic entropy that can be associated with a dimerization transition or any other magnetic transition involving $`S=1/2`$ Heisenberg chains which are weakly coupled to each other \[assuming that the (average) $`J`$ does not change at the transition\]. For example, for $`\mathrm{NaV}_2\mathrm{O}_5`$ where $`k_\mathrm{B}T_\mathrm{c}/J0.057`$, one can estimate from Fig. 9 that the magnetic entropy at $`T_\mathrm{c}`$ cannot exceed $`0.056R\mathrm{ln}2=0.32`$ J/mol K, where $`R`$ is the molar gas constant. The reason this value is the upper limit is that magnetic critical fluctuations will increase the specific heat, and hence the entropy, above $`T_\mathrm{c}`$ and thus reduce it at (and below) $`T_\mathrm{c}`$, by conservation of magnetic entropy, compared to the values for the isolated chain at the same reduced temperatures. Similarly, the $`C(T)`$ data in Fig. 8(a) allow one to estimate the minimum lattice specific heat contribution $`C^{\mathrm{lat}}(T)`$ above $`T_\mathrm{c}`$ if the $`C^{\mathrm{lat}}(T)`$ has not been determined previously from experiments and/or theory directly.

At low temperatures, the electronic specific heat coefficient $`C(T)/T`$ becomes independent of temperature (apart from logarithmic corrections, see below), as does the spin susceptibility $`\chi ^{}(t)`$, just as in a metal (Fermi liquid). Therefore it is of interest to compute a normalized ratio of these two quantities. For a metal, the relevant quantity is the Wilson-Sommerfeld ratio, which for

$`S=1/2`$ quasiparticles reads, in the notation of this paper and with $`k_\mathrm{B}`$ set to 1,

$$R_\mathrm{W}(t)=\frac{4\pi ^2\chi ^{}(t)t}{3C(t)}.$$
(98)
In a degenerate free electron gas, $`R_\mathrm{W}=1`$ and is independent of $`t`$. For exchange-enhanced metals $`1<R_\mathrm{W}10`$, for $`S=1/2`$ Kondo impurities in a metal the Wilson ratio associated with the impurities is $`R_\mathrm{W}=2`$, and for many heavy fermion metals $`R_\mathrm{W}2`$. Plotted in Fig. 10 is $`R_\mathrm{W}(t)`$ for the $`S=1/2`$ AF Heisenberg chain, where $`C(t)/t`$ and $`\chi ^{}(t)`$ were given above in Figs. 8(b) and 4, respectively. For $`t0`$, the Wilson ratio for the $`S=1/2`$ Heisenberg chain is exactly 2. With increasing $`t`$, $`R_\mathrm{W}`$ is seen to be nearly independent of $`t`$ to within $`\pm 10`$% up to $`t0.4`$, but the influence of the logarithmic corrections to both $`\chi (T)`$ and $`C(T)`$ are quantitatively important. Although the logarithmic corrections for $`\chi (T)`$ and $`C(T)`$ oppose each other in their ratio in $`R_\mathrm{W}(t)`$, the logarithmic corrections for $`\chi (t)`$ win out, giving a net $`10`$% increase in $`R_\mathrm{W}(t)`$ with increasing $`t`$ at low $`t`$. At higher $`t`$, the system crosses over to the expected local moment Heisenberg behavior where $`R_\mathrm{W}t^2`$. Thus as far as the thermodynamics is concerned, the uniform Heisenberg chain behaves at low temperatures as expected for a Fermi liquid, apart from the influence of the logarithmic corrections. This quasi-Fermi liquid behavior arises because the elementary excitations at low temperatures are $`S=1/2`$ spinons which are fermions with a Fermi surface (i.e., Fermi points in one dimension). Since the spinons carry no charge, the chain is an insulator. The deviation of the Wilson ratio from unity and the logarithmic corrections are due to spinon interactions.

#### 1 High-temperature series expansions

The HTSE for the specific heat of a spin $`S`$ AF uniform Heisenberg chain is

$$\frac{C(T)}{Nk_\mathrm{B}}=\frac{x^2}{3t^2}\left[1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{c_n(x)}{t^n}\right],$$
(100)
$$x=S(S+1),t=\frac{k_\mathrm{B}T}{J},$$
(101)
$`c_1={\displaystyle \frac{1}{2}},c_2={\displaystyle \frac{1}{15}}(38x3x^2),`$

$`c_3={\displaystyle \frac{1}{36}}(316x4x^2),`$

$`c_4={\displaystyle \frac{1}{5040}}(1921432x+1123x^2+800x^3+160x^4),`$

$$c_5=\frac{1}{21600}(4143768x+6635x^2+2624x^3+480x^4).$$
(102)
Specializing Eqs. (III B 1) to $`S=1/2`$ $`(x=3/4`$) then gives

$$\frac{C(T)}{Nk_\mathrm{B}}=\frac{3}{16t^2}\left[1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{c_n}{t^n}\right],$$
(104)
$$c_1=\frac{1}{2},c_2=c_3=\frac{5}{16},c_4=\frac{7}{256},c_5=\frac{917}{7680}.$$
(105)
The two $`C(T)`$ HTSE terms of order $`1/t^2`$ and $`1/t^3`$ in Eqs. (102) are in agreement with the general expression for the two lowest-order HTSE expansion terms for $`C(T)`$ of the $`S=1/2`$ alternating-exchange Heisenberg chain in Eq. (45) with alternation parameter $`\alpha =1`$.

In a later section, the Bethe ansatz $`C(T)`$ data will be fitted to obtain a function accurately representing the $`C(T)`$ of the $`S=1/2`$ AF uniform Heisenberg chain. In order that we are not required to change our fitting equations from those we use for fitting magnetic susceptibility data, the coefficients for the series inverted from that in Eq. (104) are required. We obtain

$$\frac{C(T)}{Nk_\mathrm{B}}=\frac{3}{16t^2}\left[1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{d_n}{t^n}\right]^1,$$
(107)
$$d_1=\frac{1}{2},d_2=\frac{9}{16},d_3=\frac{1}{8},d_4=\frac{7}{128},d_5=\frac{7}{1920}.$$
(108)
#### 2 Low-temperature logarithmic corrections

At first sight, from Fig. 8 there appear to be no singularities in the temperature dependence of the specific heat for the $`S=1/2`$ AF uniform Heisenberg chain. However, if the electronic specific heat coefficient $`C(T)/T`$ is examined in detail, one sees anomalous behavior at low temperatures. Shown as the top curve in Fig. 11(a) is a plot of the difference between the electronic specific heat coefficient and its zero-temperature value, $`\mathrm{\Delta }C(t)/Nk_\mathrm{B}t[C(t)(2/3)t]/(Nk_\mathrm{B}t)`$ for $`0t0.1`$ \[compare with Fig. 8(b)\]. From this figure, there is still nothing particularly strange about the data. However, upon further expanding the plot to study the range $`0t0.005`$ as shown in Fig. 11(b), we see that $`\mathrm{\Delta }C/Nk_\mathrm{B}t`$ is developing an infinite slope as $`t0`$. This is the signature of the existence of logarithmic corrections to the specific heat at temperatures $`t1`$, just as it was for the magnetic susceptibility.

Klümper, Lukyanov, and others have found a logarithmic correction to the low-$`t`$ limit in Eq. (91). Lukyanov’s exact asymptotic expansion for the free energy per spin in zero magnetic field is

$`f=J\mathrm{ln}2`$ $``$ $`{\displaystyle \frac{(k_\mathrm{B}T)^2}{3J}}\left[1+{\displaystyle \frac{3}{8}}g^3+𝒪(g^4)\right]`$ (109)
$``$ $`{\displaystyle \frac{3^{3/2}(k_\mathrm{B}T)^4}{10\pi J^3}}[1+𝒪(g)],`$ (111)
where $`g(t/t_0)`$ and $`t_0`$ are the same as given in Eqs. (84) and (86), respectively, and where $`g(t)`$ was plotted in Fig. 5. The specific heat at constant volume is calculated using $`C=T^2f/T^2`$, yielding

$`{\displaystyle \frac{C_{\mathrm{lt},\mathrm{g}}(T)}{Nk_\mathrm{B}}}={\displaystyle \frac{2k_\mathrm{B}T}{3J}}[1`$ $`+`$ $`{\displaystyle \frac{3}{8}}g^3+𝒪(g^4)]`$ (112)
$`+`$ $`{\displaystyle \frac{2(3^{5/2})}{5\pi }}\left({\displaystyle \frac{k_\mathrm{B}T}{J}}\right)^3[1+𝒪(g)].`$ (114)
This formula shows that the electronic specific heat coefficient $`C(T)/T`$ increases quadradically with $`T`$ at low $`T`$ (after subtracting the logarithmic corrections). The numerical prefactor of the $`t^3`$ term is $`1.98478\mathrm{}`$. If the approximate expansion for $`g()`$ in Eq. (87) is substituted into Eq. (114), one obtains

$`{\displaystyle \frac{C_{\mathrm{lt},\mathrm{log}}(T)}{Nk_\mathrm{B}}}=`$ $`{\displaystyle \frac{2k_\mathrm{B}T}{3J}}\left\{1+{\displaystyle \frac{3}{(2)^3}}+𝒪\left[{\displaystyle \frac{1}{(2)^4}}\right]\right\}`$ (115)
$`+`$ $`{\displaystyle \frac{2(3^{5/2})}{5\pi }}\left({\displaystyle \frac{k_\mathrm{B}T}{J}}\right)^3\left[1+𝒪\left({\displaystyle \frac{1}{2}}\right)\right],`$ (117)
where the prefactor 3/8 in the logarithmic correction term was found independently by Klümper, confirming

Refs.  and . The difference between $`C_{\mathrm{lt},\mathrm{log}}(T)`$ and $`C_{\mathrm{lt},\mathrm{g}}(T)`$ is plotted vs temperature in Fig. 12, where the difference becomes $`>10^{10}`$ only for $`t10^5`$.

Shown in Fig. 13 is the deviation $`\mathrm{\Delta }C/Nk_\mathrm{B}`$ ($``$) of the Bethe ansatz data from Lukyanov’s theoretical prediction in Eq. (114). For temperatures $`t10^4`$, the agreement is better than $`10^8`$. At higher temperatures, the uncertainty in the theoretical prediction due to the unknown $`𝒪(g^4)`$ and higher order correction terms becomes an important factor in the comparison. The length of the error bar on each data point in Fig. 13 has arbitrarily been set to $`(4/3)tg^4(t)`$ \[cf. Eq. (114)\]; the $`𝒪(g)`$ uncertainty in the $`T^3`$ term is negligible compared to this. Also plotted in Fig. 13 is the deviation of the numerical data from the extrapolated linear low-$`T`$ behavior ($``$). A comparison of the two data sets indicates that the $`𝒪(g^3)`$ logarithmic correction term is responsible for at least most of this latter difference for temperatures $`t0.001`$.

A more rigorous evaluation of the influence of the above logarithmic correction term is obtained by correcting for it in the plot of $`\mathrm{\Delta }C/t`$ vs $`t`$, as shown by the second curve from the top in each of Figs. 11(a) and 11(b). From the latter figure, we infer that although subtracting this correction term from the data helps to remove the zero-temperature singularity, a singularity is still present but with reduced amplitude. This means that additional logarithmic correction terms are important, within the accuracy and precision of the data. Another indication of this is shown in Fig. 14, where we have plotted $`\mathrm{\Delta }C/t^3`$ vs $`t`$. According to Eq. (114), after accounting for the logarithmic correction term(s), the result should be independent of $`t`$ at low $`t`$. Instead, both before and after accounting for the log correction term, there is a strong

upturn at low temperatures although the strength of the upturn is smaller after subtracting the influence of the log correction term.

The numerical Bethe ansatz specific heat data are sufficiently accurate and precise that we can estimate the coefficients of the next two logarithmic correction ($`g^4,g^5`$) terms in Eq. (114) from these data as follows. From Eq. (114), if we plot the numerical data as $`[C(t)/(Nk_\mathrm{B}t)(2/3)(1+3g^3/8)]/g^4`$ vs $`g`$ at low temperatures, where the $`t^3`$ term can be neglected, and fit the lowest $`t`$ data by a straight line, the $`y`$ intercept for $`g0`$ gives the coefficient of the $`g^4`$ term and the slope gives the coefficient of the $`g^5`$ term. This plot is given in Fig. 15. This type of plot places extreme demands on the accuracy of the data. Even so, we see that the data follow the required linear behavior even at the lowest temperatures. We fitted a straight line to the data from $`t=5\times 10^{25}`$ up to a maximum temperature $`t^{\mathrm{max}}`$. The fit parameters and rms deviation held nearly constant for $`t^{\mathrm{max}}=5\times 10^{15}`$ (11 data points) up to $`t^{\mathrm{max}}=5\times 10^8`$ (18 data points), but both quantities changed rapidly upon further increasing $`t^{\mathrm{max}}`$. The fit for $`t^{\mathrm{max}}=5\times 10^8`$ is shown as the straight line in Fig. 15, along with the fit parameters. From the parameters of the fit \[after accounting for the prefactor of 2/3 in Eq. (114)\], we include our estimated coefficients in Eq. (114) explicitly as

$`{\displaystyle \frac{C_{\mathrm{lt},\mathrm{g}}(T)}{Nk_\mathrm{B}}}={\displaystyle \frac{2k_\mathrm{B}T}{3J}}[1`$ $`+`$ $`{\displaystyle \frac{3}{8}}g^3+a_4g^4+a_5g^5+𝒪(g^6)]`$ (119)
$`+`$ $`{\displaystyle \frac{2(3^{5/2})}{5\pi }}\left({\displaystyle \frac{k_\mathrm{B}T}{J}}\right)^3[1+𝒪(g)],`$ (121)
$$a_4=1.5374(3),a_5=3.125(11).$$
(122)
The influences of these $`g^4`$ and $`g^5`$ logarithmic correction terms on the data in Figs. 11 and 14 are shown as the two additional data sets in each figure, where accounting for these two terms is seen to largely remove the remaining singular behavior as $`t0`$. From Fig. 14, we can now estimate that the coefficient of the $`t^3`$ term in Eq. (14) is a little larger than 2, contrary to the exact value 1.98478$`\mathrm{}`$. The magnitude of this difference is about as expected from the $`𝒪(g)`$ logarithmic correction to the $`t^3`$ term, since $`g(t0.1)0.1`$. The remaining upturn at low temperatures in Fig. 14 is due to residual logarithmic corrections which are not accounted for.

If the approximate expansion for $`g()`$ in Eq. (87) is inserted into Eq. (121), one obtains

$`{\displaystyle \frac{C_{\mathrm{lt},\mathrm{g}}}{Nk_\mathrm{B}}}`$ $`=`$ $`{\displaystyle \frac{2k_\mathrm{B}T}{3J}}\{1+{\displaystyle \frac{3}{(2)^3}}{\displaystyle \frac{9\mathrm{ln}()16a_4}{(2)^4}}`$ (123)
$`+`$ $`{\displaystyle \frac{\mathrm{ln}[18\mathrm{ln}()64a_49]+32a_5}{(2)^5}}+𝒪\left[{\displaystyle \frac{1}{(2)^6}}\right]\}`$ (124)
$`+`$ $`{\displaystyle \frac{2(3^{5/2})}{5\pi }}\left({\displaystyle \frac{k_\mathrm{B}T}{J}}\right)^3\left[1+𝒪\left({\displaystyle \frac{1}{2}}\right)\right].`$ (126)
## IV Fits to $`𝝌^{\mathbf{}}\mathbf{(}𝒕\mathbf{)}`$ and $`𝑪\mathbf{(}𝒕\mathbf{)}`$ of Heisenberg Spin Lattices

### A General $`𝝌^{\mathbf{}}\mathbf{(}𝒕\mathbf{)}`$ fit considerations

The general expression we use to fit theoretical numerical $`\chi ^{}(t)`$ data for $`S=1/2`$ Heisenberg spin lattices is

$$\chi ^{}(t)=\frac{\mathrm{e}^{\mathrm{\Delta }_{\mathrm{fit}}^{}/t}}{4t}𝒫_{(r)}^{(q)}(t),$$
(128)
$$𝒫_{(r)}^{(q)}(t)=\frac{1+_{n=1}^qN_n/t^n}{1+_{n=1}^rD_n/t^n},$$
(129)
where the orders $`q`$ and $`r`$ of the Padé approximant $`𝒫_{(r)}^{(q)}`$ are often constrained by the behavior of $`\chi ^{}(t)`$ at low $`t`$, and the fitted gap $`\mathrm{\Delta }_{\mathrm{fit}}^{}`$ is not necessarily the same as the true gap. At high $`t`$, $`\chi ^{}(t)`$ in Eqs. (IV A) approaches the Curie law $`1/(4t)`$ as required \[for a general spin $`S`$ lattice, the numerical prefactor 1/4 in Eq. (128) would be replaced by $`S(S+1)/3`$\].

The $`N_n`$ and $`D_n`$ parameters in Eq. (129) are not in general independent if one or more of the HTSE conditions in Eqs. (7) and (8) are invoked. For example, for $`n=1`$–3 one finds

$$D_1=(d_1+N_1)\mathrm{\Delta }_{\mathrm{fit}}^{},$$
(131)
$$D_2=(d_2+d_1N_1+N_2)\mathrm{\Delta }_{\mathrm{fit}}^{}(d_1+N_1)+\frac{\mathrm{\Delta }_{\mathrm{fit}}^{}{}_{}{}^{2}}{2},$$
(132)
$`D_3=(d_3`$ $`+`$ $`d_2N_1+d_1N_2+N_3)\mathrm{\Delta }^{}_{\mathrm{fit}}(d_2+d_1N_1+N_2)`$ (133)
$`+`$ $`{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{fit}}^{}{}_{}{}^{2}}{2}}(d_1+N_1){\displaystyle \frac{\mathrm{\Delta }_{\mathrm{fit}}^{}{}_{}{}^{3}}{6}}.`$ (134)
In general, one has

$$D_n=\underset{p=0}{\overset{n}{}}\underset{m=0}{\overset{np}{}}\frac{(\mathrm{\Delta }_{\mathrm{fit}}^{})^p}{p!}d_mN_{npm}.$$
(135)
A fit of experimental or theoretical $`\chi ^{}(t)`$ data by Eqs. (IV A) can be constrained by inserting one or more of Eqs. (IV A) and (135) into Eq. (129). These constraints are especially useful for high-$`t`$ extrapolations when $`\chi ^{}(t)`$ data are not available for high temperatures $`t1`$, and/or to reduce the number of fitting parameters required to obtain a fit of specified precision. In the following fits to the numerical $`\chi ^{}(t)`$ data for the dimer, the uniform chain, and finally our QMC and TMRG data for the alternating-exchange chain, the three constraints in Eqs. (IV A) on $`D_1,D_2`$, and $`D_3`$, respectively, are enforced in each case, where $`d_1,d_2`$, and $`d_3`$ for the alternating-exchange chain are given in Eq. (38).

All of the fits reported in this paper were carried out on a 400 MHz Macintosh G3 (B&W) computer with 1GB of RAM. Most fits were implemented using the program Mathematica 3.0, although a few of the simpler ones (fits to experimental data) were done using KaleidaGraph 3.08c. The fits using Mathematica sometimes required prodigious amounts of memory, e.g., 930 MB for the 28-parameter fit to the combined 2551 data point QMC and TMRG $`\chi ^{}(\alpha ,t)`$ data set for the alternating-exchange chain in Sec. V below.

### B Fit to $`𝝌^{\mathbf{}}\mathbf{(}𝒕\mathbf{)}`$ of the $`𝑺\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$ antiferromagnetic Heisenberg dimer

The spin gap of the $`S=1/2`$ Heisenberg dimer is $`\mathrm{\Delta }=J`$, where $`J`$ is the antiferromagnetic exchange constant within the dimer. The spin susceptibility and its low-temperature limit are given by Eqs. (12). The $`\chi ^{}(t)`$ is plotted in Fig. 1(a) for $`0.02t4.99`$. In order to later obtain a continuous fit function for $`\chi ^{}(\alpha ,t)`$ for the entire range $`0\alpha 1`$ of the alternating-exchange chain, it is necessary to first obtain a high accuracy fit to the exact expression (14) for the dimer by our general fitting function in Eqs. (IV A), in addition to Fit 1 obtained for the uniform chain below. The form of our fit function in Eqs. (IV A) allows both the low- and high-$`t`$ limiting forms of $`\chi ^{}(t)`$ for the dimer to be exactly reproduced. The low-$`t`$ limit in Eq. (15) requires that $`r=q`$ and that $`D_q=N_q/4`$ in the Padé approximant $`𝒫_{(r)}^{(q)}`$; we also take $`\mathrm{\Delta }_{\mathrm{fit}}=\mathrm{\Delta }`$, so the total number of fitting parameters is $`2q4`$.

We fitted the 498-point double-precision representation of $`\chi ^{}(t)`$ in Fig. 1(a) from $`t=0.02`$ to $`t=4.99`$ by Eqs. (IV A) using the above constraints. The variances of the four fits for $`q=r=4,`$ 5, 6 and 7 were $`2.5\times 10^{13}`$, $`1.17\times 10^{16}`$, $`5.3\times 10^{17}`$, and $`5.6\times 10^{19}`$, respectively, showing that Eqs. (IV A) have the potential for very high accuracy fits with a relatively small number of fitting parameters. The six $`N_n`$ ($`n=1`$–5) and $`D_4`$ parameters of the fit for $`q,r=5`$ are given in Table I, along with $`D_1,D_2`$, and $`D_3`$ computed from Eqs. (IV A) and $`D_5=N_5/4`$. The Padé approximant $`𝒫_{(5)}^{(5)}`$ in the fit function has no poles or zeros on the positive $`t`$ axis. The fit is shown by the solid curve in Fig. 1(a), and the deviation of the fit from the exact susceptibility in Eq. (14) is plotted versus $`t`$ in Fig. 1(b).

### C Fits to $`𝝌^{\mathbf{}}\mathbf{(}𝒕\mathbf{)}`$ of the $`𝑺\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$ Antiferromagnetic Uniform Heisenberg Chain

##### Fit 1: $`0.01t5`$.

Fits to the uniform chain $`\chi ^{}(t)`$ calculated by Eggert, Affleck and Takahashi for limited temperature regions were obtained previously. Here we obtain a fit (Fit 1) to the higher accuracy data of Klümper and Johnston for the temperature region $`0.01t5`$ (999 data points) using Eqs. (IV A), the results of which will be utilized later in the fit function for $`t0.01`$ for our QMC and TMRG alternating-exchange chain $`\chi ^{}(\alpha ,t)`$ data. This uniform chain fit can be accurately extrapolated to arbitrarily high $`t`$.

The requirement that $`\chi ^{}(t0)`$ is a finite non-zero value requires $`\mathrm{\Delta }_{\mathrm{fit}}^{}=0`$ and $`r=q+1`$ in Eqs. (IV A). We found that using $`q=5`$ and $`r=6`$ produces a fit sufficiently accurate for use in the fit function for our QMC and TMRG calculations for the alternating chain. The seven $`N_n`$ ($`n=1`$–5) and $`D_n`$ ($`n=4`$–6) parameters obtained for the fit with $`q=5,r=6`$ are given in the column labeled “Fit 1” in Table I, along with $`D_1`$, $`D_2`$,

and $`D_3`$ computed from Eqs. (IV A). The Padé approximant $`𝒫_{(6)}^{(5)}`$ in the fit function has no poles or zeros on the positive $`t`$ axis. The deviation of the fit from the data is plotted in Fig. 16. The variance of the fit is $`2.97\times 10^{12}`$, and the relative rms deviation of the fit from the data in the fitted $`t`$ region is 14.5 ppm. Extrapolation of the fit to higher temperatures is very accurate.

The quality of Fit 1 does not approach the limitation imposed by the absolute accuracy of the data ($`1\times 10^9`$). For an ideal fit, the variance is expected to be $`10^{18}`$ and the relative rms deviation $`0.01`$ ppm. As can be inferred from Fit 2 in the following section, the reason that Fit 1 cannot be optimized to this extent is due to the $`t=0`$ critical point and associated logarithmic divergence in the slope of $`\chi ^{}(t)`$ as $`t0`$; this divergence cannot be fitted accurately by a finite polynomial or Padé approximant. We attempted to improve the accuracy of the fit over the same temperature range $`0.01t5`$ by replacing the Padé approximant $`𝒫_{(6)}^{(5)}`$ in the fit function by $`𝒫_{(7)}^{(6)}`$, which incorporates two additional fitting parameters. The variance improved somewhat to $`2.18\times 10^{12}`$ and the relative rms deviation improved slightly to 12.2 ppm, but the Padé approximant developed a pole at $`1/t=129.23`$, and hence this fit was discarded. Although the temperature at which this pole occurs is below the fitted temperature range, as a general rule we cannot allow poles in the fit function at low temperatures because of problems that can occur when using the fit function to model experimental data which include data at temperatures lying below the fitted temperature range of the fit function. In fact, we will encounter this situation frequently in modeling experimental data later. For exam-

ple, for $`\mathrm{NaV}_2\mathrm{O}_5`$, $`t=0.01`$ corresponds to an absolute temperature $`T7`$ K, whereas the experimental data and modeling extend down to $`2`$ K.

##### Fit 2: $`0t5`$.

We can greatly improve the accuracy of the fit compared to that of Fit 1, and extend the fit to $`t=0`$, by restricting the high-temperature limit of the fit and using in the fit function one or more low-temperature logarithmic correction terms discussed in Sec. III A 2. In particular, in this section we obtain a very high precision fit (Fit 2) to the exact $`t=0`$ value and to the calculations of Klümper and Johnston over the entire temperature range $`5\times 10^{25}t5`$ of the calculations. We do not use this fit in our formulation of the fit function for the alternating-exchange chain. However, Fit 2 will be generally useful for evaluating the accuracy of other theoretical calculations of $`\chi ^{}(t)`$ for the uniform chain, such as our TMRG calculations to be presented below, and for modeling appropriate experimental $`\chi (T)`$ data whose scaled upper temperature limit is below $`t=5`$.

We initially formulated a fit function utilizing a modified Padé approximant in which the last term of the numerator and/or denominator contained the $`\chi _{\mathrm{lt},\mathrm{log}}^{}`$ expansion in Eq. (90), such that the low-temperature expansion of the fit function yielded $`\chi _{\mathrm{lt},\mathrm{log}}^{}`$ to lowest orders in $`t`$. The best fit to the data from $`t=5\times 10^{25}`$ to 2.5 (777 data points) was unsatisfactory, with a variance $`v=2.4\times 10^{11}`$ and a relative rms deviation $`\sigma _{\mathrm{rms}}=45`$ ppm. Allowing an arbitrary $`t^2`$ coefficient in place of the exact value $`\sqrt{3}/\pi ^3`$ yielded an improved fit with $`v=1.1\times 10^{12}`$ and $`\sigma _{\mathrm{rms}}=9.6`$ ppm. However, this fit was still unsatisfactory, given the high absolute accuracy ($`1\times 10^9`$) of the data. From these results it became clear that a fit function which can fit the data to much higher accuracy over such a large temperature range would indeed have to include an expression $`\chi _{\mathrm{log}}^{}(t)`$ containing logarithmic correction terms, but where the form and/or coefficient of one or more of these terms would have to be empirically determined by trial and error. This process yielded the formulation we now describe.

The $`\chi _{\mathrm{log}}^{}(t)`$ function is incorporated into our fit function in Eqs. (IV A) as follows. As in Fit 1, the finite value of $`\chi ^{}(0)`$ requires $`\mathrm{\Delta }_{\mathrm{fit}}=0`$ in Eq. (128) and $`r=q+1`$ in the Padé approximant $`𝒫_{(r)}^{(q)}(t)`$ in Eq. (129). Since the two terms highest order in $`1/t`$ in $`𝒫_{(r)}^{(q)}(t)`$ (one each in the numerator and denominator) dominate the fit as $`t0`$ and become small for $`t1`$, relative to the other terms in the numerator and denominator, respectively, we incorporate $`\chi _{\mathrm{log}}^{}(t)`$ into the last term in the numerator of a modified $`𝒫_{(r)}^{(q)}(t)`$. Trial fits showed that to obtain the optimum accuracy of the fit required $`q=8`$ and $`r=9`$.

Our final fit function for Fit 2 is

$$\chi ^{}(t)=\left(\frac{1}{4t}\right)\frac{1+\left[_{n=1}^7N_n/t^n\right]+4N_8\chi _{\mathrm{log}}^{}(t)/t^8}{1+\left[_{n=1}^8D_n/t^n\right]+N_8/t^9},$$
(137)
$$\chi _{\mathrm{log}}^{}(t)=\frac{1}{\pi ^2}\left[1+\frac{1}{2}\frac{\mathrm{ln}\left(+\frac{1}{2}\right)N_{81}}{(2)^2}+\frac{N_{82}}{(2)^3}\right],$$
(138)
$$\mathrm{ln}(t_1/t),$$
(139)
subject to the three constraints on $`D_1,D_2`$, and $`D_3`$ in Eqs. (IV A) which are required by the HTSE. Two of the four logarithmic correction terms in Eq. (138) are identical to the first two such terms in Eq. (78). By construction, the exact $`\chi ^{}(0)=1/\pi ^2`$ is fitted exactly.

We fitted all of the numerical $`\chi ^{}(t)`$ data, calculated over the range $`5\times 10^{25}t5`$ (1119 data points), by Eqs. (IV C). The 19 fitting parameters of the fit function (IV C), which are $`N_n`$ ($`n=1`$–8), $`D_n`$ ($`n=4`$–8), $`N_{81},N_{82}`$ and $`t_1`$, are given in the column labeled “Fit 2” in Table I, along with $`D_1,D_2`$, and $`D_3`$ computed from Eqs. (IV A). The data to parameter ratio is 59. The denominator of the modified Padé approximant in Eq. (137) has no zeros for any real positive $`t`$. The fit is shown in the low-temperature region $`0t0.02`$ in Fig. 4(b) \[over the larger $`t`$ range plotted in Fig. 4(a), the fit is indistinguishable from the data and is therefore not plotted there\].

The deviation of Fit 2 from the numerical data for $`10^{25}t5`$ is plotted vs $`\mathrm{log}_{10}t`$ in Fig. 17(a), and an expanded plot at the higher temperatures is shown in Fig. 17(b). Due to a logarithmic divergence in $`\chi _{\mathrm{log}}^{}(t)`$ at $`t=t_1=5.696`$, Fit 2 should not be used (e.g., for modeling experimental data) at temperatures $`t5`$. The variance of the fit is $`9.8\times 10^{17}`$, and the relative rms deviation is $`\sigma _{\mathrm{rms}}=0.087`$ ppm. These values are both much smaller than for Fit 1 above. The relatively large number of fitting parameters in Fit 2 is justified a posteriori by the extremely high quality of the fit over a temperature range spanning 25 orders of magnitude.

### D Fit to $`𝑪\mathbf{(}𝒕\mathbf{)}`$ for the $`𝑺\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$ antiferromagnetic uniform Heisenberg chain

The logarithmic corrections to the magnetic specific heat $`C(t)`$ at low temperatures, discussed above in Sec. III B 2, do not pose as serious a problem for fitting the data as for $`\chi ^{}(t)`$, because the strength of these log corrections is much smaller for $`C(t)`$ than for $`\chi ^{}(t)`$. In addition, since here we fit $`C(t)`$, and not the electronic specific heat coefficient $`C(t)/t`$, the influence of the log corrections is ameliorated by the multiplicative leading order $`t^1`$ dependence of $`C(t)`$. Even so, in order to obtain the optimum fit to the highly accurate Bethe ansatz $`C(t)`$ data, we found it necessary to take the influence of the logarithmic corrections into account.

Our fit to the Bethe ansatz $`C(t)`$ data, some of which were shown previously in Fig. 8(a), was carried out in two stages. First, the data from $`t=0.01`$ to the maximum temperature $`t=5`$ of the calculations were fitted by the Padé approximant $`𝒫_{(r)}^{(q)}`$ in Eq. (129) with a prefactor $`3/(16t^2)`$ to satisfy the HTSE in Eqs. (III B 1) to lowest order in $`1/t`$. The orders $`q`$ and $`r`$ of $`𝒫_{(r)}^{(q)}`$ were chosen to satisfy $`r=q+3`$ so that $`C(t0)t`$. To obtain a fit of the required accuracy (see the fit deviations given below) we found that $`q=6`$ and $`r=9`$ are of sufficiently high order. Due to the presence of the log corrections at very low $`t`$, we did not require the parameters $`N_6`$ and $`D_9`$ to yield the exact coefficient $`\gamma =2/3`$ in the expression $`C(t)/Nk_\mathrm{B}=\gamma t`$, in a low-$`t`$ expansion of the fit function. We also found that to obtain the best fit, only the one additional HTSE constraint (on $`D_1`$) in Eq. (131) (with $`\mathrm{\Delta }_{\mathrm{fit}}^{}=0`$) could be used. It was quite difficult to find the region in parameter space in which the absolute minimum in the variance of the fit resided; the initial starting parameters usually flowed to regions with local variance minima in them with much larger values (by two to four orders of magnitude) than the smallest variance we ultimately found. Then the deviation of the fit from all the data for $`5\times 10^{25}t5`$ was computed. The fit deviations for $`t0.01`$ were very small \[$`𝒪(10^8)`$\], but the log corrections which become most important at lower

temperatures resulted in fit deviations at $`t<0.01`$ an order of magnitude larger. We therefore fitted the fit deviation versus $`t`$ for $`0<t0.1`$ by a separate empirically determined function $`F(t)`$, so the net fit function consists of the Padé approximant fit function minus the fit function to the low-$`t`$ fit deviations. In the final fitting cycles the two functions were refined simultaneously.

Our final fit function for $`C(t)`$ in the range $`0t5`$ is

$$\frac{C(t)}{Nk_\mathrm{B}}=\frac{3}{16t^2}𝒫_{(9)}^{(6)}(t)F(t),$$
(141)
$$𝒫_{(9)}^{(6)}(t)=\frac{1+\left[_{n=1}^6N_n/t^n\right]}{1+\left[_{n=1}^9D_n/t^n\right]},$$
(142)
$$F(t)=a_1t\mathrm{sin}\left(\frac{2\pi }{a_2+a_3t}\right)\mathrm{e}^{a_4t}+a_5t\mathrm{e}^{a_6t},$$
(143)
subject to the constraint on $`D_1`$ in Eq. (131) which is required by the HTSE. By construction, the exact $`C(0)=0`$ is fitted exactly. The 20 fitting parameters, $`N_n(n=16),D_n(n=29)`$ and $`a_n(n=16)`$, are given in Table I, together with the constrained parameter $`D_1`$ computed from Eq. (131) with $`\mathrm{\Delta }_{\mathrm{fit}}^{}=0`$ and $`d_1`$ given in Eq. (108). The deviation of the fit from the data is shown in a semilog plot vs temperature in Fig. 18. The maximum deviations of $`\pm 4\times 10^8`$ occur at $`t0.3`$. The absolute rms deviation of the fit from all the data (1119 data points), which extend over the temperature range $`5\times 10^{25}t5`$, is $`1.34\times 10^8`$, and the relative rms deviation for $`0.01t5`$ (999 data points) is 0.50 ppm.

At high temperatures, our $`C(t)`$ fit function reduces by construction to the lowest order $`1/t^2`$ and $`1/t^3`$ terms of the HTSE of $`C(t)`$ in Eqs. (III B 1), so extrapolation of our $`C(t)`$ fit function to arbitrarily higher temperatures should be very accurate (see Fig. 18). In particular, even though our fit was to $`C(t)`$ and hence not optimized as a fit to the electronic specific heat coefficient $`C(t)/t`$, the magnetic entropy $`S`$ at $`t=\mathrm{}`$ computed from our $`C(t)`$ fit function is

$$\frac{S(t=\mathrm{})}{Nk_\mathrm{B}}_0^{\mathrm{}}\frac{C(t)}{Nk_\mathrm{B}t}𝑑t=\mathrm{0.693\hspace{0.17em}147\hspace{0.17em}235},$$
(144)
which is the same as the exact value $`\mathrm{ln}2=\mathrm{0.693\hspace{0.17em}147\hspace{0.17em}181}`$ to within 8 parts in $`10^8`$. This agreement reflects well on our fit function, and of course also strongly confirms the high accuracy of the Bethe ansatz $`C(t)`$ data.

### E Fit function for the $`𝑺\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$ AF alternating-exchange Heisenberg chain $`𝝌^{\mathbf{}}\mathbf{(}𝜶\mathbf{,}𝒕\mathbf{)}`$

Here we formulate a single two-dimensional ($`\alpha ,t`$) function to accurately fit numerical calculations of $`\chi ^{}(\alpha ,t)`$ for the $`S=1/2`$ alternating-exchange Heisenberg chain for the entire range $`0\alpha 1`$, and for the entire temperature range $`t0.01`$ over which our Fit 1 for $`\chi ^{}(t)`$ of the uniform chain is most accurate, subject to four general requirements as follows. (i) The HTSE of the $`\chi ^{}(\alpha ,t)`$ fit function must give the correct result to $`𝒪(1/t^4)`$, as satisfied by the fit functions for the dimer and uniform chain (Fit 1) susceptibilities above, so that the fit can be accurately extrapolated to higher temperatures. (ii) We require the $`\chi ^{}(\alpha ,t)`$ fit function to become identical with those found above for the isolated dimer and for the uniform chain (Fit 1) when $`\alpha =0`$ and $`\alpha =1`$, respectively. As discussed above in Sec. II C, at any finite temperature, $`\overline{\chi ^{}}(\delta ,\overline{t})`$ in the variables $`\delta `$ and $`\overline{t}`$ is an even (analytic) function of $`\delta `$. Therefore, as a minimum accommodation of this fact, (iii) we require that the fit function for $`\chi ^{}(\alpha ,t)`$, when transformed to the form $`\overline{\chi ^{}}(\delta ,\overline{t})`$, must have the property $`\overline{\chi ^{}}(\delta ,\overline{t})/\delta |_{\delta =0}=0`$ at all finite temperatures. This requirement is clearly the minimum necessary in order to accurately interpolate the fit vs $`\alpha `$ for $`\alpha 1`$ at each $`t`$, and to thereby accurately model the susceptibility of materials which are in or near this limit. Finally, the QMC and TMRG calculations of $`\chi ^{}(\alpha ,t)`$ to be presented below are sufficiently accurate and cover sufficiently large ranges of $`\alpha `$ and $`t`$ with sufficient resolution that (iv) we require the nonanalytic energy gap $`\mathrm{\Delta }(\alpha )`$ \[see Eqs. (II C 2) and (49)\] to be included in the fit function in order to fit the data for $`\alpha 1`$ at $`t1`$, so as to avoid the alternate necessity of including high-order power series in $`\alpha `$ and $`t`$ in the fit function. We note that according to Eq. (48) or (49), $`\overline{\mathrm{\Delta }^{}}(\delta )/\delta |_{\delta =0}=\mathrm{}`$. The major obstacle we faced in formulating the fit function for $`\chi ^{}(\alpha ,t)`$ was to simultaneously satisfy both requirements (iii) and (iv), which at first sight seem to require mutually exclusive forms for the fit function.

We found that these four requirements can all be satisfied by an extension of the form of the fit function in Eqs. (IV A) which was used above for the isolated dimer and for the uniform chain Fit 1. This extension consists of using a modified Padé approximant $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}`$ in the fit function in place of the former $`𝒫_{(r)}^{(q)}`$. The fit function is

$$\chi ^{}(\alpha ,t)=\frac{\mathrm{e}^{\mathrm{\Delta }_{\mathrm{fit}}^{}(\alpha )/t}}{4t}𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}(\alpha ,t),$$
(146)
$`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}(\alpha ,t)=`$

$$\frac{[_{n=0}^6N_n/t^n]+(N_{71}\alpha +N_{72}\alpha ^2)(\mathrm{\Delta }_0^{}/t)^y/t^7}{[_{n=0}^7D_n/t^n]+(D_{81}\alpha +D_{82}\alpha ^2)(\mathrm{\Delta }_0^{}/t)^z\mathrm{e}^{(\mathrm{\Delta }_0^{}\mathrm{\Delta }_{\mathrm{fit}}^{})/t}/t^8},$$
(147)
$$\mathrm{\Delta }_{\mathrm{fit}}^{}(\alpha )=1\frac{1}{2}\alpha 2\alpha ^2+\frac{3}{2}\alpha ^3,$$
(148)
$`\mathrm{\Delta }_0^{}(\alpha )=(1`$ $``$ $`\alpha )^{3/4}(1+\alpha )^{1/4}`$ (149)
$`+`$ $`g_1\alpha (1\alpha )+g_2\alpha ^2(1\alpha )^2,`$ (150)
$$N_0=D_0=1,$$
(151)
$$N_n(\alpha )=\underset{m=0}{\overset{4}{}}N_{nm}\alpha ^m(n=16),$$
(152)
$$D_n(\alpha )=\underset{m=0}{\overset{4}{}}D_{nm}\alpha ^m(n=17).$$
(153)
To satisfy requirement (i), $`D_1(\alpha ),D_2(\alpha )`$, and $`D_3(\alpha )`$ are determined from the $`N_1(\alpha ),N_2(\alpha )`$, and $`N_3(\alpha )`$ fitting parameters according to the three constraints in Eqs. (IV A) demanded by the HTSE. In order to satisfy requirement (ii), the $`\{N_{n0},D_{n0}\}`$ parameters are set to be identical with those determined above for the dimer, and we require $`\{N_n(1),D_n(1)\}`$ to be identical with the corresponding fit parameters determined above in Fit 1 for the uniform chain. In order to satisfy requirement (iii), the $`N_{nm}`$ and $`D_{nm}`$ coefficients must satisfy

$$\underset{m=0}{\overset{4}{}}(n2m)(N_{nm}\mathrm{or}D_{nm})=0,$$
(154)
so that no $`\delta ^1`$ term appears in the Taylor series expansions in $`\delta `$ of the transformed $`\{\overline{N}_n(\delta ),\overline{D}_n(\delta )\}`$. These various constraints on the $`\{N_{nm},D_{nm}\}`$ parameters reduce the number of independent fitting parameters within this set from 50 to 20. Together with the parameters $`N_{71},N_{72},N_{81},N_{82},y,z`$ in Eq. (147) and $`g_1,g_2`$ in Eq. (150), the total number of independent fitting parameters in the fit function is 28.

The quantity $`\mathrm{\Delta }_{\mathrm{fit}}^{}(\alpha )`$ in the exponential prefactor to $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}`$ in Eq. (146) cannot be set equal to the true nonanalytic gap $`\mathrm{\Delta }^{}(\alpha )`$, because this prefactor affects the fit at all $`t`$, and would not allow requirement (iii) above to be fulfilled. In addition, the nonanalytic critical behavior of $`\mathrm{\Delta }^{}(\alpha 1)`$ in practice only becomes manifest in $`\chi ^{}(\alpha ,t)`$ at low temperatures $`t1`$. Therefore, we separated the spin gap into an analytic part $`\mathrm{\Delta }_{\mathrm{fit}}^{}(\alpha )`$ which goes into the argument of the exponential prefactor in Eq. (146), and a nonanalytic part $`\mathrm{\Delta }_0^{}(\alpha )`$ \[satisfying requirement (iv)\] which is placed into the argument of the exponential in the last term of the denominator of $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}`$ in Eq. (147) and which therefore only becomes important at low temperatures. The first two terms of $`\mathrm{\Delta }_{\mathrm{fit}}^{}(\alpha )`$ (to order $`\alpha ^1`$) in Eq. (148) are the first two terms of the exact dimer series expansion up to $`𝒪(\alpha ^9)`$ given by Barnes, Riera and Tennant for the AF alternating-exchange chain, and the last two are included so that $`\overline{\mathrm{\Delta }_{\mathrm{fit}}^{}}(\delta )/\delta |_{\delta =0}=0`$, in accordance with requirement (iii). The nonanalytic $`\mathrm{\Delta }_0^{}(\alpha )`$ in Eq. (150) contains the behavior in Eq. (47) proposed by Barnes, Riera, and Tennant, plus two analytic terms which are included to adjust the $`\alpha `$ dependence for $`\alpha 1`$ but which make no contribution at $`\alpha =0`$ or $`\alpha =1`$. Provided that the inequality $`y,z>4/3`$ is satisfied by the powers $`y`$ and $`z`$ in Eq. (147), the last term in each of the numerator and denominator of $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}(\alpha ,t)`$, when transformed to the variables $`(\delta ,\overline{t})`$, has a partial derivative with respect to $`\delta `$ which is zero at $`\delta =0`$.

We have now shown that at $`\delta =0(\alpha =1)`$, the partial derivative of each part of $`\overline{\chi ^{}}(\delta ,\overline{t})`$ with respect to $`\delta `$ is zero (if $`y,z>4/3`$, which is confirmed in the actual fit later). Hence, the entire fit function has the property $`\overline{\chi ^{}}(\delta ,\overline{t})/\delta |_{\delta =0}=0`$ at all finite temperatures, thus satisfying requirement (iii), despite the fact that the fit function contains the nonanalytic $`\mathrm{\Delta }_0^{}(\alpha )`$ as required by requirement (iv).

At the lowest temperatures, the last term in each of the numerator and denominator of $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}`$ in Eq. (147) should dominate the fit, together with the exponential prefactor to $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}`$ in Eq. (146), so in this limit our fit function for $`0<\alpha <1`$ becomes

$$\chi ^{}(\alpha ,t0)=\frac{N_{71}\alpha +N_{72}\alpha ^2}{4(D_{81}\alpha +D_{82}\alpha ^2)}\left[\frac{\mathrm{\Delta }_0^{}(\alpha )}{t}\right]^{yz}\mathrm{e}^{\mathrm{\Delta }_0^{}(\alpha )/t}.$$
(155)
This expression has the form of Eq. (18) (with $`\gamma =yz`$) as required in the low-$`t`$ limit. In fact, the forms of the last term in each of the numerator and denominator of $`𝒫_{\mathrm{m}}^{}{}_{(8)}{}^{(7)}`$ were designed to result in the form of Eq. (18) in the low-$`t`$ limit, with $`\mathrm{\Delta }_0^{}`$ and $`t`$ entering the prefactor only as their ratio as in Eq. (57), in addition to being consistent with requirements (iii) and (iv). One might expect the fitted $`y`$ and $`z`$ powers to satisfy $`yz=\gamma =1/2`$ as in Eq. (19). However, if a fit of $`\chi ^{}(t)`$ data by Eq. (18) is not carried out completely within the low-$`t`$ limit, an effective exponent $`\gamma 1`$ is often inferred \[see, e.g., Eq. (181) and subsequent discussion, and Fig. 36 below\]. Similarly, since many of our calculated $`\chi ^{}(\alpha ,t)`$ data sets for different $`\alpha `$ in the fitted temperature range $`t0.01`$ are not, or do not contain extensive data, in the low-$`t`$ limit, we did not impose the constraint $`yz=1/2`$. On the basis of the above discussion we expect the actual fitted values of $`y`$ and $`z`$ to yield $`yz1`$. In fact, as will be seen in the next section, our fitted parameters $`y`$ and $`z`$ give $`yz=1.14`$.

## V QMC and TMRG $`𝝌^{\mathbf{}}\mathbf{(}𝜶\mathbf{,}𝒕\mathbf{)}`$ Calculations and Fit for the $`𝑺\mathbf{=}\mathrm{𝟏}\mathbf{/}\mathrm{𝟐}`$ AF Alternating-Exchange Heisenberg Chain

QMC simulations of $`\chi ^{}(\alpha ,t)`$ were carried out on $`S=1/2`$ alternating-exchange chains containing 100 spins for $`\alpha =0.05,0.1,0.15,\mathrm{},`$ 0.9, 0.92, 0.94, 0.96, 0.97, 0.98, and 0.99 in various temperature ranges spanned by $`0.01t4`$.

Complementary TMRG calculations of $`\chi ^{}(\alpha ,t)`$ of $`S=1/2`$ alternating-exchange chains were carried out for $`\alpha =0.80`$, 0.82, …, 0.96, 0.97, 0.98, 0.99, 0.995 and 1, where the number of states kept was $`m=150`$ or 256. The calculations were carried out for temperatures given by $`1/t=0.1,0.2,\mathrm{},(1/t)^{\mathrm{max}}`$, with $`(1/t)^{\mathrm{max}}500`$ increasing with increasing $`\alpha `$, and comprised a total of 22 370 ($`\alpha ,t`$) parameter combinations. The details of the calculational method are given in Refs.  and . It should be noted that the TMRG calculations by their nature are explicitly in the thermodynamic limit.

The reason for doing TMRG calculations for the uniform chain ($`\alpha =1`$) was to enable comparison of the results with the values computed with the Bethe ansatz which have a high absolute accuracy of $`1\times 10^9`$. This comparison was done using the above very accurate and precise Fit 2 for the Bethe ansatz data. The relative deviation of the TMRG data from Fit 2 is shown in Fig. 19(a), and an expanded plot for the higher temperature region $`t0.01`$ is shown in Fig. 19(b). This comparison indicates that the accuracy of the TMRG calculations for both $`m=150`$ and 256 in the range $`t0.01`$ is better than 0.1 %, which is the same as the estimate made previously for $`m=80`$. However, the accuracy of these calculations deteriorates rapidly at lower $`t`$, to about 3 % at the lowest temperatures $`t0.002`$ for $`m=150`$.

Since the TMRG calculations extend close to the $`t=0`$ limit for most of the above-stated $`\alpha `$ values, the spin gaps can be estimated from these data. Comparisons with previous work can then be made of the dependence of the spin gap on $`\alpha `$. An important question, not yet

answered in previous work, is the approximate $`\alpha `$ value at which the asymptotic critical region is entered upon approaching the uniform limit. Performing these estimates and comparisons will be postponed to the following sections. In the present section, we present the QMC and TMRG $`\chi ^{}(\alpha ,t)`$ data and obtain a fit to these combined data by the fit function formulated in the previous section.

Some of the results for $`t2`$ are shown as the filled symbols without error bars in Fig. 20(a) (the error bars are smaller than the data symbols); an expanded plot of data for $`t0.4`$ is shown in Fig. 20(b). \[A log-log plot of the TMRG $`\chi ^{}(\alpha ,t)`$ data at low $`t`$ is shown below in Fig. 27.\] Also shown in both figures as the two bounding solid curves with no data points are the fits we obtained above to $`\chi ^{}(t)`$ for the dimer and uniform chain (Fit 1), respectively. The data points plotted for a given $`\alpha `$ value are the subset below the upper temperature limits of the figures, of the subset of available data points which were fitted by our fit function as described below.

We fitted a combined QMC and TMRG $`\chi ^{}(\alpha ,t)`$ data set containing 2551 selected data points over the temperature range $`0.01t10`$. The 802 QMC data points covered the ranges $`0.01t4`$ and $`0.05\alpha 0.99`$.

The average estimated absolute accuracy of these QMC data is $`1.7\times 10^4`$. The best estimated accuracy among these QMC data is $`7.7\times 10^6`$ and the worst is $`1.5\times 10^3`$, with the better accuracies occurring at the highest temperatures. The 1749 TMRG data points covered the ranges $`0.01t10`$ and $`0.8\alpha 0.995`$. We did not use all 22 370 TMRG data points in the available data set, because this would have weighted the region $`\alpha 1`$ too heavily in the fit, and in any case a large fraction of these are for temperatures below our low-temperature fitting limit of $`t=0.01`$. We used the low-temperature data to determine the spin gaps as described in the following section.

We fitted this $`\chi ^{}(\alpha ,t)`$ data set by Eqs. (IV E), with the constraints on the parameters discussed above. Obtaining a reliable 28-parameter two-dimensional fit to these data over the full above-cited ranges of $`t`$ and $`\alpha `$, with no poles in the fit, posed a very difficult challenge. The particular choice of starting parameters and the detailed sequence of refinements were found to be important to avoiding poles in the final fit. Since there are a total of 28 parameters in the fit function for 2551 data points, the data to parameter ratio is 91. The number of fitting parameters seems large, until it is realized that we are simultaneously fitting $`\chi ^{}(t)`$ data for 29 different $`\alpha `$ values, so on average a $`\chi ^{}(t)`$ data set for a given $`\alpha `$ value is fitted by a single parameter. A weighting function was not included during the variance minimization, because we were interested in obtaining a fit which treated all the data points the same on an absolute scale; this choice optimizes the fit for use in modeling experimental data.

The parameters of the fit are given in Table II, where we have also included the constrained parameters for completeness and for ease of implementation of our fit function by the reader. From Eqs. (IV A), the constrained parameters $`D_2`$ and $`D_3`$ contain products of the third-order (in $`\alpha `$) polynomial $`\mathrm{\Delta }_{\mathrm{fit}}^{}`$ with itself and/or with the fourth-order $`N_1`$ fitting polynomial, so $`D_2`$ and $`D_3`$ are of seventh and tenth-order, respectively. The two-dimensional fit is shown as the set of solid curves in Fig. 20. The variance of the fit is $`v=3.77\times 10^8`$. The absolute rms deviation $`\sqrt{v}1.9\times 10^4`$ is about the same as the average estimated accuracy of the QMC data noted above, indicating that the fit function is appropriate and that the fit is a reliable representation of

the data. The fit deviations from the 802 QMC and 1749 TMRG data are shown separately in Figs. 21(a) and 21(b), respectively. A comparison of the two figures shows that the TMRG data are, on average, significantly more precise at a given temperature.

After the parameters in the present $`\chi ^{}(\alpha ,t)`$ fit function were finalized, as a check on the accuracy of the fit function for $`\alpha `$ values close to the uniform limit, we carried out QMC $`\chi ^{}(t)`$ simulations for alternating-exchange chains of length $`L=400`$ and 800, factors of four and eight longer than the chains for which QMC data were combined with TMRG data to determine the fit function, respectively. The simulations were carried out for $`\alpha =0.98,0.985,0.99`$, and 0.995 at temperatures $`0.01t4`$. Overall, the fit function was found to be in extremely good agreement with the QMC data. For $`0.4t4`$, the $`\chi ^{}(\alpha ,t)`$ fit function agreed with the simulation data to within about $`\pm 5\times 10^5`$ or better. The deviations of the fit function from the data for $`0.01t0.4`$ are shown in Fig. 22, along with the error bars on the QMC data. As can be seen from the figure, the only significant deviation of the fit function from the QMC data in this $`t`$ range is at the lowest temperature $`t=0.01`$ for each of the four $`\alpha `$

values. Because the fit deviations at this temperature remain upon increasing the length of the simulated chain from $`L=400`$ to $`L=800`$, these fit deviations are most likely due to inaccuracies in the fit function, as expected at this lowest fitted temperature.

For compounds showing spin-Peierls or other types of second-order spin dimerization transitions, it is more appropriate to scale $`\chi `$ by $`1/J`$ and $`T`$ by $`J`$, where $`J`$ is the average of $`J_1`$ and $`J_2`$, in which case the appropriate alternation parameter is $`\delta `$ rather than $`\alpha `$. It is straightforward to convert our $`\chi ^{}(\alpha ,t)`$ fit function to the form $`\overline{\chi ^{}}(\delta ,\overline{t})`$, where $`\overline{t}k_\mathrm{B}T/J`$, using Eq. (41). We have done this and plot the $`\overline{\chi ^{}}(\delta ,\overline{t})`$ fit function versus temperature for a series of $`\delta `$ values in Fig. 23(a). An appealing monotonic progression of $`\overline{\chi ^{}}(\delta ,\overline{t})`$ with increasing $`\delta `$ is seen in Fig. 23(a); an expanded plot at lower temperatures is shown in Fig. 23(b). This formulation of the fit function allows accurate estimates to be made of the temperature-dependent spin gap in compounds exhibiting spin-dimerization transitions, provided that the nearest neighbor $`S=1/2`$ AF alternating-exchange Heisenberg model is appropriate to them. An illustration of the procedure and the results to be gained will be given later when we model the $`\chi (T)`$ data for NaV2O5.

## VI Spin Gap from TMRG $`𝝌^{\mathbf{}}\mathbf{(}𝜶\mathbf{,}𝒕\mathbf{)}`$

According to Eq. (23), if highly precise $`\chi ^{}(t)`$ data in the low-$`t`$ limit are available, the spin gap $`\mathrm{\Delta }^{}`$ can in principle be computed directly from the derivative of these data with respect to inverse temperature. However, in general the maximum temperature of the low-$`t`$ limit region is ill defined since it depends on how precise and accurate the data are and the accuracy to which $`\mathrm{\Delta }^{}`$ is to be determined. Therefore, in practice one could define a temperature-dependent effective spin gap $`\mathrm{\Delta }_{\mathrm{eff}}^{}`$ from Eq. (23) as

$$\mathrm{\Delta }_{\mathrm{eff}}^{}(t)=\frac{\mathrm{ln}(\chi ^{}\sqrt{t})}{(1/t)},$$
(156)
and then try to ascertain $`\mathrm{\Delta }^{}`$ from the extrapolated zero-temperature limit $`\mathrm{\Delta }^{}=lim_{t0}\mathrm{\Delta }_{\mathrm{eff}}^{}(t)`$. Using Eq. (22) would be less desirable and precise because a fit of this type typically averages $`\mathrm{\Delta }_{\mathrm{eff}}^{}(t)`$ over a rather large temperature range.

An overview of $`\mathrm{\Delta }_{\mathrm{eff}}^{}(\alpha ,t)`$ determined from our TMRG $`\chi ^{}(\alpha ,t)`$ data for $`0.8\alpha 0.995`$ using Eq. (156) is shown in Fig. 24(a). At the lowest temperatures, and for $`\alpha `$ not too close to 1, the $`\mathrm{\Delta }_{\mathrm{eff}}^{}(\alpha ,t)`$ data do approach a constant value with decreasing $`t`$, confirming the applicability of Eqs. (16) and prior assumptions and hence Eqs. (1) and (156) to the alternating chain, and the approximate values of $`\mathrm{\Delta }^{}(\alpha )`$ can be estimated from the figure. Closer inspection reveals that $`\mathrm{\Delta }_{\mathrm{eff}}^{}(\alpha ,t)`$ shows a weak maximum before decreasing by $`\frac{1}{2}`$% to $`\mathrm{\Delta }^{}`$ as $`t0`$, as illustrated in Fig. 24(b) for $`\alpha =0.8`$. For this among other reasons, we will not use Eq. (156) to extract the spin gaps from our TMRG $`\chi ^{}(\alpha ,t)`$ data. On the other hand, we need to know whether such behavior is expected, since it could conceivably arise from systematic errors in the TMRG calculations. Therefore, in the next section we study the $`\mathrm{\Delta }_{\mathrm{eff}}^{}(\alpha ,t)`$ expected at low temperatures for the alternating-exchange chain. As part of this study, we formulate and discuss the fit function which we will use in Sec. VI B to extract $`\mathrm{\Delta }^{}(\alpha )`$ from our TMRG $`\chi ^{}(\alpha ,t)`$ data at low temperatures.

### A Effective spin gap $`𝚫_{\mathrm{𝐞𝐟𝐟}}^{\mathbf{}}\mathbf{(}𝚫^{\mathbf{}}\mathbf{,}𝒕\mathbf{)}`$ for the alternating-exchange chain

From our definition of $`\mathrm{\Delta }_{\mathrm{eff}}^{}`$ in Eq. (156), a discussion of how this quantity varies with $`t`$ at low $`t`$ requires an independent estimate of $`\chi ^{}(\alpha ,t)`$ for the alternating-exchange

chain, which must include at least the leading order correction to the low-$`t`$ limit in Eqs. (16). As a first attempt, we used the general expression for $`\chi ^{}(t)`$ in Eqs. (II B), which requires as input the one-magnon dispersion relation $`\epsilon (k)`$ for the alternating chain. For this we used the explicit $`\epsilon (\mathrm{\Delta }^{},k)`$ for $`0\alpha 1`$ in Eqs. (3) that we presented and discussed previously. The resultant $`\chi ^{}(t)`$ is plotted for eleven $`\mathrm{\Delta }^{}`$ values in Fig. 25(a), where the results are designated by $`\chi ^{(1)}`$ in the figure. Although $`\chi ^{(1)}(\mathrm{\Delta }^{},t)`$ is exact in both the low- and high-$`t`$ limits, the results are only qualitatively correct at intermediate temperatures, as can be seen by comparing Fig. 25(a) with the QMC and TMRG data and fit in Fig. 20. Troyer, Tsunetsugu and Würtz obtained a very good fit of $`\chi ^{(1)}(\mathrm{\Delta }^{},t)`$ to QMC $`\chi ^{}(t)`$ simulation data over a large temperature range for the $`S=1/2`$ two-leg Heisenberg ladder with spatially isotropic exchange; however, they assumed a $`\epsilon (\mathrm{\Delta }^{},k)`$ in the fit function which was later found to be inaccurate over much of the Brillouin zone.

We formulated an approximation \[designated as $`\chi ^{(2)}`$\] which is more accurate in the low-temperature range, and which we will use in the next section as a fit function to fit our TMRG $`\chi ^{}(t)`$ data at low $`t`$ to extract $`\mathrm{\Delta }^{}(\alpha )`$. The function $`\chi ^{(2)}`$ was obtained by summing the susceptibilities of isolated dimers with a distribution of singlet-trip-

let energy gaps given by our one-parameter dispersion relation $`\epsilon (\mathrm{\Delta }^{},k)`$ for $`0\mathrm{\Delta }^{}1`$ in Eqs. (3), which takes into account the interdimer interactions. Thus from Eq. (14) we simply obtain

$$\chi ^{(2)}(\mathrm{\Delta }^{},t)=\frac{1}{\pi t}_0^\pi \frac{dk}{3+\mathrm{e}^{\epsilon (\mathrm{\Delta }^{},k)/t}}.$$
(157)
Note that we make no assumptions here about the form of the function $`\mathrm{\Delta }^{}(\alpha )`$, since only $`\mathrm{\Delta }^{}`$ appears in the expression. This $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ is exact in both the low- and high-$`t`$ limits, as is $`\chi ^{(1)}(\mathrm{\Delta }^{},t)`$, and both reproduce $`\chi ^{}(t)`$ for the isolated dimer ($`\mathrm{\Delta }^{}=1`$) exactly, but $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ is more accurate at intermediate temperatures for $`\alpha 1`$ as shown in Fig. 25(b). In addition, by comparing $`\chi ^{(1)}(\mathrm{\Delta }^{},t)`$ and $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ with the TMRG $`\chi ^{}(\alpha ,t)`$ calculations at low $`t`$, we found that the low-$`t`$ corrections to the low-$`t`$ limit in Eqs. (16) are much more accurately given by $`\chi ^{(2)}(t)`$ than by $`\chi ^{(1)}(t)`$. We will therefore not discuss $`\chi ^{(1)}(t)`$ further here.

At low temperatures, the approximation $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ is expected to accurately describe the leading-order $`t`$ corrections to the low-$`t`$ limit only as long as the average number of magnons $`n_\mathrm{m}`$ occupying a state near the minimum in the one-magnon band is much less than unity. Using the expression for the boson occupation number

for this case,

$$n_\mathrm{m}=\frac{1}{\mathrm{e}^{\mathrm{\Delta }^{}/t}1},$$
(158)
yields $`t/\mathrm{\Delta }^{}=0.22`$ and 0.42 for $`n_\mathrm{m}=0.01`$ and 0.1, respectively. Thus, when fitting our low-$`t`$ TMRG $`\chi ^{}(\alpha ,t)`$ data by the fit function $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ in the following section, we expect $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ to be sufficiently accurate only for $`t/\mathrm{\Delta }^{}0.4`$. For this reason, our fits will be limited to this maximum scaled temperature.

We have computed $`\mathrm{\Delta }_{\mathrm{eff}}(\mathrm{\Delta }^{},T)/\mathrm{\Delta }`$ from $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ in Eq. (157), using the definition in Eq. (156), and plot the results vs $`k_\mathrm{B}T/\mathrm{\Delta }`$ in Fig. 26. For the dimer ($`\mathrm{\Delta }^{}=1`$), one finds analytically that $`\mathrm{\Delta }_{\mathrm{eff}}(T)/\mathrm{\Delta }=12k_\mathrm{B}T/\mathrm{\Delta }`$ to lowest order in $`T`$. On the other hand, for $`0<\mathrm{\Delta }^{}<1`$, the initial dependence is positive and quadratic in $`T`$, and a maximum is seen in $`\mathrm{\Delta }_{\mathrm{eff}}(T)/\mathrm{\Delta }`$, which for $`\mathrm{\Delta }^{}0.4`$ occurs at $`t/\mathrm{\Delta }^{}k_\mathrm{B}T/\mathrm{\Delta }0.14`$ with a height of $`0.5`$%. This height is quantitatively consistent with the data in Fig. 24(b) derived from the TMRG $`\chi ^{}(t)`$ for $`\alpha =0.8`$. Thus the weak maximum seen in that figure is not a spurious effect.

### B Fits to the low-$`𝒕`$ TMRG $`𝝌^{\mathbf{}}\mathbf{(}𝜶\mathbf{,}𝒕\mathbf{)}`$ data

We were tempted to fit $`\mathrm{\Delta }_{\mathrm{eff}}^{}(\alpha ,t)`$ derived from the low-$`t`$ TMRG $`\chi ^{}(\alpha ,t)`$ data, as discussed above, to obtain the spin gaps $`\mathrm{\Delta }^{}(\alpha ,t)`$. However, this procedure would have weighted the $`\chi ^{}(\alpha ,t)`$ data in an ill advised way. We therefore decided to do conventional fits of the low-$`t`$ $`\chi ^{}(\alpha ,t)`$ data by the fit function $`\chi ^{(2)}(\mathrm{\Delta }^{},t)`$ in Eq. (157). For a given $`\alpha `$, this is a one-parameter ($`\mathrm{\Delta }^{}`$) fit function and the fits are therefore stringent tests of both the appropriateness of the fit function and the precision and accuracy of the data. Because the temperature dependence of the accuracy of the calculations is unknown except for the uniform chain data (see Fig. 19), we assumed that all data for a given $`\alpha `$ in a given fitted temperature range have the same accuracy. Thus in the nonlinear least-squares fits for each $`\alpha `$ we minimized the square of the relative rms deviation of the fit from the data

$$\sigma _{\mathrm{rms}}^2=\frac{1}{N_\mathrm{p}}\underset{i=1}{\overset{N_\mathrm{p}}{}}\frac{[\chi ^{(2)}(t_i)\chi ^{}(t_i)]^2}{[\chi ^{}(t_i)]^2},$$
(159)
where $`N_\mathrm{p}`$ is the number of data points fitted, which was usually between 250 and 1500.

Due to the presence of the spin gap $`\mathrm{\Delta }^{}`$ in the exponential of the fit function, $`\sigma _{\mathrm{rms}}`$ is extremely sensitive to the precise value of $`\mathrm{\Delta }^{}`$ when low-$`t`$ fits are carried out. For example, close to the optimum $`\mathrm{\Delta }^{}`$ fit value, a change in $`\mathrm{\Delta }^{}`$ by only 0.0001 ($`0.1\%`$) can change $`\sigma _{\mathrm{rms}}`$ by up to $`300`$ %. Thus a few percent accuracy in the $`\chi ^{}(t)`$ data at low $`t`$ is sufficient to allow $`\mathrm{\Delta }^{}`$ to be determined for a given fit to a precision better than 0.0001. For a given $`\alpha `$, the obtained $`\mathrm{\Delta }^{}`$ was found to be insensitive, typically to within $`0.000`$2, to the $`t`$ range of the fit, as long as the maximum fitted temperature satisfied $`t/\mathrm{\Delta }^{}0.4`$, consistent with the above discussion of the boson occupation number. This lack of sensitivity of the value of the fitted $`\mathrm{\Delta }^{}`$ to the precise fit range demonstrated that the fit function $`\chi ^{(2)}(t)`$ is an appropriate one. Depending on the $`\alpha `$ value and the $`t`$ range of the fit, $`\sigma _{\mathrm{rms}}`$ was typically between 0.1 % and several percent.

The $`\mathrm{\Delta }^{}(\alpha )`$ values obtained from the fits are listed in Table III, together with the estimated accuracies in parentheses. Note that a quoted accuracy is associated with variations in $`\mathrm{\Delta }^{}`$ in fits to a specific set of data for a given $`\alpha `$ over various temperature ranges as discussed above, and does not include possible systematic errors due to, e.g., the finite fixed number of states kept in the TMRG calculations. Also included in Table III are literature data which will be compared with the present results in the next section. Log-log plots of the low-$`t`$ data and fits are shown in Fig. 27, where on the scale of this figure, for most $`\alpha `$ values the data and fit are identical (cannot be distinguished) within the fitted temperature range. Extrapolations of the fits to higher and lower temperatures are also shown for comparison with the data.

## VII Comparisons of the Calculations with Previous Work

### A Spin gap

Our $`\overline{\mathrm{\Delta }^{}}(\delta )`$ spin gap data determined by fitting our TMRG $`\chi ^{}(t)`$ data by Eq. (157) are plotted in Fig. 28(a),

along with the results of previous workers listed in Table III. The solid curve is the function $`\overline{\mathrm{\Delta }^{}}=2\delta ^{3/4}`$ in Eq. (48) proposed by Barnes, Riera, and Tennant (BRT). The overall behavior of the data in Fig. 28(a) is well described by this function, but significant deviations of the data from the curve occur as illustrated in the expanded plot for $`\delta 0.1`$ in Fig. 28(b). The error bars are included with each plotted data point in Fig. 28(b), except for the data of Ref.  which were not available, but they are all small and not clearly seen. Our values for $`\delta 0.1`$ are significantly smaller than those of Uhrig et al., where the differences are far outside the combined limits of error, and are larger than those of Augier et al.

As will be seen explicitly in Sec. VIII C below, our $`\overline{\chi ^{}}(\delta ,\overline{t})`$ fit function allows $`\delta (T)`$ to be determined for real materials by using the fit function to model experimental $`\chi (T)`$ data. However, if one would like to determine the spin gap $`\mathrm{\Delta }(T)`$ from the derived $`\delta (T)`$, an expression is needed for $`\mathrm{\Delta }(\delta )`$ over the entire range $`0\delta 1`$ in order to be generally useful and applicable. At present, the only extant expression is that of BRT in Eqs. (II C 2). As seen in Fig. 28 and in Table IV below, this expression is only an approximation that fits neither BRTs’ $`\mathrm{\Delta }^{}(\alpha )`$ data for $`0.1\alpha 0.9`$ nor our TMRG spin gap data for $`0.8\alpha 0.995`$ to within the respective error bars. To formulate a more flexible expression, we modify BRTs’ formula to read

$$\overline{\mathrm{\Delta }^{}}(\delta )\frac{\mathrm{\Delta }(\delta )}{J}=2\delta ^{y(\delta )},$$
(161)
so the $`\delta `$-dependent power $`y`$ is

$$y(\delta )=\frac{\mathrm{ln}[\mathrm{\Delta }(\delta )/2J]}{\mathrm{ln}\delta }.$$
(162)
The numerical prefactor “2” in Eq. (161) must be retained in order to reproduce the exact $`\overline{\mathrm{\Delta }^{}}(\delta =1)=2`$. Shown in Fig. 29(a) is a semilog plot of $`y`$ versus $`\delta `$ for the same numerical data as in Fig. 28. This plot \[and Fig. 29(b) below\] explicitly shows, from BRTs’ data, that the exponent deviates significantly from the value 3/4 even for $`\delta 1`$. The plot also clearly differentiates the various numerical data for small $`\delta `$ by the different groups, and shows that one of our two data points from the TMRG for $`\delta =0.0025`$ (the one derived from $`m=256`$ data at high $`t`$) is not in agreement with the trend of the remainder of our data. This data point will not be included in the plot and fit to be discussed in the next paragraph.

Our $`y(\delta )`$ data at small $`\delta `$ are in agreement with both the magnitude and trend of BRTs’ data at larger $`\delta `$. The

$`y(\delta )`$ for these two data sets from Fig. 29(a), with the exception of one of our two data points for $`\delta =0.0025`$ just noted above, are plotted together on an expanded vertical scale in Fig. 29(b) where a rather smooth behavior of $`y(\delta )`$ is seen over the combined range of the two calculations $`0.0025\delta 1`$. With the behavior in Fig. 29(b) in mind, we formulated a five-parameter fit function for these two combined $`y(\delta )`$ data sets that yields the correct limits $`\overline{\mathrm{\Delta }^{}}(\delta 0)=0`$ and $`\overline{\mathrm{\Delta }^{}}(\delta 1)=2`$, with the property $`lim_{\delta 0}y(\delta )=`$ const, given by

$`y(\delta )=y(1)`$ $`+`$ $`n_1\mathrm{tanh}\left[{\displaystyle \frac{\mathrm{ln}\delta }{m_1}}\mathrm{ln}\left({\displaystyle \frac{\mathrm{ln}\delta }{m_2}}\right)\right]`$ (163)
$`+`$ $`n_2\mathrm{tanh}^2\left[{\displaystyle \frac{\mathrm{ln}\delta }{m_1}}\mathrm{ln}\left({\displaystyle \frac{\mathrm{ln}\delta }{m_2}}\right)\right].`$ (164)
An unweighted fit of this expression to all the data in Fig. 29(b) yielded the parameters

$`y(1)=0.74922,n_1=0.00776,n_2=0.00685,`$

$$m_1=3.3297,m_2=2.2114,$$
(165)
so that

$$\underset{\delta 0}{lim}y(\delta )=y(1)n_1+n_2=0.7346.$$
(166)
The fit is plotted as the solid curve in Fig. 29(b). As can be seen from the figure, our data are fitted to within our error bars. In addition, when the $`y(\delta )`$ fit function in Eqs. (164) and (165) is inserted into Eq. (161), the predicted values of $`\overline{\mathrm{\Delta }^{}}(\delta )`$ are in agreement with each of the values of BRT at larger $`\delta `$ to within 0.0001, which is sufficient for modeling experimental data. The $`\delta =0`$ limit of $`y(\delta )`$ in Eq. (166) is in agreement with the theoretical effective value $`y(0)=0.72(3)`$, which was obtained without the log correction term by Singh and Weihong from an eleventh-order dimer series expansion of the triplet dispersion relation. We will use Eqs. (27) to compute $`\mathrm{\Delta }(T)`$ from the experimentally derived $`\delta (T)`$ for $`\mathrm{NaV}_2\mathrm{O}_5`$ in Sec. VIII C below.

In order to test the critical behavior prediction $`\overline{\mathrm{\Delta }^{}}=A\delta ^{2/3}/|\mathrm{ln}\delta |^{1/2}`$ in Eq. (49), which need only hold in the asymptotic critical regime $`\delta 0`$ in contrast to the fit function for $`0\delta 1`$ in Eqs. (27), in Fig. 30 is plotted $`\mathrm{\Delta }/J`$ vs $`\delta ^{2/3}/|\mathrm{ln}\delta |^{1/2}`$ in the region $`\delta 0.06`$ for the same data and symbols as in Fig. 28. A proportionality appears to be developing in our data for $`\delta 0.005`$, as shown by the straight line with slope $`A=3.3`$ passing through the origin of the figure, suggesting that the asymptotic critical regime begins with decreasing $`\delta `$ below $`\delta 0.005(\alpha 0.99)`$. High-accuracy $`\overline{\mathrm{\Delta }^{}}(\delta )`$ data for $`\delta 0.001`$ are needed to test this conjecture. From Fig. 30, the slope 3.3 of the line drawn is evidently a lower limit of the prefactor $`A`$ within the actual asymptotic critical regime.

### B Numerical $`𝝌^{\mathbf{}}\mathbf{(}𝜶\mathbf{,}𝒕\mathbf{)}`$ results

Barnes and Riera previously carried out exact diagonalizations of Hamiltonian (II C) for $`S=1/2`$ alternating chains of length up to 16 spins using the Lanczos technique. Their computed $`\chi ^{}(t)`$ values for $`\alpha =0.2,`$ 0.4, 0.6, 0.7, and 0.8 were extrapolated to the bulk limit and the results are shown as the symbols in Fig. 31(a). Our fit function as in Fig. 20 for the same $`\alpha `$ values is plotted as the solid curves in Fig. 31(a), which are seen to be in good overall agreement with the calculations of Barnes and Riera. The deviations of the data of Barnes and Riera from our fit function are plotted vs temperature in Fig. 31(b). The average deviation of their data from our fit function is very small for each data set: $`0.41,+0.33,0.40,0.26`$, and $`+0.79\times 10^4`$ for $`\alpha =0.2`$, 0.4, 0.6, 0.7, and 0.8, respectively. The absolute rms deviations $`\sigma _{\mathrm{rms}}`$ of their data from our fit function for $`\alpha =0.2`$, 0.4, 0.6, 0.7, and 0.8 are (in units of $`10^4`$) 1.73, 1.43, 0.73, 0.78, and 3.76, respectively. We conclude that their data are in good quantitative agreement with our data and fit, with the exception of their data point for $`\alpha =0.8`$ at their lowest temperature $`t=0.05`$.

### C Bulaevskii Theory

Bulaevskii calculated $`\chi ^{}(t)`$ analytically in the Hartree-Fock approximation. He first obtained an integral equation for the magnon spectrum $`E(k)`$:

$`\epsilon (k){\displaystyle \frac{E(k)}{J_1}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\sqrt{1+\alpha ^22\alpha \mathrm{cos}k}`$ (167)
$`+`$ $`{\displaystyle \frac{C_1+\alpha C_2(\alpha C_1+C_2)\mathrm{cos}k}{\sqrt{1+\alpha ^22\alpha \mathrm{cos}k}}}],`$ (168)
$`C_1(t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^\pi }dk{\displaystyle \frac{1\alpha \mathrm{cos}k}{\sqrt{1+\alpha ^22\alpha \mathrm{cos}k}}}\mathrm{tanh}{\displaystyle \frac{\epsilon (k)}{2t}},`$ (169)
$`C_2(t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^\pi }dk{\displaystyle \frac{\alpha ^2\alpha \mathrm{cos}k}{\sqrt{1+\alpha ^22\alpha \mathrm{cos}k}}}\mathrm{tanh}{\displaystyle \frac{\epsilon (k)}{2t}},`$ (171)
where $`k`$ is measured in units of $`2\pi /d`$. $`d1`$ is the lattice repeat distance along the chain, which is twice the average distance between spins. He then expressed $`\chi ^{}(t)`$ in terms of $`\epsilon (k)`$:

$`\chi ^{}(t)`$ $``$ $`{\displaystyle \frac{F(t)}{2+(1+\alpha )F(t)}},`$ (172)
$`F(t)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi t}}{\displaystyle _0^\pi }{\displaystyle \frac{\mathrm{d}k}{\mathrm{cosh}^2[\epsilon (k)/(2t)]}}.`$ (174)
At $`t=0`$ and $`\alpha 0`$, from Eqs. (LABEL:EqEpsK2) we obtain

$`C_1(\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\{(1+\alpha )\mathrm{E}\left[{\displaystyle \frac{4\alpha }{(1+\alpha )^2}}\right]`$ (176)
$`+(1\alpha )\mathrm{K}\left[{\displaystyle \frac{4\alpha }{(1+\alpha )^2}}\right]\},`$
$`C_2(\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\{(1+\alpha )\mathrm{E}\left[{\displaystyle \frac{4\alpha }{(1+\alpha )^2}}\right]`$ (179)
$`(1\alpha )\mathrm{K}\left[{\displaystyle \frac{4\alpha }{(1+\alpha )^2}}\right]\},`$
where K$`(y)`$ and E$`(y)`$ are, respectively, the complete elliptic integrals of the first and second kinds. The dispersion relations versus $`\alpha `$ at $`t=0`$ are obtained by inserting Eqs. (LABEL:EqEpsK3) into (168) and a selection of results is shown in Fig. 32. From Eqs. (168) and (LABEL:EqEpsK3), at $`t=0`$ the spin-gap $`\mathrm{\Delta }_{k=0}(\alpha )`$ at $`k=0`$ is given by

$$\mathrm{\Delta }_{k=0}(\alpha )=\frac{1\alpha }{2}\left\{1+\frac{2(1\alpha )}{\pi }\mathrm{K}\left[\frac{4\alpha }{(1+\alpha )^2}\right]\right\}.$$
(180)
This expression gives the actual spin-gap for $`0<\alpha 0.79`$. However, for $`0.79\alpha <1`$, the minimum in the dispersion relation does not occur at $`k=0`$, as illustrated in an expanded plot of $`E(k)`$ for $`\alpha =0.9`$ in the inset to Fig. 32. The wave vector $`k_\mathrm{G}`$ at which the minimum spin gap $`\mathrm{\Delta }_\mathrm{B}`$ occurs is plotted versus $`\alpha `$ in Fig. 33(b). The $`\mathrm{\Delta }_\mathrm{B}`$ from Bulaevskii’s theory at $`t=0`$ is plotted versus $`\alpha `$ as the solid curve in Fig. 33(a), and a few representative values are given in Table IV. The predictions of Bulaevskii’s theory are in very good agreement with those of Barnes, Riera, and Tennant for $`\alpha 0.3`$, but the agreement becomes progressively worse as $`\alpha `$ increases further. From Eqs. (168) and (LABEL:EqEpsK2), $`E(k)`$ is temperature dependent. In addition, in the range $`0.79\alpha <1`$ for which $`k_\mathrm{G}0`$ at $`t=0`$, we find that $`k_\mathrm{G}`$ depends on $`t`$, as shown in Fig. 34. From Fig. 34, $`k_\mathrm{G}0`$ at $`t0.083`$, 0.122, 0.131, 0.125, and 0.111 for $`\alpha =0.8`$, 0.85, 0.9, 0.95, and 0.99, respectively.

We computed $`\chi ^{}(t)`$ by inserting $`\epsilon (k)`$ in Eq. (168) into Eqs. (LABEL:EqEpsK2), numerically solving the latter two simultaneous equations for $`C_1`$ and $`C_2`$ at each $`t`$, and then inserting the resulting $`\epsilon (k)`$ into Eqs. (LABEL:EqChiBul). The progression of $`\chi ^{}(t)`$ with increasing $`\alpha `$ from 0.001 to 0.99 is shown in Fig. 35. As noted by Bulaevskii, the values of $`\chi ^{}`$ at the maxima are too large and the temperatures at which these occur are too small by $`5`$–10% (compare Fig. 35 with Fig. 20).

At low temperatures $`0.033t1/4`$, Bulaevskii fitted $`\chi ^{}(\alpha ,t)`$, calculated from Eqs. (LABEL:EqChiBul), by the two-parameter form

$$\chi ^{}(\alpha ,t)=\frac{a(\alpha )}{t}\mathrm{e}^{\mathrm{\Delta }(\alpha )/(J_1t)},$$
(181)
and obtained values of $`a_\mathrm{B}`$ and $`\mathrm{\Delta }_{\mathrm{B},\mathrm{Fit}}/J_1`$ for $`0\alpha 0.9`$ which are reproduced in Table IV; $`\mathrm{\Delta }_{\mathrm{B},\mathrm{Fit}}(\alpha )/J_1`$ is plotted as the open squares in Fig. 33(a). Note that the temperature exponent in the prefactor to the exponential is $`\gamma =1`$, contrary to the $`\gamma =1/2`$ in Eq. (19) which is expected in the low-$`t`$ limit for any 1D $`S=1/2`$ Heisenberg spin system with a spin gap (and with a nondegenerate one-magnon band with a parabolic minimum). We have confirmed that over the temperature range fitted by Bulaevskii, one indeed obtains $`\gamma 1`$ for the best fit of Eq. (181) to numerical calculations of $`\chi ^{}(\alpha ,t)`$. We infer that the discrepancy between Bulaevskii’s $`\gamma =1`$ and the expected $`\gamma =1/2`$ arises because the fits were not carried out completely within the low-$`t`$ limit. This issue is discussed in more detail below.

Equation (181) together with Bulaevskii’s table of {$`a_\mathrm{B}(\alpha ),\mathrm{\Delta }_{\mathrm{B},\mathrm{Fit}}(\alpha )/J_1`$} values were subsequently used extensively in the analysis of experimental $`\chi (T)`$ data for compounds exhibiting spin-Peierls transitions to determine the alternation parameter $`\alpha `$ at low temperatures $`TT_\mathrm{c}`$ where the experimental spin-gap is nearly independent of $`T`$. However, from Table IV and Fig. 33(a), the $`\mathrm{\Delta }_{\mathrm{B},\mathrm{Fit}}(\alpha )`$ values of Bulaevskii are in generally poor agreement with the actual spin gaps $`\mathrm{\Delta }_\mathrm{B}(\alpha )`$ of his theory and with those \[$`\mathrm{\Delta }_{\mathrm{BRT}}(\alpha )`$\] calculated for the same $`\alpha `$ values by Barnes, Riera, and Tennant. Therefore, one should consider the $`\mathrm{\Delta }_{\mathrm{B},\mathrm{Fit}}`$ parameters as fitting parameters only, with no direct relation to the actual spin gap.

According to Eq. (14) for $`\chi ^{}(t)`$ of the isolated dimer which is a zero-dimensional spin system, the form of $`\chi ^{}(t)`$ in Eq. (181) with $`\gamma =1`$ is correct for $`\alpha =0`$ and $`t0`$. On the other hand, for one-dimensional spin systems such as the two-leg spin ladder (and the alternating-exchange chain) at temperatures $`k_\mathrm{B}T\mathrm{\Delta }`$ and $`k_\mathrm{B}T`$ one-magnon bandwidth, Eqs. (16) apply, with $`\gamma =1/2`$, assuming that the triplet one-magnon dispersion relation $`E(k)`$ is parabolic at the minimum. In this case one expects $`\gamma =1/2`$ at sufficiently low $`t`$ for any finite $`\alpha `$. Thus, in the temperature region of validity of Eq. (18), a plot of the left-hand-side of Eq. (21) vs ln$`t`$ should give a straight line with slope $`\gamma `$. Shown in Fig. 36 are such plots, obtained using our $`\chi ^{}(\alpha ,t)`$ calculated from Bulaevskii’s theory as described above, for $`\alpha =0.001`$ to 0.99. For $`\alpha =0.001`$, a crossover is clearly evident from $`\gamma =1`$ to $`\gamma =1/2`$ with decreasing $`t`$. The other curves also exhibit signs of a crossover, with $`\gamma 1/2`$ at the lowest temperatures, with the exception of the curve for $`\alpha =0.8`$. For this $`\alpha `$ value, which is just above the value $`\alpha 0.79`$ at which $`k_\mathrm{G}`$ becomes nonzero at $`t=0`$ \[see Fig. 33(b)\], the $`\gamma `$ at the lowest $`t`$ is intermediate between the values of 1/2 and 1, and the assumption of a parabolic form for $`E(k)`$ at the band minimum is evidently not satisfied (see Fig. 32). In fact, Troyer, Tsunetsugu and Würtz calculated the low-$`t`$ limit of $`\chi ^{}(t)`$ for 1D systems with general dispersion relation $`\epsilon (k)=\mathrm{\Delta }^{}+c^{}|ka|^n`$, where $`k`$ is the deviation of the wave vector from that at the band minimum. They found the same form $`\chi ^{}(t)=(A_n/t^\gamma )\mathrm{exp}(\mathrm{\Delta }^{}/t)`$ as for the parabolic case $`n=2`$, but where $`\gamma =1(1/n)`$. Thus, e.g., $`\gamma =2/3`$ and 3/4 for $`n=3`$ and 4, respectively. This range of $`\gamma `$ is consistent with the slope of the data at the lowest temperatures for $`\alpha =0.8`$ in Fig. 36.

The predictions of Bulaevskii’s theory for $`\chi ^{}(t)`$ from Fig. 35 are compared with our $`\chi ^{}(\alpha ,t)`$ fit function (solid

curves) for $`\alpha =0.2`$, 0.4, 0.6, 0.8 and 0.99 (as in Fig. 20) in Fig. 37, where the Bulaevskii prediction for each of these $`\alpha `$ values is shown as the corresponding dashed curve. The disagreement between the two calculations becomes progressively more severe as temperature decreases and as the uniform chain limit is approached with increasing $`\alpha `$. Therefore, the accuracies of the $`\alpha `$ and $`J_1`$ values previously extracted from experimental data at low $`T`$ for compounds with $`\alpha 1`$ using Bulaevskii’s theory are unclear. Our $`\chi ^{}(\alpha ,t)`$ fit function now provides a much more accurate and reliable means of extracting exchange constants and spin gaps from experimental $`\chi (T)`$ data.

## VIII Magnetic Susceptibility of NaV2O5

Crystals of $`\mathrm{Na}_{0.996(3)}\mathrm{V}_2\mathrm{O}_{5.00(6)}`$ were grown at the Max-Planck-Institut für Festkörperforschung, Stuttgart, in a Pt crucible in flowing Ar atmosphere by a self-flux method from a 5:1:1 mixture of $`\mathrm{NaVO}_3`$, $`\mathrm{V}_2\mathrm{O}_3`$ and $`\mathrm{V}_2\mathrm{O}_5`$. The flux was dissolved by boiling the solidified melt in distilled water. X-ray powder diffraction patterns collected with a STOE diffractometer yielded the lattice parameters $`a=11.3187(8)`$ Å, $`b=3.6111(3)`$ Å, and $`c=4.8007(5)`$ Å. Chemical analyses on two independent representative samples of the batch were performed with a standard AAS analysis technique for V and Pt and ICP emission spectroscopy for the Na content. The oxygen content was determined by measuring with IR spectroscopy the amount of CO generated when the sample is fused in a graphite crucible at 2700 $`{}_{}{}^{\mathrm{o}}\mathrm{C}`$ in vacuo. Platinum impurities above the level of sensitivity of the analysis (500 ppm with respect to V) could not be detected.

At Ames Laboratory, single crystals of $`\mathrm{NaV}_2\mathrm{O}_5`$ were grown out of the ternary melt. Powders of $`\mathrm{V}_2\mathrm{O}_5`$ and V2O3 were prepared by oxidizing and reducing NH4VO3 at 600$`^{}`$C and 900 C, respectively. The resulting $`\mathrm{V}_2\mathrm{O}_5`$ is reacted with Na2CO3 at 550 C yeilding NaVO3. About 10 grams of NaVO3, $`\mathrm{V}_2\mathrm{O}_5`$, and V2O3 in the molar ratio 32:1:1 were placed in a Pt crucible and sealed in an evacuated quartz tube. The melt was then slowly cooled from 800 to 660 C over 50 h and the remaining liquid was decanted. Small amounts of solidified melt remaining on the crystals were dissolved with hot water. Typical dimensions of the ribbon-shaped crystals grown in this manner are $`0.5\times 1.5\times 11`$ mm3 with the $`c`$ axis perpendicular to the plane of the ribbon, the $`b`$ axis along the length of the ribbon and the $`a`$ axis along the width of the ribbon, with lattice parameters $`a11.303`$ Å, $`b3.611`$ Å, and $`c4.752`$ Å. The crystal denoted as AL1 has a mass of 8.2 mg and approximate planar dimensions $`1.5\times 2.5`$ mm2.

The magnetic susceptibility $`\chi (T)M(T)/H`$ of the crystals was measured using Quantum Design SQUID magnetometers at Stuttgart and Ames. The measurements on eight crystals of NaV2O5 in Stuttgart were carried out in a field $`H=1`$ T along the V ladder ($`b`$) axis direction in various temperature ranges between 2 and 750 K. Measurements of the anisotropy of $`\chi (T)`$ along the $`a`$, $`b`$, and $`c`$ axis directions were carried out from 2 to 300 K in Ames on crystal AL1 in $`H=2`$ T.

The results for two of the crystals up to 750 K are shown in Fig. 38. The data illustrate the variabilities we have observed between measurements along the same axis on different crystals. Above $`50`$ K, the two data sets are nearly parallel, with the difference between them being $`34\times 10^5`$ cm3/mol; we have no explanation for this difference, and no comments have been made in the literature about such variabilities and/or their origins in $`\chi (T)`$ along the same axis in different crystals that we are aware of. The data from $`T_\mathrm{c}33`$–34 K up to 300 K are in approximate agreement with the single crystal data of Isobe, Kagumi and Ueda taken in this $`T`$ range along the same axis in $`H=5`$ T. A variable Curie-Weiss-like contribution $`\chi ^{\mathrm{CW}}(T)`$ to $`\chi (T)`$ occurs below $`20`$ K which is attributed to paramagnetic defects, impurities, inclusions and/or intergrowths in the crystals. The “Fit” shown in the figure will be discussed later in Sec. VIII B.

The experimental data are analyzed with the general expression

$$\chi (T)=\chi _0+\chi ^{\mathrm{CW}}(T)+\chi ^{\mathrm{spin}}(T),$$
(183)
$$\chi _0=\chi ^{\mathrm{core}}+\chi ^{\mathrm{VV}},$$
(184)
$$\chi ^{\mathrm{CW}}(T)=\frac{C_{\mathrm{imp}}}{T\theta },$$
(185)
$$\chi ^{\mathrm{spin}}(T)=\frac{Ng^2\mu _\mathrm{B}^2}{J}\overline{\chi ^{}}\left(\frac{k_\mathrm{B}T}{J}\right),$$
(186)
where $`\chi _0`$ is the sum of a temperature independent and (nearly) isotropic orbital diamagnetic core contribution and a usually anisotropic and temperature independent orbital paramagnetic Van Vleck contribution. We estimate $`\chi ^{\mathrm{core}}`$ using the values $``$5, $``$7, $``$4, and $`12\times 10^6`$ cm3/mol for Na+1, V+4, V+5, and O-2, respectively, yielding the isotropic value

$$\chi ^{\mathrm{core}}=7.8\times 10^5\frac{\mathrm{cm}^3}{\mathrm{mol}\mathrm{NaV}_2\mathrm{O}_5}.$$
(187)
The second term in Eq. (183) is the above-noted Curie-Weiss impurity and/or defect contribution and the last term is the intrinsic spin susceptibility, each of which may or may not be anisotropic. For a Heisenberg spin system, $`\overline{\chi ^{}}`$ is isotropic, and therefore so is $`\chi ^{\mathrm{spin}}`$ apart from anisotropy in the $`g`$ factor. The impurity Curie-Weiss term $`\chi ^{\mathrm{CW}}(T)`$ can be anisotropic if the impurities are defects or intergrowths in the crystals which have atomic coordination principal axes which are fixed with respect to the crystal axes rather than being randomly oriented. We model our $`\chi (T)`$ data according to Eq. (183) in terms of the $`\overline{\chi ^{}}(t)`$ in Eq. (186), which are (fit functions to) theoretical susceptibility calculations presented in previous sections. Before moving on to do that, we first experimentally examine the anisotropy in $`\chi (T)`$ of $`\mathrm{NaV}_2\mathrm{O}_5`$ and its implications in the next section.

### A Anisotropy of the magnetic susceptibility

The magnetic susceptibilities of $`\mathrm{NaV}_2\mathrm{O}_5`$ crystal AL1 along the $`a`$, $`b`$, and $`c`$ axes are plotted vs temperature in Fig. 39(a), where the $`a`$ and $`c`$ axes are perpendicular to the V chains which run along the $`b`$ axis, and the $`c`$ axis is perpendicular to the trellis layers that the V chain/ladders reside in. The data are similar to the aniso-

tropic $`\chi (T)`$ data reported by Isobe, Kagami and Ueda, although the anisotropies we measure at both room temperature and at low temperatures are somewhat larger than they reported. The anisotropies at low temperatures are seen more clearly if the respective impurity term $`\chi ^{\mathrm{CW}}(T)`$ in Eq. (185) is subtracted from each data set, as shown in Fig. 39(b). The impurity Curie constant $`C_{\mathrm{imp}}`$ and Weiss temperature $`\theta `$ for each direction of the applied field were determined by the requirement that $`\chi (T)`$ become independent of $`T`$ for $`T0`$. The fitted values of $`C_{\mathrm{imp}}`$ were found to be slightly anisotropic and are given in Table VI below. The values of $`C_{\mathrm{imp}}`$ are equivalent to the contribution of only 0.07 mol% of $`S=1/2`$ impurities with $`g=2`$; if the impurity spin is actually greater than 1/2, the concentration of paramagnetic impurities could be much less than this estimate. From a comparison of Figs. 39(a) and 39(b), $`\chi ^{\mathrm{CW}}(T)`$ is seen to make a negligible contribution to the measured $`\chi (T)`$ above $`100`$ K. Since in the presence of a spin gap $`\chi ^{\mathrm{spin}}=0`$ at the lowest temperatures, from Eqs. (38) and Fig. 39(b) we obtain

$`\chi _b^{\mathrm{VV}}`$ $`=`$ $`18.7\times 10^5{\displaystyle \frac{\mathrm{cm}^3}{\mathrm{mol}}},\chi _c^{\mathrm{VV}}=13.3\times 10^5{\displaystyle \frac{\mathrm{cm}^3}{\mathrm{mol}}},`$ (188)
$`\chi _a^{\mathrm{VV}}`$ $`=`$ $`20.0\times 10^5{\displaystyle \frac{\mathrm{cm}^3}{\mathrm{mol}}}(TT_\mathrm{c}).`$ (190)
From Fig. 39(b), the anisotropies of $`\chi (T)`$ are seen to be quite temperature-dependent upon heating through $`T_\mathrm{c}=33.4`$ K. These results are surprising, because $`\chi ^{\mathrm{spin}}`$ is expected to be isotropic (apart from the small anisotropy due to the anisotropic $`g`$ factor), with $`\chi ^{\mathrm{spin}}(T0)=0`$ because of the spin gap, and the anisotropic $`\chi ^{\mathrm{VV}}`$ values are expected to be temperature independent for our $`S=1/2`$ system over the temperature range of our measurements. Thus one expects the difference $`\chi _\alpha (T)\chi _\beta (T)(\alpha ,\beta =a,b,c)`$ to be nearly independent of temperature compared with the magnitude of either, where a subscript refers to the crystallographic axis along which the magnetic field is applied.

To be more quantitative, we define the anisotropy in the intrinsic susceptibility as

$$\mathrm{\Delta }\chi _{\alpha \beta }(T)[\chi _\alpha \chi _\alpha ^{\mathrm{CW}}](T)[\chi _\beta \chi _\beta ^{\mathrm{CW}}](T),$$
(191)
which eliminates extrinsic anisotropy in the Curie-Weiss impurity contribution from the values calculated from the experimental data. For example, according to this definition, $`\mathrm{\Delta }\chi _{ac}(T)`$ is the difference between the uppermost and lowermost data sets in Fig. 39(b). The three $`\mathrm{\Delta }\chi _{\alpha \beta }(T)`$ anisotropies are plotted in Fig. 40. It seems to us that the only reasonable explanation for the strong temperature-dependent anisotropies in Fig. 40 for two of the three data sets is that one or more of the $`\chi _\alpha ^{\mathrm{VV}}`$ susceptibilities is strongly temperature dependent near $`T_\mathrm{c}`$, contrary to our initial expectations. Such a temperature dependence may be associated with the crystallographic and charge-ordering transitions which occur at or near the same temperature as the spin dimerization transition, as discussed in the Introduction.

One can make rather strong general statements about the magnetic susceptibility anisotropies and their temperature dependences as follows. Defining the Van Vleck susceptibility anisotropy $`\mathrm{\Delta }\chi _{\alpha \beta }^{\mathrm{VV}}=\chi _\alpha ^{\mathrm{VV}}\chi _\beta ^{\mathrm{VV}}`$ and similarly the spin susceptibility anisotropy $`\mathrm{\Delta }\chi _{\alpha \beta }^{\mathrm{spin}}=\chi _\alpha ^{\mathrm{spin}}\chi _\beta ^{\mathrm{spin}}`$, from Eqs. (38) one obtains an expression for the anisotropy $`\mathrm{\Delta }\chi _{\alpha \beta }(T)`$ for a spin system in which the only anisotropy in $`\chi ^{\mathrm{spin}}`$ arises from anisotropy in the $`g`$ factor, given by

$$\mathrm{\Delta }\chi _{\alpha \beta }(T)=\mathrm{\Delta }\chi _{\alpha \beta }^{\mathrm{VV}}+\frac{N\mu _\mathrm{B}^2(g_\alpha ^2g_\beta ^2)}{J}\overline{\chi ^{}}\left(\frac{k_\mathrm{B}T}{J}\right).$$
(192)
The reduced spin susceptibility $`\overline{\chi ^{}}(\overline{t})`$ is necessarily positive, and it is isotropic for a Heisenberg spin system as noted above. Thus, if $`\chi _\alpha ^{\mathrm{VV}}`$ and $`\chi _\beta ^{\mathrm{VV}}`$ and therefore $`\mathrm{\Delta }\chi _{\alpha \beta }^{\mathrm{VV}}`$ are independent of temperature, the slope $`\mathrm{\Delta }\chi _{\alpha \beta }(T)/T`$ must have the same sign as the difference $`g_\alpha ^2g_\beta ^2`$. As discussed in the next subsection, for $`\mathrm{NaV}_2\mathrm{O}_5`$, this difference has been reported to be positive for

$`\alpha \beta =ac`$ and $`bc`$ and near zero for $`\alpha \beta =ab`$, consistent with the slopes in Fig. 40. However, in a simple ionic crystalline electric field model and with a positive spin-orbit coupling parameter for V one would predict that a $`\chi _\alpha ^{\mathrm{VV}}`$ should increase with the negative deviation of $`g_\alpha `$ from the free electron value $`g=2`$. Thus, a particularly visible and puzzling discrepancy is that since $`(2g_a)(2g_b)<(2g_c)`$ according to the reported $`g_\alpha `$ values below, on this basis one strongly expects $`\chi _a^{\mathrm{VV}}\chi _b^{\mathrm{VV}}<\chi _c^{\mathrm{VV}}`$; thus two of the three $`\chi _\alpha ^{\mathrm{VV}}`$ values should be about the same and smaller than the third one. Qualitatively contrary to this expectation, for $`TT_\mathrm{c}`$ we observe in Eq. (190) that $`\chi _a^{\mathrm{VV}}\chi _b^{\mathrm{VV}}>\chi _c^{\mathrm{VV}}`$.

We will not emphasize or further discuss these puzzling discrepancies with expectation with respect to their possible influence on our theoretical modeling of our $`\chi (T)`$ data in Secs. VIII B and VIII C, since at present there is no way to model, e.g., a temperature dependent Van Vleck susceptibility which changes rapidly near $`T_\mathrm{c}`$, but the anisotropic susceptibility results and the above discussion should be kept in mind. In the following two subsections the reported anisotropies in the $`g`$ factor as measured using electron spin resonance (ESR) and in the Van Vleck susceptibility as deduced from nuclear magnetic resonance (NMR) measurements will be discussed, respectively, in light of our anisotropic $`\chi (T)`$ data.

##### Anisotropy in the $`g`$ factor from ESR.

Many ESR measurements have recently been reported for NaV2O5. Each study found a signal with $`g2`$ which was attributed to bulk V species, and the $`g`$ values found in the various studies were the same within the errors, e.g.,

$$g_a=g_b=1.972(2),g_c=1.938(2).$$
(193)
The powder-average value is $`g=\sqrt{(g_a^2+g_b^2+g_c^2)/3}=1.961(2)`$. The $`g`$ values were found to be independent of $`T`$ down to 20 K, which is below $`T_\mathrm{c}`$. From all these measurements, there is no indication that the $`S=1/2`$ Heisenberg Hamiltonian is not appropriate to the spin system in NaV2O5. Unfortunately, given the sensitivity of the ESR technique, we cannot be certain that these ESR results are representative of the bulk spin species in NaV2O5, because no quantitative measurements of the concentration of spin species observed in these measurements were reported. Although the (uncalibrated) ESR intensity versus temperature measurements approximately mirror the bulk susceptibility behavior in most (but not all) of these studies, it is still possible that the signal arises from a minority spin species that is coupled to the bulk spin system. An interesting related issue which has not been discussed in the literature is why the presumed bulk $`S=1/2`$ species in NaV2O5 are observable to low temperatures $`T0.03J/k_\mathrm{B}`$ by ESR, where the AF exchange constant is $`J/k_\mathrm{B}580`$ K (see below), whereas the bulk Cu+2 spins 1/2 in the high-$`T_\mathrm{c}`$ cuprates are not observable by ESR up to 1100 K, which is $`0.7J/k_\mathrm{B}`$ where $`J/k_\mathrm{B}1600`$ K is only a factor of 2.8 larger.

In Ref.  the authors estimated the $`\chi ^{\mathrm{VV}}`$ values using the reported anisotropic $`g`$ values obtained from ESR measurements, obtaining $`\chi _a^{\mathrm{VV}}=\chi _b^{\mathrm{VV}}=2.4\times 10^5`$ cm3/mol and $`\chi _c^{\mathrm{VV}}=6.6\times 10^5`$ cm3/mol, which were stated to be in agreement with the values from their $`K`$-$`\chi `$ analysis discussed in the following subsection. These values do not agree with our $`T=0`$ values in Eq. (190). In addition, from the $`\chi ^{\mathrm{VV}}`$ values of Ohama et al., one obtains $`\mathrm{\Delta }\chi _{ca}^{\mathrm{VV}}=\mathrm{\Delta }\chi _{cb}^{\mathrm{VV}}=4.2\times 10^5`$ cm3/mol, which are similar in magnitude but opposite in sign to our data in Eq. (190). If the strong change in each of $`\mathrm{\Delta }\chi _{ac}`$ and $`\mathrm{\Delta }\chi _{bc}`$ below $`T_\mathrm{c}`$ in Fig. 40 is due to a respective $`\mathrm{\Delta }\chi _{\alpha \beta }^{\mathrm{VV}}`$ which is strongly temperature dependent in this temperature range, an effect similar to that reported to occur from NMR measurements discussed in the next subsection, it is hard to understand why this change is not reflected in a distinct change in the reported temperature dependent anisotropy of the $`g`$ values at $`T_\mathrm{c}`$ if these $`g`$-value measurements are recording the characteristics of the bulk phase.

##### Anisotropy in the Van Vleck susceptibility from NMR.

From a so-called $`K`$-$`\chi `$ analysis using NMR paramagnetic nuclear resonance shift $`K(T)`$ data, combined with $`\chi (T)`$ measurements, under certain assumptions $`\chi ^{\mathrm{VV}}`$ can be obtained if $`K`$ is proportional to $`\chi `$, with $`T`$ as an implicit parameter. In this way, $`\chi ^{\mathrm{VV}}`$ values have been obtained by Ohama and coworkers for NaV2O5 using 23Na (Ref. ) and 51V (Ref. ) NMR measurements on the same aligned powder sample. The former 23Na study yielded $`\chi _b^{\mathrm{VV}}=23\times 10^5`$ cm3/mol below $`T_\mathrm{c}`$ and $`16\times 10^5`$ cm3/mol above $`T_\mathrm{c}`$, corresponding to a decrease of $`7\times 10^5`$ cm3/mol at $`T_\mathrm{c}`$. Their low temperature value is quite similar to our value in Eq. (190).

The 51V NMR study, carried out above $`T_\mathrm{c}`$, yielded $`\chi _b^{\mathrm{VV}}=2(1)\times 10^5`$ cm3/mol, roughly an order of magnitude smaller than obtained in the authors’ first study (no comment was made about this discrepancy), and in addition gave $`\chi _a^{\mathrm{VV}}=1(1)\times 10^5`$ cm3/mol and $`\chi _c^{\mathrm{VV}}=4(1)\times 10^5`$ cm3/mol. These values are significantly smaller than our values. We note that a $`K`$-$`\chi `$ analysis on the $`d^1\mathrm{V}^{+4}`$ compound VO2 yielded $`\chi ^{\mathrm{VV}}=6.5\times 10^5`$ cm3/mol.

### B Modeling the susceptibility of $`\mathrm{𝐍𝐚𝐕}_\mathrm{𝟐}𝐎_\mathrm{𝟓}`$ above $`𝑻_𝐜`$

Turning now to the experimental $`\chi (T)`$ data in Fig. 38, we have $`T^{\mathrm{max}}370`$ K. Assuming the validity of the Hamiltonian (1), Eq. (67) for the uniform Heisenberg chain yields the exchange constant $`J/k_\mathrm{B}580`$ K. Then the $`g_b`$ value in Eq. (193) and our $`\chi _0`$ values at $`T=0`$ in Table VI below, together with Eqs. (70) and (38), predict that the measured $`\chi ^{\mathrm{max}}40\times 10^5`$ cm3/mol, which is similar to the measured values of $`44`$ and $`48\times 10^5`$ cm3/mol for the two crystals in Fig. 38, respectively. We therefore proceeded to try to fit the data by the uniform chain model. The “Fit”, shown as the solid curve in Fig. 38, is a plot of Eqs. (38), with $`\chi ^{}(t)`$ being the susceptibility of the uniform chain (Fit 2 above) and with the parameters

$`\chi _0=8\times 10^5{\displaystyle \frac{\mathrm{cm}^3}{\mathrm{mol}}},C_{\mathrm{imp}}=0,`$ (194)
(195)
$`g=1.972,{\displaystyle \frac{J}{k_\mathrm{B}}}=580\mathrm{K}.(\mathrm{`}\mathrm{`}\mathrm{Fit}\mathrm{"})`$ (196)
This “Fit” is not really a fit, since we just set the $`g`$ and $`J`$ values to those estimated above and then set $`\chi _0`$ so that the calculated curve is in the vicinity of the data, because no small change in the parameters can bring the theory in agreement with the data. It is clear that adjusting $`\chi _0`$ further will not improve the agreement, nor will including a nonzero impurity Curie constant $`C_{\mathrm{imp}}`$. However, the shapes of the curve and the data are similar, so the agreement can be improved considerably (not shown) by simultaneously decreasing $`\chi _0`$ to $`10\times 10^5\mathrm{cm}^3/\mathrm{mol}`$, which is not possible physically according to Eqs. (38) because it would require the Van Vleck susceptibility to be negative, and increasing $`g`$ to the unphysically large value of $`2.4`$, while keeping $`J`$ constant. These results are in disagreement with the conclusion of Isobe and Ueda who found that the Bonner-Fisher prediction fitted their powder susceptibility data from 50 to 700 K very well assuming $`g=2`$. We can only note that their $`\chi (T)`$ data have not been quantitatively reproduced in either their or others’ subsequent measurements on NaV2O5, including ours, and that the Bonner-Fisher prediction is not accurate at temperatures below $`J/(4k_\mathrm{B})145`$ K as discussed in the Introduction.

Lohmann et al. and Hemberger et al. also previously concluded that the $`\chi (T)`$ of NaV2O5 is not described (below 250 K) by the prediction for the $`S=1/2`$ Heisenberg chain, based on their fits by the Bonner-Fisher prediction to their $`\chi (T)`$ deduced from ESR measurements up to 650 K. They suggested that additional exchange couplings may be required to explain the observed $`\chi (T)`$. We consider this possibility here by modeling the influence of possible interchain spin coupling. Because there are no accurate and generally applicable numerical calculations for this case that we are aware of, we utilize the following simple molecular field theory (MFT) prediction for the spin susceptibility

$$\frac{1}{\chi ^{}(t)}=\frac{1}{\chi _{\mathrm{chain}}^{}(t)}+\frac{z^{}J^{}}{J},$$
(197)
where $`\chi _{\mathrm{chain}}^{}(t)`$ is the reduced spin susceptibility of the isolated quantum $`S=1/2`$ uniform Heisenberg chain (our Fit 2 above). The parameter $`z^{}`$ is the effective number of spins on other chains to which a spin in a given chain is coupled with effective (or average) exchange constant $`J^{}`$. To be consistent with our sign convention for the intrachain exchange constant $`J`$, $`J^{}`$ is positive for AF interactions and negative for ferromagnetic (FM) interactions. Equation (197) is very accurate when $`|z^{}J^{}/J|1`$.

We fitted the $`\chi (T)`$ data above 50 K for the two crystals in Fig. 38 by Eqs. (38) and (197), where we fixed $`g_b=1.972`$ and $`C_{\mathrm{imp}}=0`$ and allowed $`\chi _0`$, $`J`$ and the product $`z^{}J^{}`$ to vary. Very good fits were obtained, for which the fitting parameters are given in Table V. The fits are plotted as the solid curves in Fig. 41. For the parameters of the two crystals taken together, the fitted $`J/k_\mathrm{B}`$ = 584(9) K is the same as deduced above (580 K) from the temperature of the maximum in $`\chi (T)`$, and the fitted $`\chi _0=1.4(16)\times 10^5`$ cm3/mol is similar to our results at low temperatures in Fig. 39(b). A moderately large and negative (FM) interchain coupling $`z^{}J^{}/J=1.26(5)`$ was obtained. This coupling is sufficiently strong that long-range magnetic ordering might be expected, but which is not observed, possibly due to magnetic frustration effects. If the present mean-field interchain coupling analysis is correct, this interchain coupling should be evident in the magnon dispersion relations observable by inelastic magnetic neutron scattering measurements. Indeed, moderately strong dispersions of 1.4 meV in each of two bands perpendicular to the chains have in fact been observed by Yosihama et al. in such measurements on single crystals. It remains to be seen whether the magnitude and sign of the interchain exchange coupling that we infer in the mean-field analysis are consistent with the dispersion relations deduced from the neutron scattering data.

An alternative and/or additional mechanism which can produce a strong deviation of the measured $`\chi (T)`$ of a uniform chain compound from that predicted for Heisenberg uniform and alternating chains is the spin-phonon interaction. At low $`T`$ this interaction can lead to a spin-Peierls transition and can strongly modify $`\chi (T)`$ above $`T_\mathrm{c}`$ from that expected for the Heisenberg chain. Sandvik, Singh, and Campbell carried

out a detailed QMC investigation of a spin-Peierls model in which the spin 1/2 interactions were modified by the presence of dynamical (quantum mechanical) dispersionless Einstein phonons. For particular values of the spin-phonon coupling constant and phonon frequency, they found that the effective exchange constant $`J_{\mathrm{eff}}`$ decreases strongly with increasing $`T`$, and at $`T=0`$ is about 27.3% larger than the bare $`J`$. Perhaps surprisingly, they found however that if the bare $`g`$ factor is reduced by $`7`$% and the bare $`J`$ by $`18`$% in the $`\chi ^{}(t)`$ predicted for the Heisenberg model, this model was then in good agreement with their QMC simulations for temperatures above $`T_\mathrm{c}`$. A recent important extensive study of many finite-temperature properties of the same model using QMC simulations was carried out by Kühne and Löw. They found that for not too low temperatures, the susceptibilities for various Einstein phonon frequencies and spin-phonon coupling constants can all be scaled onto a universal curve, given by that for the uniform Heisenberg chain, using only a suitably defined effective exchange constant $`J_{\mathrm{eff}}>J`$. Contrary to the result of Ref. , they found that a rescaling of the $`g`$ factor was not necessary.

Our experimental results for $`\mathrm{NaV}_2\mathrm{O}_5`$ are not consistent with either of these theoretical studies, because as discussed below Eq. (195) above, to force-fit the Heisenberg chain $`\chi (T)`$ prediction onto the data requires an unphysically large negative value of $`\chi _0`$, as well as an unphysically large increase in $`g`$.

On the other hand, our observed $`\chi (T)`$ does not agree with the Heisenberg chain model (with a temperature-independent $`J`$), and in the next section we simultaneously model the data both above and below $`T_\mathrm{c}`$ within the context of the Heisenberg chain model with a temperature-dependent $`J`$, where we find that $`J(T)`$ above $`T_\mathrm{c}`$ is very similar in form to that deduced in the calculations of Refs.  and . Thus it may be the case that the spin-phonon interaction is indeed important to determining $`\chi (T)`$ in NaV2O5, but where the effects on $`\chi (T)`$ are somewhat different than calculated in the models. In particular, the theoretical predictions may be substantially modified if phonon spectra appropriate to real materials were to be used in the calculations instead of dispersionless Einstein phonons.

### C Simultaneous modeling of the susceptibility of NaV2O5 below and above $`𝑻_𝐜`$

Previous modeling of $`\chi (T)`$ of NaV2O5 to extract the spin gap has usually been done at the lowest temperatures without reference to the magnitude of $`\chi `$ above $`T_\mathrm{c}`$. Here we utilize our fit to the $`\chi ^{}(t)`$ for the Heisenberg chain to extract $`J`$ above $`T_\mathrm{c}`$ from the experimental data. Clearly, since the measured $`\chi (T)`$ above $`T_\mathrm{c}`$ cannot be modeled within this framework using a temperature-independent $`J`$ as shown in the previous section, it follows that if we are to remain within this framework, $`J`$, which is then evidently an effective exchange constant incorporating additional physics of the material, must be temperature dependent. Then with $`J(T)`$ fixed, we derive the $`T`$-dependent spin gap $`\mathrm{\Delta }(T)`$ and exchange alternation parameter $`\delta (T)`$ near and below $`T_\mathrm{c}`$ directly from the measured $`\chi (T)`$ data, which has not, to our knowledge, been carried out before for any system showing a spin-dimerization transition, using our $`\chi ^{}(\alpha ,t)`$ fit function for the alternating chain.

The specific procedure we adopted for modeling the $`\chi _b(T)`$ measurement on each crystal consists of the following six steps, where we fixed $`g_b=1.972`$ in steps 3–5.

(1) The $`\chi (T)`$ from 2 to 10 K is fitted by Eqs. (38), setting $`\chi ^{\mathrm{spin}}=0`$ because of the presence of the spin gap, thereby obtaining the parameters $`\chi _0,C_{\mathrm{imp}}`$, and $`\theta `$.

(2) Using these $`\chi _0,C_{\mathrm{imp}}`$ and $`\theta `$ parameters, we solve for $`J(T)`$ for $`T60`$ K, or for $`T=50`$ K only, using our “Fit 1” function for $`\overline{\chi ^{}}(\overline{t})`$ of the Heisenberg chain, which is one end-point function of our $`\overline{\chi ^{}}(\delta ,\overline{t})`$ fit function, and fit the resulting $`J(T)`$ by a polynomial in $`T`$ for extrapolation below $`T_\mathrm{c}`$; we used the extrapolation function $`J(T)=J(0)+aT^2+bT^3`$.

(3) With this $`J(T)`$, or using $`J(50`$ K) only, we fitted $`\chi (T)`$ from 2 to 20 K, now including $`\chi ^{\mathrm{spin}}(T)`$ for the alternating-exchange chain \[i.e., using our alternating chain $`\overline{\chi ^{}}(\delta ,\overline{t})`$ fit function\] assuming a $`T`$-independent $`\delta `$ (and $`\mathrm{\Delta }`$), and obtain a new set of $`\chi _0,C_{\mathrm{imp}}`$ and $`\theta `$ parameters \[in addition to $`\delta (0)`$\].

(4) Steps 2 and 3 are repeated until convergence is achieved, which takes in practice only one additional iteration. Note that we implicitly assume that $`\chi _0,C_{\mathrm{imp}}`$ and $`\theta `$ are independent of $`T`$, i.e., that the transition(s) at $`T_\mathrm{c}`$ do not affect them.

(5) The experimentally determined molar spin susceptibility $`\chi ^{\mathrm{spin}}(T)`$ is now computed by inserting the final $`\chi _0,C_{\mathrm{imp}}`$ and $`\theta `$ fit parameters into Eq. (183). Then using the fitted $`J(T)`$ or $`J(50`$ K), the $`\delta (T)`$ is computed using our fit function $`\overline{\chi ^{}}(\delta ,\overline{t})`$ for the alternating-exchange chain by finding the root for $`\delta `$, at each data point temperature $`T`$, of

$$\chi _b^{\mathrm{spin}}(T)=\frac{N_\mathrm{A}g_b^2\mu _\mathrm{B}^2}{J_b(T)}\overline{\chi ^{}}[\delta ,\frac{k_\mathrm{B}T}{J_b(T)}].$$
(198)
(6) In a separate step not associated with the fitting procedure in steps 1–5, the spin gap $`\mathrm{\Delta }(T)`$ is computed from $`\delta (T)`$ determined in step 5 using an independently known function $`\overline{\mathrm{\Delta }^{}}(\delta )\mathrm{\Delta }(\delta )/J`$ and our $`J(T)`$ or $`J(50`$ K). We used our $`\overline{\mathrm{\Delta }^{}}(\delta )`$ fit function in Eqs. (27) for this purpose.

We measured $`\chi _b(T)`$ for nine different crystals from four different batches of NaV2O5 and now present illustrative results obtained in each of the above modeling steps 2 to 4 (final iteration), 5 and 6 for three representative crystals. We will follow in graphical form the data modeling through successive steps for these three crystals to show how differences in one property between the crystals may or may not propagate through the next step(s)

of the analysis, but we present the fitting parameters for all of the crystals in Table VI.

The measured $`\chi (T)`$ data below 50 K for the three crystals are shown in Fig. 42(a), where the fits below 20 K in step 4 are shown as the solid curves with parameters in Table VI. Crystals E097A and AL1 are seen to have much lower levels of paramagnetic impurities than E083EF, as reflected in the impurity Curie constant, i.e., the magnitude of the impurity Curie-Weiss upturn at low $`T`$. By subtracting the $`\chi _0`$ and impurity Curie-Weiss terms from the data, the spin susceptibility $`\chi ^{\mathrm{spin}}(T)`$ is obtained for each crystal as shown in Fig. 42(b). These data show good consistency below $`T_\mathrm{c}`$ for the three crystals, despite the differences in the $`\chi _0`$ values, the magnitudes of the Curie-Weiss impurity term and in the $`\chi (T)`$ above $`T_\mathrm{c}`$. The $`J(T)`$ determined for the three crystals in step 2 are shown up to 300 K in Fig. 43. $`J`$ is found to decrease by $`10`$–20 % upon increasing $`T`$ from 60 to 300 K, which when $`T`$ is scaled by $`J`$ is similar to the fractional decrease predicted by Sandvik et al. due to the spin-phonon interaction. It is noteworthy that crystal E083EF, with by far the highest level of paramagnetic defects and/or impurities, also has the largest $`J(T)`$ and the largest change in $`J`$ with $`T`$.

Figures 44(a) and 44(b) show the spin dimerization parameter $`\delta (T)`$ and spin gap $`\mathrm{\Delta }(T)`$ determined for each of the three crystals in the final modeling steps 5 and 6, respectively. Several features of these data are of note. First, there is a rather large variation in the dimerization parameter, $`\delta (0)=0.028`$–0.040, between the three crystals, despite the fact that $`T_\mathrm{c}=33`$–34 K is nearly the same for the different crystals; the most impure crystal E083EF has the smallest $`\delta (0)`$, as might have been expected. Despite this variability, these $`\delta (0)`$ values are all significantly smaller than the three values reported for various samples by different groups as determined using different techniques, which are listed in Table VII along

with other related information.   On the other hand, the

corresponding range of $`\mathrm{\Delta }(0)/k_\mathrm{B}=103(2)`$ K for the three crystals is fractionally much smaller than that of $`\delta (0)`$. We infer that some of the discrepancies between the $`\mathrm{\Delta }(0)`$ values in Table VII reported for NaV2O5 by different groups may arise from differences in, e.g., the types of measurements which are used to determine $`\mathrm{\Delta }(0)`$ and in the different analyses of those data, rather than from different $`\mathrm{\Delta }(0)`$ values in the samples. The variability in $`\delta (0)`$ between the crystals in Fig. 44(a), compared with the lack of much variability in $`\mathrm{\Delta }(0)`$ in Fig. 44(b), evidently arises because $`\delta `$ must be combined with $`J`$ to obtain $`\mathrm{\Delta }`$, and the variations in the first two parameters must largely cancel. Thus, not surprisingly, the low-$`T`$ $`\chi ^{\mathrm{spin}}(T)`$ is governed by the spin gap $`\mathrm{\Delta }`$ and not by $`\delta `$ or $`J`$ separately.

The $`\delta (T)`$ data for our best crystals show very sharp, nearly vertical increases with decreasing $`T`$ at $`T_\mathrm{c}`$. We cannot extract a precise critical exponent $`\beta `$ from our $`\delta (T)`$ data due to the large temperature-dependent background above $`T_\mathrm{c}`$, to be discussed shortly. However, rough fits below $`T_\mathrm{c}`$ by the expression $`\delta (T)(1T/T_\mathrm{c})^\beta `$ gave $`\beta `$ values consistent with the values $`\beta =0.25(10)`$ (Ref. ) from infrared reflectivity measurements, 0.34(8) (Ref. ) from sound velocity measurements along the chain axis and 0.35(8) (Ref. ) from thermal expansion measurements along that axis. We note that these values are a factor of two larger than the value of $`0.15`$ (Ref. ) inferred from x-ray diffuse scattering measurements.

The data in Figs. 44(a) and 44(b) clearly show the existence of spin dimerization fluctuations and a spin pseudogap above $`T_\mathrm{c}`$ for each crystal, respectively, irrespective of the crystal quality as judging from the Curie-Weiss impurity term in the low-$`T`$ $`\chi (T)`$, with magnitudes just above $`T_\mathrm{c}`$ of about 20 % and 40 % of $`\delta (0)`$ and $`\mathrm{\Delta }(0)`$, respectively. This is a robust result, which was obtained

for each of the nine crystals we measured, which does not depend on the precise value of $`J`$ \[and resultant $`\chi ^{\mathrm{spin}}(T,\delta =0)`$\] or the details of how $`J`$ is determined above $`T_\mathrm{c}`$, or even on the detailed formulation of the $`\chi ^{}(\alpha ,t)`$ fit function for the alternating-exchange chain. For example, setting $`J`$ to be a constant, equal to the value at 50 K, yields nearly the same $`\mathrm{\Delta }(T)`$ near $`T_\mathrm{c}`$ as determined using a $`T`$-dependent $`J`$. Similarly, deleting the impurity Curie-Weiss term in the fit to the data above $`T_\mathrm{c}`$ changes the derived $`\chi _0`$ and $`J(T)`$ or $`J`$(50 K) somewhat as well as the detailed temperature dependence of the pseudogap $`\mathrm{\Delta }(T)`$ above $`T_\mathrm{c}`$ but has little influence on the magnitude of $`\mathrm{\Delta }`$ near $`T_\mathrm{c}`$. Further, in a previous version of the QMC and TMRG $`\chi ^{}(\alpha ,t)`$ fit function (not otherwise discussed in this paper), we did not enforce the requirement (iii) in Sec. IV E that the transformed $`\overline{\chi ^{}}(\delta ,\overline{t})`$ satisfy $`\overline{\chi ^{}}(\delta ,\overline{t})/\delta |_{\delta =0}=0`$, and the same fluctuation effects above $`T_\mathrm{c}`$ were found using that fit function as using the present one, although these fluctuations were somewhat reduced in magnitude compared to the present results. Finally, these fluctuations are observable directly in the measured $`\chi (T)`$ data in Fig. 42(a) as a rounding of the susceptibility curves above $`T_\mathrm{c}`$.

From Fig. 44, the fluctuation effects persist up to high temperatures $`T>50`$ K, although the fluctuation amplitudes decrease with increasing $`T`$. Precursor effects above $`T_\mathrm{c}`$ have been reported in x-ray diffuse scattering measurements up to $`90`$ K, in ultrasonic sound velocity and optical absorption measurements up to $`70`$ K, and in specific heat measurements up to $`40`$–50 K, so it is not surprising that spin dimerization parameter fluctuations in Fig. 44(a), and a spin pseudogap in Fig. 44(b) reflecting fluctuations in the spin gap, are found above $`T_\mathrm{c}`$.

### D Specific heat of NaV2O5

In order to correlate the magnetic effects discussed above in NaV2O5 with thermal effects, we have carried out specific heat vs temperature $`C_\mathrm{p}(T)`$ measurements on the same crystal AL1, and a crystal E097 from the same batch as E097A, for which $`\chi (T)`$ data were presented and modeled above. The results from 2 K to 50 K for crystals E097 and AL1 are shown in Fig. 45(a). Over this temperature range, the $`C_\mathrm{p}(T)`$ data for the two crystals agree extremely well, except in the range 33.0–34.2 K, i.e., in the vicinity of the transitions as will be discussed shortly. The shapes of the specific heat anomalies at $`T_\mathrm{c}`$ are not mean-field-like specific heat jumps as observed in, e.g., conventional superconductors, but instead are $`\lambda `$-shaped anomalies. Thus, any attempt to define a (mean-field) “specific heat jump at $`T_\mathrm{c}`$” is fraught with ambiguity. These shapes are retained in plots of $`C_\mathrm{p}(T)/T`$ vs $`T`$ as shown in Fig. 46(a). This $`\lambda `$ shape has been observed previously, and variously attributed to fluctuation effects or a possible smeared-out first order transition. In view

of the coupled structural, charge-ordering and spin dimerization transitions at $`T_\mathrm{c}`$ in NaV2O5 as discussed in the Introduction, their relative contributions to the specific heat anomalies are not clear, if indeed their contributions can be uniquely distinguished.

Expanded plots of $`C_\mathrm{p}(T)`$ and $`C_\mathrm{p}(T)/T`$ versus $`T`$, shown in Figs. 45(b) and 46(b), respectively, reveal a sharp high peak at 33.4 K for crystal AL1, which is slightly split by $`0.1`$ K in spite of the fact that the overall height of the anomaly is much larger than pre-

viously reported for any crystal of NaV2O5. Two peaks are also observed for crystal E097, at 33.4 K and 33.8 K, which are more widely separated than for crystal AL1. From Fig. 46(b), the entropy under the anomaly(ies) for each crystal is about the same (see below). Comparing these results with the $`\delta (T)`$ and $`\mathrm{\Delta }(T)`$ data in Fig. 44, the larger splitting of the $`C_\mathrm{p}(T)`$ peak for crystal E097 does not result in any major difference in the magnetic order parameter properties between the two crystals, although the transition onset is slightly rounded for crystal E097A compared to AL1. By using the Fisher relation, $`[\chi (T)T]/TC(T)`$ where $`C(T)`$ is the magnetic contribution to the specific heat, one obtains results which show the same features near $`T_\mathrm{c}`$ as does the specific heat, as shown in Fig. 45(c). Thus, careful scrutiny of the magnetic properties can reveal the fine detail observed in the specific heat near $`T_\mathrm{c}`$. In particular, this comparison suggests that both anomalies in the specific heat near $`T_\mathrm{c}`$ for each crystal are associated with and/or reflected by magnetic effects.

The splitting of the transition into two apparent transitions that we report here was previously observed in thermal expansion, but not seen in their specific heat, measurements of a crystal by Köppen et al. The detailed origin of the transition splitting, and more fundamentally whether the splitting is instrinsic to ideal crystallographically ordered NaV2O5, remain to be clarified. An essential feature that any explanation must account for is that the temperature splitting between the two transitions in a crystal varies from crystal to crystal.

##### Modeling.

In this section we will only consider the model utilized above for analyzing our $`\chi (T)`$ data, in which $`\mathrm{NaV}_2\mathrm{O}_5`$ consists, effectively, of isolated $`S=1/2`$ uniform or (below $`T_\mathrm{c}`$) alternating-exchange Heisenberg chains, where the (average) exchange constant $`J`$ shows, at most, only a smooth and relatively small change below $`T_\mathrm{c}`$. For reasons which will become clear below, unfortunately we cannot use our specific heat data to extract detailed information about the magnetic subsystem in $`\mathrm{NaV}_2\mathrm{O}_5`$. However, other types of important information about the thermodynamics will be derived using various of the theoretical results presented and discussed previously in this paper.

There have been two reports deriving the spin gap from $`C_\mathrm{p}(T)`$ data at $`T15`$ K. We first discuss the limits of this type of analysis. Using $`J(0)/k_\mathrm{B}=600`$ K, $`\delta (0)=0.040`$, and $`\mathrm{\Delta }/k_\mathrm{B}=100`$ K (see Table VI), Eqs. (64) and (65) predict that the magnetic specific heat $`C(T)`$ in the dimerized phase at low temperatures $`T(\mathrm{\Delta }/k_\mathrm{B},T_\mathrm{c})`$ is

$`C(T)=`$ $`1.0{\displaystyle \frac{\mathrm{J}}{\mathrm{mol}\mathrm{K}}}\left({\displaystyle \frac{100}{T}}\right)^{3/2}`$ (200)
$`\times \left[1+{\displaystyle \frac{T}{100}}+{\displaystyle \frac{3}{4}}\left({\displaystyle \frac{T}{100}}\right)^2\right]\mathrm{e}^{100/T},`$
with $`T`$ in units of K. Equation (200) predicts that $`C(15\mathrm{K})=0.026`$ J/mol K, which is about 40 times smaller than the observed $`C_\mathrm{p}(15\mathrm{K})1`$ J/mol K (which must therefore be due to the lattice contribution) and hence is unresolvable at such low temperatures. Within this model, we must therefore conclude that the previous estimates of the spin gap based upon modeling the low temperature specific heat were most likely artifacts of modeling the lattice specific heat. This can happen if one does not utilize the fact that the prefactor to the activated exponential term of the magnetic contribution $`C(T)`$ is not an independently adjustable parameter, but is rather determined by the spin gap itself as we have previously demonstrated and emphasized in Sec. II C 3.

A related question is whether the entropy associated with the transition(s) at $`T_\mathrm{c}`$ can be associated solely with the magnetic subsystem. The minimum possible estimate of the entropy of the transition is obtained from the $`C_\mathrm{p}(T)/T`$ vs $`T`$ data in Fig. 46(b) by drawing a horizontal line from the data at the $`C_\mathrm{p}(T)/T`$ minimum at $`35.0`$ K, just above $`T_\mathrm{c}`$, to the data that the line intersects with below $`T_\mathrm{c}`$ at $`30.6`$ K, and then computing the area between the line and the peak(s) above the line. In this way we obtain a value of 0.397 J/mol K for crystal E097 and 0.375 J/mol K for crystal AL1. On the other hand, the maximum magnetic entropy of the $`S=1/2`$ uniform chain subsystem at $`T_\mathrm{c}`$, using rough values $`J/k_\mathrm{B}=600`$ K

and $`T_\mathrm{c}=34`$ K, is roughly $`S(T_\mathrm{c})(2R/3)(k_\mathrm{B}T_\mathrm{c}/J)=0.31`$ J/mol K. Thus, the specific heat $`\lambda `$ anomaly at $`T_\mathrm{c}`$ cannot arise solely from the magnetic subsystem, since the minimum possible entropy of the transition is significantly greater than the maximum possible magnetic entropy at $`T_\mathrm{c}`$. At the least, the remaining entropy must therefore be due to the crystallographic and/or charge-ordering transitions which occur at or close to the spin dimerization transition temperature as discussed in the Introduction.

A potentially definitive and effective way to proceed from this point would be to quantitatively determine the magnetic contribution $`C(T)`$ to the measured specific heat $`C_\mathrm{p}(T)`$ from $`\chi (T)`$ at and near $`T_\mathrm{c}`$, using a relationship between $`\chi (T)`$ and $`C(T)`$ such as the Fisher relation cited above, and then compare this result with $`C_\mathrm{p}(T)`$. From a comparison of Figs. 45(b) and 45(c), it seems clear that such a relation must exist, at least for temperatures near $`T_\mathrm{c}`$, but the relationship between $`\chi (T)`$ and $`C(T)`$ near spin dimerization transitions has not yet been worked out theoretically.

In the absence of such a formulation, we proceed to estimate the change in the specific heat associated with the transition(s). In order to do this modeling, we must fit $`C_\mathrm{p}(T)`$ to higher temperatures than we have been discussing so far. The $`C_\mathrm{p}(T)`$ data from 2 to 100 K for $`\mathrm{NaV}_2\mathrm{O}_5`$ crystal E097 are shown as the open circles in Fig. 47. As noted above, except in the immediate vicinity of $`T_\mathrm{c}`$ the $`C_\mathrm{p}(T)`$ data for crystal AL1 are nearly identical with those for crystal E097 up to at least 50 K, so it will suffice to model the data for crystal E097. The four modeling steps and the assumptions we employed are as follows.

(1) We assume that critical and other order parameter fluctuations associated with the transition(s) at $`T_\mathrm{c}`$ make a negligible contribution to $`C_\mathrm{p}(T)`$ over some specified high temperature ($`TT_\mathrm{c}`$) range. By subtracting the known magnetic contribution $`C(T)`$ due to isolated Heisenberg chains \[obtained using our fit function for $`C(k_\mathrm{B}T/J)`$\] in this temperature range, we obtain the background lattice contribution $`C^{\mathrm{lat}}(T)`$ in the high temperature region. Also, since we have shown that $`C(T)`$ is negligible for $`T15`$ K, the measured $`C_\mathrm{p}(T)`$ in this $`T`$ range is assumed to be identical to $`C^{\mathrm{lat}}(T)`$ at these temperatures (we again neglect the possible but unknown specific heats associated with possible order parameter fluctuations in this range). Thus we obtain background lattice specific heats $`C^{\mathrm{lat}}(T)`$ in high and low temperature ranges which are assumed unaffected by the transition(s) and associated order parameter fluctuations.

(2) We interpolate between the $`C^{\mathrm{lat}}(T)`$ determined in step 1 in the low- and high-temperature ranges to obtain, in the intermediate temperature range, what $`C^{\mathrm{lat}}(T)`$ would have been in the absence of the transition(s) and associated order parameter fluctuations.

(3) We add the $`C(T)`$ for isolated chains, used in step 1, back to the $`C^{\mathrm{lat}}(T)`$ derived in step 2 over the full temperature range of the measurements. This is the total background specific heat that would have occurred in the absence of the transition(s) and associated order parameter fluctuations. Then we subtract the total calculated background specific heat from the measured $`C_\mathrm{p}(T)`$ data. This difference $`\mathrm{\Delta }C(T)`$ should hopefully be a reasonable estimate of the change in the specific heat associated with the transition and order parameter fluctuations, including all lattice, charge and spin contributions. $`\mathrm{\Delta }C(T)`$ must go to zero, by construction, at temperatures above the lower end of the high temperature region fitted in step 2.

(4) Finally we integrate $`\mathrm{\Delta }C/T`$ with $`T`$ up to and beyond $`T_\mathrm{c}`$ to obtain the change in entropy $`\mathrm{\Delta }S(T)`$ associated with the transition and order parameter fluctuations. $`\mathrm{\Delta }S(T)`$ must become constant, by construction, at temperatures above the lower end of the high temperature region fitted in step 2.

In the following we will present and discuss the results in each of the four steps of our modeling program described above.

##### Step 1.

Here we first use our $`C(T)`$ fit function for the numerical $`C(T)`$ data, which was given in Eqs. (18), to extract $`C^{\mathrm{lat}}(T)`$ in the high-temperature region above $`T_\mathrm{c}`$. For consistency with our analysis of the susceptibility in Sec. VIII C, we use the temperature-dependent $`J(T)`$ derived in that section for crystal E097A (see Fig. 43) when computing $`C(T)`$. The background $`C(T)`$ thus estimated for crystal E097, i.e., the values which would have been observed if no transition(s) at $`T_\mathrm{c}`$ or associated order parameter fluctuations had occurred, is shown in Fig. 48. Comparison of these data with the measured $`C_\mathrm{p}(T)`$ data in Figs. 45 and 47 shows that this $`C(T)`$ is a small, but non-negligible ($`1`$%), fraction of $`C_\mathrm{p}(T)`$ above $`T_\mathrm{c}`$. On the other hand, $`C(T)`$ is much larger than the observed $`C_\mathrm{p}(T)`$ at low temperatures, because in this temperature range $`CT`$ whereas $`C_\mathrm{p}(T)C^{\mathrm{lat}}(T)T^3`$. The $`C^{\mathrm{lat}}(T)=C_\mathrm{p}(T)C(T)`$ in the high temperature (60–100 K) region is shown in Fig. 49, together with

$`C^{\mathrm{lat}}(T)C_\mathrm{p}(T)`$ in the low temperature (2–15 K) region.

##### Step 2.

In this step we must interpolate $`C^{\mathrm{lat}}(T)`$ between the low- and high-temperature regions, i.e., in a broad temperature range spanning the transition region. The best way to do this would be to determine $`C^{\mathrm{lat}}(T)`$ directly from $`C_\mathrm{p}(T)`$ measurements on a suitably chosen reference compound, but such measurements have not yet been done. At first sight, a physically realistic possibility might be to interpolate the low and high temperature $`C^{\mathrm{lat}}(T)`$ data using the Debye specific heat function; however, this method is questionable because the Debye temperature $`\mathrm{\Theta }_\mathrm{D}`$ in real materials can be rather strongly temperature dependent within the temperature range of interest here. The Debye function for the molar lattice specific heat at constant volume $`C^{\mathrm{Debye}}(T)`$ is given by

$$C^{\mathrm{Debye}}(T)=9rR\left(\frac{T}{\mathrm{\Theta }_\mathrm{D}}\right)^3_0^{\mathrm{\Theta }_\mathrm{D}/T}\frac{x^4\mathrm{e}^x}{(\mathrm{e}^x1)^2}𝑑x,$$
(201)
where $`r`$ is the number of atoms per formula unit ($`r=8`$ here) and $`R`$ is the molar gas constant. We attempted to fit our $`C^{\mathrm{lat}}(T)`$ data for the temperature ranges 2–15 K and 40–100 K to 80–100 K by Eq. (201). The fits parametrized the data very poorly. We obtained a more reasonable fit by allowing $`r`$ to be a fitting parameter, yielding a fitted value $`r4`$, but the data were still poorly fitted, due to too much curvature in the Debye function in the high temperature region. Therefore, we were led to interpolating between the low- and high-temperature regions using a simple polynomial interpolation function.

To obtain the background lattice specific heat interpolation function, we fitted the combined $`C^{\mathrm{lat}}(T)`$ data (a total of 141 data points) in the low and high temperature ranges 2–15 K and 60–101 K, respectively, by polynomials of the form

$$C^{\mathrm{lat}}(T)=\underset{n=3}{\overset{n^{\mathrm{max}}}{}}c_nT^n.$$
(202)
The minimum summation index $`n=3`$ is set by the expected Debye low-temperature $`T^3`$ behavior of the lattice specific heat. The maximum value $`n^{\mathrm{max}}`$ was varied to see how the fit parameters and variance changed. In addition, for checking the final fits we fitted the $`C^{\mathrm{lat}}(T)`$ data in the 2–15 K low-$`T`$ range together with $`C^{\mathrm{lat}}(T)`$ data in a high-temperature range varying from 40–101 K to 90–101 K. We found that the most stable fits were for $`n^{\mathrm{max}}=7`$ and 8. For both values, the fit did not visibly change when the lower limit of the upper temperature range of the fitted data was varied from 60 to 70 K. We chose to use the fit for $`n^{\mathrm{max}}=7`$ because in this case the fit was also stable for lower limits of 50 and 80 K. This stability allows one to be confident that the interpolation of the fit between the fitted low- and high-temperature ranges is an accurate representation of the background lattice specific heat in the interpolated intermediate temperature range. The fit for the temperature ranges 2–15 K and 60–101 K is shown as the solid curve in Fig. 49. The absolute rms deviation of this fit from the fitted data is quite small, $`\sigma _{\mathrm{rms}}=0.046`$ J/mol K. The curve over the full temperature range 2–101 K represents the background lattice specific heat $`C^{\mathrm{lat}}(T)`$ expected in the absence of any transitions or order parameter fluctuations.

##### Step 3.

Adding the magnetic background specific heat contribution $`C(T)`$ obtained in step 1 to the lattice background specific heat contribution $`C^{\mathrm{lat}}(T)`$ obtained in step 2 gives the total background specific heat, which is plotted as the solid curve in Fig. 47. We reiterate that this background is interpreted as the specific heat that would have been observed had the transition(s) at $`T_\mathrm{c}`$

and associated order parameter fluctuations not occurred. The difference $`\mathrm{\Delta }C`$ between the measured $`C_\mathrm{p}(T)`$ and the total background specific heat is plotted versus temperature in Fig. 50(a). As would have been qualitatively anticipated, $`\mathrm{\Delta }C`$ is negative below about 16 K due to the loss of magnetic specific heat at low temperatures arising from the opening of the spin gap at $`T_\mathrm{c}`$. This negative $`\mathrm{\Delta }C`$ does not arise from a problem in our polynomial interpolation $`C^{\mathrm{lat}}(T)`$ fit function or from our $`C(T)`$ function; these functions are both positive for all $`T>0`$. Since the magnetic background contribution is proportional to $`T`$ and the lattice background contribution \[which is assumed not to change below 15 K due to the occurrence of the transition(s)\] is proportional to $`T^3`$ at low $`T`$, opening a spin gap at $`T_\mathrm{c}`$ must necessarily lead to a negative $`\mathrm{\Delta }C`$ at sufficiently low temperatures since the magnetic contribution then becomes exponentially small there.

##### Step 4.

Finally, we can compute the change $`\mathrm{\Delta }S`$ in the total entropy of the system versus temperature due to the transition(s) and associated order parameter fluctuations by integrating $`\mathrm{\Delta }C(T)`$ from step 3 according to $`\mathrm{\Delta }S(T)=_0^T[\mathrm{\Delta }C(T)/T]𝑑T`$. The result is shown in Fig. 50(b). The entropy change is negative below about 22 K, due to the loss of magnetic entropy at low temperatures associated with the loss of magnetic specific heat as just discussed. From conservation of magnetic entropy, this lost entropy must reappear at higher temperatures.

By construction, step 2 requires that $`\mathrm{\Delta }C(T>60`$ K) = 0 and consequently $`\mathrm{\Delta }S(T>60`$ K) = const. This requirement is not desirable, but we had to enforce it to ensure that the $`\mathrm{\Delta }C(T)`$ and $`\mathrm{\Delta }S(T)`$ derived at lower temperatures were accurate. Since the effects of the order parameter fluctuations are likely to continue to be present at temperatures higher than 60 K, the $`\mathrm{\Delta }C(T)`$ and $`\mathrm{\Delta }S(T)`$ at temperatures at and near 60 K in Fig. 50 are lower limits.

The net change in the entropy at 60 K in Fig. 50(b) due to the occurrence of the transition(s) at $`T_\mathrm{c}34`$ K and associated order parameter fluctuations above and below $`T_\mathrm{c}`$ is $`\mathrm{\Delta }S(60`$ K) = 2.28 J/mol K. This is far larger than the maximum possible change $`\mathrm{\Delta }S_{\mathrm{mag}}^{\mathrm{max}}=0.556`$ J/mol K in the magnetic entropy at this temperature obtained from Fig. 48, where this value is just the maximum possible entropy of the magnetic subsystem at this temperature, confirming our qualitative conclusion above based on very rough arguments. In particular, our quantitative analysis indicates that at least 76% of the entropy change at 60 K must arise from the lattice and charge degrees of freedom, and only a minor fraction ($`<24`$ %) from the magnetic degrees of freedom. Similarly, at $`T_\mathrm{c}=33.7`$ K, we obtain $`\mathrm{\Delta }S=1.38`$ J/mol K and $`\mathrm{\Delta }S_{\mathrm{mag}}^{\mathrm{max}}=0.311`$ J/mol K, yielding $`\mathrm{\Delta }S_{\mathrm{mag}}^{\mathrm{max}}/\mathrm{\Delta }S23`$ % at $`T_\mathrm{c}`$.

As a closing remark for this section, it is clear from Fig. 50 and the discussion in the above two paragraphs that $`\mathrm{\Delta }C`$ and $`\mathrm{\Delta }S`$ do not saturate to their respective high temperature limiting values until a temperature of at least 60 K is reached, which is almost twice $`T_\mathrm{c}`$. The present analysis of the thermal behavior of $`\mathrm{NaV}_2\mathrm{O}_5`$ thus lends strong support to our independent analysis and interpretation of our magnetic susceptibility data for this compound in Sec. VIII C.

## IX Summary and Concluding Discussion

We have shown that the high-accuracy numerical Bethe ansatz calculations of the magnetic susceptibility $`\chi ^{}(t)`$ for the $`S=1/2`$ uniform Heisenberg chain by Klümper and Johnston are in excellent agreement with the theory of Lukyanov over 18 decades of temperature at low temperatures. An independent high precision empirical fit to these data was obtained over 25 decades of temperature which we found useful to determine the accuracy of our TMRG $`\chi ^{}(t)`$ calculations. The magnetic specific heat data for the uniform chain at very low temperatures was also compared with the theoretical predictions of Lukyanov, and extremely good agreement was found over many decades in temperature. We formulated an empirical fit function for these data which is highly accurate over a temperature range spanning 25 orders of magnitude; the infinite temperature entropy calculated using this fit function is within 8 parts in $`10^8`$ of the exact value. We used both of the above fit functions to model our respective experimental data for $`\mathrm{NaV}_2\mathrm{O}_5`$ in later sections of the paper. We expect that they will be useful to other theorists and experimentalists as well.

We have carried out extensive QMC simulations and TMRG calculations of $`\chi ^{}(\alpha ,t)`$ for the spin $`S=1/2`$ antiferromagnetic alternating-exchange Heisenberg chain for reduced temperatures $`tk_\mathrm{B}T/J_1`$ from 0.002 to 10 and alternation parameters $`\alpha J_2/J_1`$ from 0.05 to 1, where $`J_1(J_2)`$ is the larger (smaller) of the two alternating exchange constants. An accurate global two-dimensional ($`\alpha ,t`$) fit to these combined data was obtained, constrained by the fitting parameters for the accurately known $`\chi ^{}(t)`$ for the $`\alpha `$ parameter end points, the dimer ($`\alpha =0`$) and the uniform chain ($`\alpha =1`$), resulting in an accurate fit function over the entire range $`0\alpha 1`$ of the alternation parameter. Our fit function incorporates the first four terms of the exact high-temperature series expansion in powers of $`1/t`$, which allows accurate extrapolation to arbitrarily high temperatures. This function should prove useful for many applications including the modeling of experimental $`\chi (T)`$ data as we have shown.

Our $`\chi ^{}(\alpha ,t)`$ fit function for the alternating chain can be easily transformed (as we have done) into an equivalent fit function $`\overline{\chi ^{}}(\delta ,\overline{t})`$ in the two variables $`\delta (J_1J_2)/(2J)`$ and $`\overline{t}k_\mathrm{B}T/J`$, where the average exchange constant is $`J=(J_1+J_2)/2`$. This is a more appropriate function for analyzing experimental $`\chi (T)`$ data for $`S=1/2`$ Heisenberg chain compounds showing dimerization transitions (such as a spin-Peierls transition) which result in an alternating-exchange chain with a small value of $`\delta `$ at low temperatures. Once $`J`$ has been determined by fitting our function for $`\delta =0`$ to the experimentally determined spin susceptibility $`\chi ^{\mathrm{spin}}(T)`$ data above the transition temperature, the alternation parameter $`\delta `$ is uniquely determined by our fit function at each temperature below the transition temperature from the value of $`\chi ^{\mathrm{spin}}`$ at that temperature. One can then find the spin gap $`\mathrm{\Delta }(T)`$ using an independently known $`\overline{\mathrm{\Delta }^{}}(\delta )`$.

Our QMC and TMRG data and fit for $`\chi ^{}(\alpha ,t)`$ are in good agreement with previous calculations based on exact diagonalization of the nearest neighbor Heisenberg Hamiltonian for short chains with $`\alpha =0.2`$, 0.4, 0.6, 0.7, and 0.8, extrapolated to the thermodynamic limit, by Barnes and Riera. However, the numerical and analytical theoretical predictions of Bulaevskii, which have been used extensively in the past by experimentalists to model their $`\chi (T)`$ data for weakly-dimerized chain compounds, are found to be in poor agreement with our results and should be abandoned for such use in favor of our fit function. Similarly, the previously used fit function for the Bonner-Fisher calculation of $`\chi ^{}(t)`$ for the uniform chain ($`\alpha =1`$) should be replaced by one of our two fit functions for the most accurate calculation to date of $`\chi ^{}(t)`$ for the uniform chain.

An important theoretical issue in the study of the alternating exchange chain is how the spin gap $`\overline{\mathrm{\Delta }^{}}(\delta )`$ evolves as the uniform chain limit is approached ($`\delta 0,\alpha 1`$). We formulated a fit function for the temperature dependence of our TMRG susceptibility $`\chi ^{}(\alpha ,t)`$ calculations at low temperatures, which was used to extract the dependence $`\overline{\mathrm{\Delta }^{}}(\delta )`$ in this regime. We find that the asymptotic critical regime is not entered until, at least, $`\delta 0.005`$ ($`\alpha 0.99`$). We compared our spin gap data with many literature data. We formulated a fit function for our spin gap data together with those of Barnes, Riera, and Tennant which quite accurately covers the entire range $`0\delta 1`$.

In the remainder of this paper, we showed how the above theoretical results could be used to obtain detailed information about real systems. As a specific illustration, we carried out a detailed analysis of our experimental $`\chi (T)`$ and specific heat $`C_\mathrm{p}(T)`$ data for $`\mathrm{NaV}_2\mathrm{O}_5`$ crystals. This compound shows a transition to a spin dimerized state below the transition temperature $`T_\mathrm{c}34`$ K. We used one of our two $`\chi ^{}(t)`$ fit functions for the uniform Heisenberg chain to model the $`\chi (T)`$ above $`T_\mathrm{c}`$, where we found that the experimental $`\chi (T)`$ is not in quantitative agreement with the prediction for the uniform Heisenberg chain. A model incorporating a mean-field ferromagnetic interchain coupling between quantum $`S=1/2`$ Heisenberg chains fits the experimental data very well with reasonable parameters. It remains to be seen whether the inelastic neutron scattering measurements of the magnon dispersion relations are consistent with our derived intrachain and interchain exchange constants.

In an alternate description, we modeled the deviation in the measured $`\chi (T)`$ of NaV2O5 above 60 K $`>T_\mathrm{c}`$ from the Heisenberg chain model (with fixed exchange constant $`J`$) as due to a temperature-dependent $`J`$. We found that this $`J`$ decreases with increasing $`T`$ up to 300 K in a manner very similar to $`J_{\mathrm{eff}}(T)`$ predicted by Sandvik, Singh and Campbell and Kühne and Löw for the spin-Peierls chain. Our $`J(T)`$ cannot however be compared directly with their $`J_{\mathrm{eff}}(T)`$ because the two quantities are defined differently. They found that by defining an appropriate effective exchange constant $`J_{\mathrm{eff}}`$, their resulting susceptibility $`\chi (k_\mathrm{B}T/J_{\mathrm{eff}})`$ is universal at the higher temperatures for various Einstein phonon frequencies and spin-phonon coupling constants. This function agrees well with the $`\chi (k_\mathrm{B}T/J)`$ for the $`S=1/2`$ AF uniform Heisenberg chain at these temperatures. As we discussed, these $`\chi (T)`$ calculations are not applicable to NaV2O5, possibly because the calculations do not incorporate realistic phonon spectra.

Below $`T_\mathrm{c}`$, we used the $`J(T)`$ extrapolated from above 60 K and our global $`\chi ^{}(\alpha ,t)`$ fit function for the alternating Heisenberg chain to determine the temperature-dependent alternation parameter $`\delta (T)`$, and then the spin gap $`\mathrm{\Delta }(T)`$ from $`\delta (T)`$, directly from the $`\chi (T)`$ data. We find that the $`\mathrm{\Delta }(0)/k_\mathrm{B}`$ values for nine single crystals of NaV2O5 are in the range 103(2) K. This result is in agreement, within the errors, with many previous analyses of data from various types of measurements for this compound by other groups. However, our values of $`\delta (0)=0.034(6)`$ for various crystals are significantly smaller than previous estimates. We note that the two estimates with $`\delta (0)0.1`$ in Table VII were obtained using Bulaevskii’s theory for the alternating-exchange chain, which we have shown is not accurate at low temperatures in the relevant alternation parameter range.

The dispersion of two one-magnon branches perpendicular to the chains observed in the neutron scattering measurements has been recently explained quantitatively by Gros and Valenti assuming that a zig-zag charge ordering transition occurs at $`T_\mathrm{c}`$. They also predict that $`\delta (0)0.034`$. This is within our range of $`\delta (0)`$ values in spite of the fact that we assumed that $`J(T)`$ is either constant or increases slightly with decreasing $`T`$ below $`T_\mathrm{c}`$, contrary to their prediction that $`J`$ decreases below $`T_\mathrm{c}`$. Gros and Valenti made no predictions for $`\chi (T)`$, $`\delta (T)`$, $`\mathrm{\Delta }(T)`$ or $`C(T)`$, so comparisons with our results for these quantities are not possible. We note that Klümper, Raupach, and Schönfeld obtained a good fit to the $`\chi (T)`$ data below $`T_\mathrm{c}`$ for the spin-Peierls compound CuGeO3 within the context of a spin-Peierls model containing frustrating second-neighbor interactions and static spin-phonon coupling.

We discovered that $`\mathrm{\Delta }(T)`$ \[and $`\delta (T)`$\] of NaV2O5 does not go to zero at $`T_\mathrm{c}`$, indicating the existence of a spin pseudogap above $`T_\mathrm{c}`$ with a large magnitude just above $`T_\mathrm{c}`$ of $`40`$% of $`\mathrm{\Delta }(0)`$; the pseudogap is present up to at least 50 K with a magnitude decreasing with increasing $`T`$ above $`T_\mathrm{c}`$. To our knowledge, this pseudogap has not been reported previously, and there are as yet no theoretical predictions for the magnitude or temperature dependence of this pseudogap. The pseudogap is strongly reminiscent of the spin pseudogap derived by one of us using $`\chi (T)`$ measurements above the transition temperature of inorganic quasi-one-dimensional charge density wave compounds, as predicted theoretically by Lee, Rice, and Anderson long before those observations were made. Similar to that case, in the present system one may think of the pseudogap as the rms fluctuation in the spin gap above $`T_\mathrm{c}`$, with an associated reduction in the magnon density of states at low energy. In this interpretation, the pseudogap in NaV2O5 should be observable in high resolution quasielastic neutron scattering and other spectroscopic measurements probing the low energy magnetic excitations.

Finally, we carried out an extensive modeling study of our specific heat data for NaV2O5 crystals, using the same model that we used to analyze our susceptibility data. The most important part of this study is that we have been able to determine a limit on the relative contributions of the magnetic and lattice/charge degrees of freedom to the entropy associated with the transition(s) at $`T_\mathrm{c}`$. We find that at least 77 % of the change in the entropy at $`T_\mathrm{c}`$ must arise from the lattice and/or charge degrees of freedom, to which the spin degrees of freedom must of course be coupled, and that the spin degrees of freedom themselves contribute less than 23 % of this entropy change. Our results also indicate that order parameter fluctuation effects are important in the specific heat up to at least 60 K, strongly confirming the above similar and independent conclusion based on our modeling of our magnetic susceptibility data for the same crystals.

## Acknowledgments

We thank E. Brücher and C. Lin for help with the sample preparation, and C. Song for assistance with Laue x-ray diffraction measurements. We are grateful to S. Eggert and T. Barnes for providing the numerical $`\chi ^{}(t)`$ calculation results for the uniform chain in Ref.  and the alternating chain in Ref. , respectively, and to D. Poilblanc and G. S. Uhrig for sending us their $`\overline{\mathrm{\Delta }^{}}(\delta )`$ data in Refs.  and , respectively. We are grateful to M. Greven and X. Zotos for helpful discussions, and to A. A. Zvyagin for helpful correspondence. One of us (D.C.J.) thanks the Max-Planck-Institut für Festkörperforschung, Stuttgart, where this work was started, for kind hospitality. Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. The work at Ames was supported by the Director for Energy Research, Office of Basic Energy Sciences. The QMC program was written in C++ using a parallelizing Monte Carlo library developed by one of the authors. The QMC simulations by M.T. were performed on the Hitachi SR2201 massively parallel computer of the University of Tokyo and on the IBM SP-2 of the Competence Center for Computational Chemistry of ETH Zürich. X.W. acknowledges Swiss National Funding Grant No. 20-49486.96. A.K. acknowledges financial support by the Deutsche Forschungsgemeinschaft under Grant No. Kl 645/3 and support by the research program of the Sonderforschungsbereich 341, Köln-Aachen-Jülich.