Multiplicative logarithmic corrections to quantum criticality in three-dimensional dimerized antiferromagnets

# Multiplicative logarithmic corrections to quantum criticality in three-dimensional dimerized antiferromagnets

Yan Qi Qin Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
B. Normand Department of Physics, Renmin University of China, Beijing 100872, China    Anders W. Sandvik Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA
Zi Yang Meng Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
July 4, 2019
###### Abstract

We investigate the quantum phase transition in an dimerized Heisenberg antiferromagnet in three spatial dimensions. By performing large-scale quantum Monte Carlo simulations and detailed finite-size scaling analyses, we obtain high-precision results for the quantum critical properties at the transition from the magnetically disordered dimer-singlet phase to the antiferromagnetically ordered Néel phase. This transition breaks O() symmetry with in dimensions. This is the upper critical dimension, where multiplicative logarithmic corrections to the leading mean-field critical properties are expected; we extract these corrections, establishing their precise forms for both the zero-temperature staggered magnetization, , and the Néel temperature, . We present a scaling Ansatz for , including logarithmic corrections, which agrees with our data and indicates exact linearity with , implying a complete decoupling of quantum and thermal fluctuation effects even arbitrarily close to the quantum critical point. We also demonstrate the predicted -independent leading and subleading logarithmic corrections in the size-dependence of the staggered magnetic susceptibility. These logarithmic scaling forms have not previously been identified or verified by unbiased numerical methods and we discuss their relevance to experimental studies of dimerized quantum antiferromagnets such as TlCuCl.

###### pacs:
75.10.Jm, 75.40.Cx, 75.40.Mg

## I Introduction

Antiferromagnetic insulators exhibit a multitude of fundamental phenomena in the neighborhood of the phase transitions separating their magnetically ordered ground states from different types of quantum paramagnetic phase. These quantum phase transitions (QPTs) occur at temperature as a consequence of non-thermal parameters (examples include magnetic fields, applied pressure, and dopant concentration) that act to change the effect of quantum mechanical fluctuations Sachdev2008 (); Kaul2013 (). At finite temperatures, a further dimension is opened in the presence of both quantum and classical (thermal) fluctuations, and the rich physics arising from their interplay includes all the properties of the quantum critical (QC) regime rs ().

Experimentally, the material in which the most detailed study of intertwined classical and quantum critical behavior has been performed is TlCuCl. This compound is composed of antiferromagnetically coupled pairs of Cu ions (), which tend to form dimer singlets and have antiferromagnetic interdimer couplings in all three spatial dimensions () Matsumoto2004 (). At ambient pressure and zero field, TlCuCl is a nonmagnetic insulator with a gap of 0.63 meV to triplet spin excitations. As a consequence of this small gap, an applied magnetic field of 5.4 T is sufficient to drive the system to an ordered antiferromagnetic state, through a QPT in the Bose-Einstein universality class rgrt (). A relatively small applied hydrostatic pressure, kbar, is also sufficient to create an antiferromagnetically ordered state rrs2004 (), through a QPT in the three-dimensional (3D) O(3) universality class due to spontaneous breaking of the SU(2) spin symmetry (which is further reduced in TlCuCl by a weak unixial anisotropy, making the universality class 3D XY).

The elementary excitations on the ordered side of the zero-field quantum critical point (QCP) are gapless spin waves, the Goldstone modes associated with spontaneous breaking of spin-rotational symmetry. On the disordered side, quantum fluctuations, towards spin-singlet formation on the dimers, suppress the long-range antiferromagnetic order, restoring the symmetry and ensuring that all excitations are gapped. This evolution of the excitation spectrum in TlCuCl has been measured in Ref. Ruegg2008 (). At finite temperatures on the ordered side, a classical phase transition occurs at the Néel temperature, , where the long-ranged magnetic order is “melted” not by quantum but by thermal fluctuations. At finite temperatures around the QCP, the combination of quantum and thermal fluctuations creates the QC regime, where the only characteristic energy scale of the system is the temperature itself and many universal properties emerge rs (). The phase diagram of TlCuCl under pressure and the restoration of classical critical scaling around were the subject of a recent investigation Merchant2014 ().

QPTs in dimerized quantum spin models have been studied numerically by a number of authors, primarily by quantum Monte Carlo (QMC) simulations. Early investigations of the bilayer square-lattice antiferromagnet Sandvik1994 () and other two-dimensional geometries Troyer1996 (); Matsumoto2001 () have been followed more recently by high-precision studies of a range of critical properties Wang2006 (); Meng2008 (); Wenzel2009 (); Sandvik2010 (); Sen2015 (). In three dimensions, the focus of investigations has been on the field-induced transition Nohadani2005 (), on the effects of dimensionality Troyer1997 (); Yao2007 (), and on physical observables at the coupling-induced QPT Jin2012 (); Oitmaa2012 (); Kao2013 ().

A minimal model of a dimerized quantum antiferromagnet has only two coupling constants, on and between the dimer units, and therefore only one control parameter, . The geometry considered in the present study is the double cubic lattice shown in Fig. 1(a). In this system at large , intradimer singlet correlations dominate the physics and the ground state is magnetically disordered, while at small the interdimer correlations establish long-ranged magnetic order. The order parameter of the Néel phase is the staggered magnetization, , and along with the ordering temperature, , it can be driven continuously to zero by increasing , as illustrated in Fig. 1(b). By standard arguments of dimensionality and symmetry, the dynamical exponent of this system is rs (); Chakravarty88 () and the QCP belongs to the O(3) universality class, which is at the upper critical dimension, , of all O() models Zinn-Justin2002 (); rs (). At , mean-field critical scaling behavior alone is not sufficient to capture the physics of fluctuations around the QCP and multiplicative logarithmic corrections to the physical quantities (thermodynamic functions) are expected.

The theoretical importance of multiplicative logarithmic corrections to mean-field scaling behavior lies not only in the statistical physics of condensed matter systems but also in high-energy physics Zinn-Justin2002 (); rkbc (). The general problem of a quantum field theory with an -component field is encapsulated by an “O() theory,” a Lagrangian containing a dynamic (quadratic gradient) term and a potential term with quadratic () and quartic () contributions. On changing the sign of the quadratic term, the system is driven through a QPT separating a phase with from one with (a “Mexican hat” potential). As noted above, the low-energy properties of the 3D dimerized antiferromagnet of SU() quantum spins with Heisenberg interactions correspond to a field theory with and (including the time dimension); and correspond respectively to Ising and XY spin interactions.

Beyond the upper critical dimension , the scaling behavior of the O() theory is straightforward, with the critical exponents being exactly those given by mean-field theory Aizenman1981 (); rbg (); rkbc (), namely , , , , and . As we discuss below, this may be taken as an expression of the independence of quantum and thermal fluctuations when the phase space is sufficiently large. For , the situation is complex and these exponents take anomalous values. However, exactly at the upper critical dimension, , the leading scaling behavior coincides with that of the mean-field theory, but modified by multiplicative logarithmic corrections rkl (); rkle (); Zinn-Justin2002 (). While the leading exponents are -independent, a measure of -dependence resides in the logarithmic corrections and these must be taken into account to establish the universality class of the transition Janssen2004 (). Because the established results for the form of these corrections are based on perturbative techniques applied to low-energy theories, it is desirable to verify them using unbiased numerical methods applied directly to the lattice Hamiltonians, and this is what we achieve here.

Despite the insight into general QPT phenomena obtained from simulations using this type of minimal model for 3D dimerized systems Nohadani2005 (); Troyer1996 (); Yao2007 (); Jin2012 (); Oitmaa2012 (), the question of logarthmic corrections has to date been addressed only briefly and inconclusively tbk07 (); Kao2013 (). Experimentally, the feasibility of observing logarithmic corrections in systems such as TlCuCl remains a challenging open issue Merchant2014 (). In this paper we provide a systematic numerical study. We employ large-scale QMC simulations to investigate the critical behavior of the order parameter and Néel temperature on the double cubic lattice [Fig. 1(a)] for small values of unattainable in all previous studies. State-of-the-art QMC techniques and finite-size-scaling analysis using very large systems (sizes exceeding spins) allow us to detect and characterize the multiplicative logarithmic corrections in the universal scaling relations for the QC regime at the upper critical dimension, here for the O(3) QCP. In fact our results constitute hitherto unavailable exact numerical verification of the logarithmic forms predicted both by perturbative renormalization-group calculations rw (); rwr (); rblzj () based on the -component theory at and by additional considerations exploiting the zeros of the partition function rkl (); rkle (); Kenna2004 (); rkbc (). To the best of our knowledge, no systematic numerical calculations have been performed beyond rkl (); rfkp ().

As will become clear, our numerical results demonstrate to high precision the validity of the detailed theoretical predictions for the expected universality class. Both size-dependent scaling and the order parameter in the thermodynamic limit show evident deviations from pure mean-field behavior, which are accounted for by logarithmic corrections whose exponents are in very close ageement with the predicted values where available. In the case of the Néel temperature, we are not aware of any previous scaling predictions including logarithms. Here we test an Ansatz based on the known scaling behavior for the relevant energy scales rs (); Fisher1989 () and the logarithmic corrections in corresponding classical systems rkbc (). In addition to probing the asymptotic behavior, our results also provide direct insight into the range of validity of logarithmically modified critical scaling forms as one moves away from the QCP, which will be essential in evaluating the experimental relevance of logarithmic corrections.

The paper is organized as follows. In Sec. II we introduce the model and the numerical method, describing the measurement of physical observables in our QMC simulations. In Sec. III we begin the presentation of our numerical results with the precise determination of , the position of the QCP, using finite-size-scaling techniques. Section IV discusses the observation of clear logarithmic corrections in the staggered magnetic susceptibility, , at the QCP as a function of the system size, . We present our results for the sublattice magnetization, , at in Sec. V, discussing in detail its extrapolation to the thermodynamic limit, where we investigate the presence of logarithmic corrections to the leading mean-field behavior. In Sec. VI we present a scaling Ansatz for the Néel temperature, apply finite-size scaling to extract it as a function of , and again investigate corrections to mean-field behavior. We compute the characteristic velocity of spin excitations, demonstrate the precise linearity of and , and discuss the physical interpretation of this behavior. We summarize our results in Sec. VII and comment further on their theoretical and experimental consequences.

## Ii Model and methods

As a representative 3D dimerized lattice with an unfrustrated geometry, we choose to study the double cubic model shown in Fig. 1(a). This system consists of two interpenetrating cubic lattices with the same antiferromagnetic interaction strength, , connected pairwise by another antiferromagnetic interaction, . The QPT occurs when the coupling ratio is increased, changing the ground state from a Néel-ordered phase of finite staggered magnetization to a dimer-singlet (“quantum disordered”) phase, as illustrated in Fig. 1(b). An advantage of this geometry over cases where the dimerization is imposed within a single lattice, such as the simple cubic lattice Jin2012 (), is that all symmetries of the cubic lattices are retained, facilitating the consideration of quantities such as the spin stiffness or the velocity of spin excitations.

The Hamiltonian is given by

 H=∑⟨i,j⟩Jij→Si⋅→Sj, (1)

where is an spin operator residing on a double cubic lattice of sites with periodic boundary conditions. The sum is taken only over nearest-neighbor sites, where every site has six neighbors on the same cubic lattice with coupling strength and one neighbor on the opposite cubic lattice with [Fig. 1(a)]. We set as the unit of energy and use as the control parameter.

To study this system, we use the stochastic series expansion (SSE) QMC technique Sandvik1991 (); Sandvik1999 (); Evertz2003 (); Sandvik2010 () to obtain unbiased results, i.e. numerically exact within well-characterized statistical errors, for physical quantities in systems of finite side-length . Here we present results up to at temperatures with up to . We then perform detailed analyses by finite-size scaling Fisher1972 () to extract information in the thermodynamic limit both in the ordered state and at the QCP, as detailed in the separate sections to follow. Here we define the physical quantities of interest and discuss some technical aspects of their calculation within the SSE QMC method.

Because spin-rotation symmetry is not broken in simulations of finite-size systems, one may measure the squared order parameter and take its square root as a post-simulation step. The staggered magnetization is given by

 ms=√1NS(QAF), (2)

where

 S(q)=1NN∑i,je−iq⋅(ri−rj)⟨→Si⋅→Sj⟩ (3)

is the magnetic structure factor, with denoting the real-space position of the spin on lattice site , and is the wave vector of antiferromagnetic order, with the fourth denoting the phase factor between the two cubic lattices. We consider only the -component of the magnetization and average it over the time dimension of the QMC configurations, computing the expectation value of the squared order parameter in the form

 m2sz=1β∫β0dτm2sz(τ), (4)

where

 msz(τ)=1NN∑i=1e−iQAF⋅riSzi(τ), (5)

with

 Szi(τ)=eτHSzie−τH (6)

the time-evolved spin operator at imaginary time . In an SSE simulation, the integral in Eq. (4) is transformed into a discrete sum with no approximations and the relation compensating for the rotational averaging of the single measured component of the order parameter,

 ms=√3⟨m2sz⟩, (7)

is applied post-simulation. The magnetic susceptibility is defined as

 χ(q)=1N∑ij∫β0dτ⟨Szi(τ)Szj(0)⟩e−iq⋅(ri−rj). (8)

In the SSE approach, the squared order parameter [Eq. (4)] is readily evaluated at any because the QMC configurations are constructed in the basis Sandvik1997 (); Sandvik2010 (). The dynamical spin-spin correlation function contained in Eq. (8) can also be obtained easily by applying an operator string connecting the states at different imaginary times, with the integral computed analytically to give a direct, formally exact QMC estimator not requiring post-simulation integration Sandvik1999 (); Sandvik2010 ().

The Binder ratio Binder1981 () is the ratio of the fourth moment of a quantity to the square of its second moment. For our purposes, the relevant quantity is

 R2=⟨m4sz⟩⟨m2sz⟩2, (9)

which is dimensionless and satisfies the crucial property of being size-independent at the QCP in the limit of large system sizes. The spin stiffness, or helicity modulus, of the system is defined as

 ραs=1N∂2F(ϕα)∂2ϕα∣∣ϕα→0,α=x,y,z, (10)

where is the free energy and is the angle of a twist imposed between all spins in planes perpendicular to the axis. In an SSE simulation, the most efficient way to extract the spin stiffness is to take the derivative in Eq. (10) directly in the QMC expression for at , giving

 ραs=3⟨w2α⟩4β, (11)

where

 wα=1L(N+α−N−α) (12)

is the winding number Pollock1987 (); Sandvik1997 () of the spin in spatial direction and and are the numbers of occurrences of the operators and on bond in the -direction within imaginary time . As noted above, is the same in all three directions due to the symmetry of the double cubic lattice, and the average may be taken over all of these. The spin stiffness follows the scaling law in spatial dimensions Sandvik2010 () and, because the dynamic exponent is here, the quantity , or equivalently , is also size-independent at the QCP, up to a logarithmic correction at the upper critical dimension.

Finally, the spin-wave velocity, , can be obtained reliably by monitoring the fluctuations of the spatial and temporal winding numbers Kaul2008 (); Jiang2011 (); Kao2013 (); Sen2015 (). For a fixed system size, , the inverse temperature is adjusted to the value where the system has equal winding-number fluctuations in the spatial and temporal directions,

 ⟨w2α(β∗)⟩=⟨w2τ(β∗)⟩,α=x,y,z. (13)

The temporal winding number is the net magnetization (number of up minus number of down spins) of the system, , which is easily obtained in the basis Sen2015 (). When the condition (13) is met, the spin-wave velocity is given by the ratio

 c=Lβ∗(L). (14)

The isotropy of the lattice is an advantage also in this case. For each value of , an extrapolation is performed to obtain in the thermodynamic limit. For further details of these procedures, we refer the reader to the recent extensive tests of this method conducted in Ref. Sen2015 ().

## Iii Determination of the QCP

The key to an accurate characterization of logarithmic corrections is a high-precision determination of the location of the QCP. For this we employ the Binder ratio, , and the appropriately scaled spin stiffness, , which both have scaling dimension zero and therefore should approach constant values at when , up to possible logarithmic corrections. We stress that the scaling forms for the approach of both quantities to the critical point are valid for a four-dimensional (4D) theory, with the temperature (imaginary time) providing the fourth dimension, and are applicable on a “critical contour” where the inverse temperature is taken to infinity symmetrically with the spatial dimension of the system. This form is appropriate for a system with dynamic exponent (in general . All values of yield the same results in the limit , and in principle the contour is optimal when ; the spin-wave velocity, , is a number of order unity discussed in detail in Sec. VI, and here we use .

Away from , and approach different constant values with increasing system size. In the Néel state, approaches due to diminishing fluctuations in the magnitude of the rotationally-invariant order parameter, of which we measure only the -component in Eq. (9). In the quantum disordered phase, approaches a higher value dictated by Gaussian fluctuations, which from the symmetries of the double cubic model is . The spin stiffness falls from non-zero values in the Néel phase to zero in the disordered phase. When calculated as functions of , the curves and obtained for different system sizes should cross at the QCP, up to corrections that are well understood from the theory of finite-size scaling. We analyze these corrections to obtain an unbiased value of the critical coupling, , in the thermodynamic limit Sandvik2010 ().

Figure 2(a) shows as a function of in the neighborhood of the critical coupling ratio for various system sizes. We have performed simulations for all even-length sizes , 8, 10, …, 40, but here we present only the , 32, …, 40 data for clarity. Analogous curves for the scaled spin stiffness, , are shown in Fig. 2(b), again only for system sizes , 32, …, 40. In both cases, the system sizes are sufficiently large that the crossing points exhibit only a very weak dependence on on the scale used in the figure, and both sets of data may be used independently to show that the QCP is located at . A detailed analysis is required to obtain the most precise results attainable, free of any finite-size effects, and we first discuss the general scaling behavior of before presenting our numerical results.

### iii.1 Scaling Forms for Critical-Point Estimators

To describe the evolution of the crossing points with , we perform a systematic extrapolation of the finite-size data to the thermodynamic limit by extracting the crossing points between data sets for all pairs of system sizes, and , based on polynomial interpolations. Figure 2(c) shows the crossing points and obtained in this manner.

For any quantity probing a singularity in the thermodynamic limit, one may define a size-dependent critical point . In general, this quantity is expected to shift by an amount proportional to with respect to the true infinite-size critical point, (hereafter denoted for simplicity by ), i.e. for large ,

 gc(L)=gc+aL−1/ν, (15)

where is the standard correlation-length exponent. However, with a definition based on crossing points of a dimensionless quantity computed for two different sizes, the leading corrections cancel and the convergence is faster,

 gc(L)=gc+aL−(1/ν+ω), (16)

where is the dominant irrelevant exponent. In practice, with data fits to a rather limited range of available system sizes, the corrections to Eq. (15) contained in Eq. (16) will have exponents and prefactors that deviate from their asymptotic values due to the neglected corrections of higher order in and therefore these should be considered as “effective” quantities.

The above forms are applicable in the absence of logarithmic corrections, but such corrections are the primary focus of our study and are expected at . Kenna has derived the modified form of Eq. (15) for classical systems with logarithmic corrections rkbc (); rkb (),

 gc(L)=gc+aL−1/νln^λL, (17)

where the exponent of the logarithm for the 4D O universality class is . For the crossing points, Eq. (16) is modified to

 gc(L)=gc+aL−(1/ν+ω)ln^cL, (18)

as shown in Appendix A, where if the subleading term has no multiplicative logarithmic correction, but is altered by an unknown amount if it does. Under the circumstances, with a number of unknowns and with simulation data only for a restricted range of system sizes, we fit our data using not Eq. (18) but instead the purely algebraic form of Eq. (16) with and , the effective value of the subleading exponent over the fitting range, treated as a different fitting parameter for the separate quantities and .

### iii.2 Numerical Determination of gc

We take both pairs of fitting parameters and in Eq. (16) to be independently free for the two datasets and , but impose the constraint that the curves have the same . As shown in Figs. 2(a) and 2(b), we obtain good fits to both functions, meaning with a reduced (per degree of freedom, hereafter denoted ) close to , to the data for all system sizes (). These allow us to conclude that , where the numbers in parenthesis denote the expected errors (one standard deviation) in the preceding digit, i.e. the relative error on is approximately one part in . If we allow independent parameters and in the fits to the two data sets, both estimates of the critical point are statistically consistent with this , albeit with somewhat larger error bars.

We note here that the values we find for the subleading exponent, for the data and for the data, lie far from a common asymptotic value. Thus indeed should be considered as an effective exponent accounting for crossover effects in system size, neglected higher-order irrelevant fields, and the expected weak logarithmic corrections. However, the good match obtained between the two extrapolated -estimators, especially when approaching the infinite-size value from different directions, would not be expected in the presence of any corrections not taken sufficiently into account by the fitting functions. Thus we believe the error bar on quoted above to be completely representative of all statistical and systematic uncertainties, in the sense that any remaining systematic errors due to the fitting form should be smaller than the statistical errors. The tests we perform on the critical scaling behavior around the QCP in the subsequent sections also support this statement.

## Iv Size-Dependent Logarithmic Corrections at the QCP

The critical O() theory, by which is meant the theory at the upper critical dimension and at the critical point, obeys many fundamental and universal properties, some of which depend on while others are -independent. In Ref. Kenna2004 () it was shown that the zeros of the partition function (Lee-Yang zeros) YangLee1952 (), and hence the thermodynamic functions, obey a finite-size scaling theory, which was derived by renormalization-group methods. These perturbative arguments demonstrate that there are multiplicative logarithmic corrections in the system-size dependence of derivable thermodynamic functions, which are closely linked to those of the Lee-Yang zeros and, furthermore, are independent of for odd . This leads to the key practical observation that size-dependent logarithmic corrections in physical observables such as the magnetic susceptibility and the specific heat at the critical point follow a universal, -independent form when . Here we provide a non-perturbative calculation of the magnetic susceptibility, [Eq. (8)], for systems of finite at the QCP, , of the -dimensional O() transition to test the predicted logarithmic corrections.

The universal form of the magnetic susceptibility at the critical point in a finite-size system is given by Kenna2004 ()

 χ(QAF,L)=aL2[lnL]1/2[1+bln(lnL)lnL], (19)

with non-universal but -independent parameters and . We used this expression with a fixed value , which is within the standard deviation of the value found in Sec. III, to investigate the logarithmic corrections to the -dependence of . We calculate the susceptibility at the ordering wave vector, , for systems of all even sizes from to 40. As in Sec. III, the scaling predictions under test are valid for a 4D theory and again we use the critical contour with .

Figure 3(a) shows our results for the critical magnetic susceptibility normalized by . In the absence of logarithmic corrections, would be constant and the curve would be a flat line. Instead we observe that a reasonable account of the data for our larger system sizes () requires a fit of the form , as anticipated in Ref. Kenna2004 (). However, for a more quantitative fit over the full size range available, we find [Fig. 3(a)] that it is necessary to include the predicted subleading logarithmic correction in Eq. (19).

To examine the sensitivity of these results to the exponent of the multiplicative logarithm in Eq. (19), we replace this predicted exponent by a variable . We determine this exponent by calculating the goodness of fit as a function of . As Fig. 3(b) makes clear, the best fits are indeed obtained close to , in complete consistency with Eq. (19).

From the fact that our exact numerical data confirm not only the leading but also the subleading corrections to scaling, we conclude that obvious logarithmic corrections can be observed in the size-dependence of the thermodynamic functions at the QCP. This result also demonstrates that our determination of is sufficiently precise to study logarithmic corrections without significant distortions arising from uncertainties in its value.

## V Sublattice Magnetization

Physical condensed-matter systems at continuous QPTs are generally in the thermodynamic limit, and size-scaling measurements of the type easily performed in QMC simulations (Sec. IV) are not a realistic experimental option. However, as discussed in Sec. I, multiplicative logarithmic corrections are expected in a range of physical quantities close to the QCP. The primary physical observables in the quantum antiferromagnet are the zero-temperature staggered magnetization, , and the Néel temperature, . Calculating these quantities in the thermodynamic limit is significantly more challenging than studies of size-dependence, as careful extrapolations of finite-size data are required. Here and in Sec. VI we describe and then implement appropriate measures for extrapolating to infinite system size and, for , to zero temperature, thereby revealing the logarithmic corrections to both and .

We compute the staggered magnetization according to Eq. (7) for coupling ratios as close to as . For a given value of , we calculate the squared quantity over a range of system sizes. As shown in Fig. 4, the staggered magnetization clearly decreases with increasing and converges towards a fixed limit, suggesting a controlled extrapolation even very close to the QCP. Definitive extrapolation to the thermodynamic limit in this regime is a complex issue, and a discussion of several technical points is in order before analyzing our results.

### v.1 Extrapolation scheme

First, most of the simulations in this section are performed at a temperature , such that the extrapolation includes both system size and temperature. Because the system is ordered for sufficiently low and the order parameter converges quickly to a non-zero value below the ordering temperature, one may use the form with arbitrary prefactor and exponent to study the magnetization as a function of . These conditions for obtaining a () extrapolation of the order parameter for contrast with the need to follow the contour (; when ) in Secs. III and IV for studies of the QCP. Although large values of and should improve the convergence, in practice one must consider the balance between computation time and convergence rate, and the choice , works well in most cases. However, for coupling ratios very close to the QCP, the temperature may be a significant fraction of and thus could be far from its zero-temperature value. We have therefore performed additional simulations at to verify that the extrapolation does remain well controlled and fully representative of the thermodynamic limit in temperature as well as in system size.

Second, in contrast to Sec. III, where we used the non-trivial power-law scaling forms (16) known to be appropriate for extrapolating the location of a critical point, the ground-state order parameter inside the Néel phase can be extrapolated by using simple polynomial fits. To obtain , one may extrapolate the squared quantity and then take its root afterwards or take the square root for each system size before extrapolating (the procedure followed in Fig. 4); the corresponding polynomials are

 m2sz(g,L) = a(g)+b(g)L−2+c(g)L−3+…, (20) √m2sz(g,L) = a′(g)+b′(g)L−2+c′(g)L−3+… (21)

Here the leading -dependence in the extrapolation of a non-vanishing order parameter, at fixed inside the ordered phase, is known cardy () to be due to the dimension-dependent power-law decay of the transverse correlation function. Here , and the resulting leading dependence is shown clearly in Fig. 4. Because these procedures use a polynomial of finite order to approximate physical behavior containing, in principle, an infinite number of corrections, the extrapolated value of obtained with the forms (20) and (21) will not be exactly the same, but for reliable fits they should agree within statistical errors.

We stress that no logarithmic corrections are expected in this case, meaning that on grounds of principle they should not be present in the asymptotic large- corrections to a non-zero-valued order parameter. This non-critical behavior contrasts with the case of the shift in the critical point discussed in Sec. III, where logarithmic corrections should in principle be present, although we concluded that their effects are not detectable in practice.

Finally, however, non-trivial corrections may still be expected in the dependence of for close to the QCP, where the order parameter is small, in the form of crossover behavior from near-critical at small system sizes to asymptotic ordered-state scaling at large . Quite generally, no analytic functional forms are available for describing such crossovers, and great care is required to ensure that the asymptotic region, where Eqs. (20) and (21) are valid, has been reached. As , successively larger systems are required for this, and here we find that reliable extrapolations are no longer possible beyond because of the limits on system size set by the available computer resources.

In fits to the forms (20) and (21), it is necessary to select the order of the polynomial and the range of system sizes to include. The size of the error bars on the QMC data points has a significant influence here, because deviations from the leading correction are easier to detect with smaller error bars. We have performed a systematic study using fits of orders to 6, including different system-size ranges. We characterize the quality of the fits using the standard reduced measure, and for a “good” fit we require that the optimal value must fall within three standard deviations of its mean, i.e. we demand that

 χ2r−1=χ2nL−np−1≤3√2nL−np, (22)

where is the number of data points (system sizes) and is the number of fitting parameters. For a given and largest system size , we use all available system sizes down to a smallest size for which the above condition is still satisfied. We then study the behavior as a function of for different , and compare the values of obtained from extrapolations based on Eqs. (20) and (21). To estimate the error bars on the extrapolated , we performed additional polynomial fits with Gaussian noise (whose standard deviation is equal to the corresponding QMC error bars) added to the finite-size data. The standard deviation of the distribution of extrapolated values defines the statistical error.

Figure 5 shows results for several values of approaching the QCP, with cases where the fits are relatively straightforward (further from ) shown in panels (a) and (b) and more challenging cases (closer to ) shown in panels (c) and (d). The upper panel in each group corresponds to the square root being taken after the extrapolation [Eq. (20)], while the lower corresponds to fitting the square root for each system size [Eq. (21)]. In panels (a) and (b), the extrapolated values are observed to be very stable with respect to the range of system sizes and the order of the polynomial, whereas panels (c) and (d) manifest some of the crossover behavior expected close to , showing considerable variation as the maximum system size is increased. There are also significant differences between the two fitting procedures, until the largest system sizes where the extrapolations stabilize; we take the fact that the two types of fits give consistent results for these largest systems at all values of as an indication that the extrapolations are reliable. We have not been able to achieve good convergence based on system sizes up to for values closer to than those shown in Figs. 5(c) and (d). In these most challenging cases, our results show that it is better to use the fitting form of Eq. (21), extrapolating the square root of the staggered magnetization for each system. Regarding the quality of the fits obtained by varying the polynomial order, , in Fig. 5, we find that extra terms in the fit scarcely justify the additional degrees of freedom lost in the determination of . All of the results presented below were obtained by extrapolating with polynomials of order .

### v.2 Thermodynamic Limit

With all of the above considerations, we are able to obtain reliable and high-precision extrapolations of the staggered magnetization in the thermodynamic limit for values of as close to the QCP as . In Fig. 6 we show all of our data for on logarithmic axes. If these data satisfied mean-field scaling alone, with no discernible logarithmic corrections, one would expect a curve of the form

 ms(g)=a|g−gc|1/2, (23)

but this (green line in Fig. 6) is manifestly unable to describe the data. For the zero-temperature order parameter, perturbative renormalization-group considerations applied to the O() field theory at the upper critical dimension predict the form

 ms(g)=a|g−gc|β|ln(|g−gc|/b)|^β, (24)

where is the mean-field exponent and the exponent of the multiplicative logarithmic correction is given by Kenna2004 (). A fit to this form, using for (blue curve in Fig. 6), yields excellent agreement with the data all the way to our smallest values of ; the fitting parameters are and . We note that the fit is very insensitive to the precise value of , and for further analysis we fix this to .

To test the predicted exponent in Eq. (24), we treat it as a free parameter and fit our data using different numbers of values, including all points closest to and studying the behavior as points further away from are added one by one. Figure 7 shows and as functions of the number of data points fitted. With the exception of cases including the two points furthest away from , all the fits appear reasonable, with . However, by the properties of the distribution, a fit should be considered statistically acceptable only if a criterion analogous to Eq. (22) is satisfied, i.e. the largest number of data points for which remains less than three times its standard deviation () marks the boundary between good and poor fits. At this point we obtain , which lies well within one standard deviation of the predicted value . If more points are excluded, the fitted exponent evolves slowly [Fig. 7(b)] while remaining statistically well compatible with the predicted value. Because the fitting error increases, less weight should be placed on results including less data, and taking an error-weighted average over all the points below the cut-off line, , in Fig. 7 yields . We take this as complete confirmation of the predicted value.

As important as finding clear logarithmic corrections to scaling is that we have demonstrated their presence over a significant region around the QCP; indeed, most of the points we have computed are well described by Eq. (24). Including the multiplicative logarithmic correction converts an inadequate description of the data into an excellent one (Fig. 6) as far inside the Néel phase as , where the order parameter is already at of its maximum possible value (, at which point no quantum fluctuation effects remain). This improvement is clearer still in the inset of Fig. 6, which shows the results on linear axes. Under the assumption that data points at large no longer fall on the fitted curve because they lie outside the region controlled by the QCP, we can determine the size of the critical region based on a threshold maximum deviation of the data from the curve. Although the choice of threshold value is somewhat arbitrary, the region indicated by the yellow shading in Fig. 4 reflects a threshold of approximately , which lies well above achievable experimental uncertainties. We comment in Sec. VII on the utility of our results for the case of TlCuCl.

## Vi Néel temperature

We turn next to the scaling form of the Néel temperature, , near the QCP. Unlike the order parameter, as far as we are aware there is no prediction from perturbative field-theoretical calculations including logarithmic corrections for the scaling form of finite- critical points at the upper critical dimension. Close to a QCP, the general power-law form without logarithmic corrections is discussed in Ref. rs (), but one may also expect a multiplicative logarithmic term as in the other quantities we have discussed. We first derive the exponent of the logarithm for the O() transition in 3+1 dimensions, based on the known scaling properties of related quantities. We then present our QMC calculations of for the Heisenberg model on the double cubic lattice and test our prediction.

### vi.1 Scaling hypothesis

In a path-integral construction in imaginary time, the size of the system in the time dimension is proportional to the inverse temperature, . This can be considered as a length, , on a parallel with the spatial lengths, , of a -dimensional system. If the spatial lengths are taken to infinity in all directions, what remains is a single finite length, , for the effective system, and finite-size scaling in this length corresponds to finite- scaling in the original quantum system Chakravarty88 ().

Without logarithmic corrections, by analogy with the finite- shift of the critical point discussed in Sec. III.1, the same type of shift as in Eq. (15) can be expected because . Thus

 gc(T)−gc(0)∼L−1/ντ, (25)

as a consequence of the finite temporal size, and the scaling behavior is , as discussed in detail in Ref. rs (). In the case of spatial finite-size scaling, with all lengths finite, the shifted critical point (sometimes called the pseudo-critical point) is not a singular point, but the singularity develops as . By contrast, in the finite- case in spatial dimensions, the shifted point is a true (classical) phase transition, although from a scaling perspective this difference is not relevant.

In order to discuss logarithmic corrections, it is useful to first express using a macroscopic, zero-temperature energy scale of the system that vanishes as rs (). For the spin system considered here, the only such energy scale is the spin stiffness, . According to Ref. Fisher1989 (), the scaling form of this quantity in the ordered phase when is

 ρs∼(gc−g)ν(d−1). (26)

Consistency with the result then gives the scaling of the critical temperature for ,

 T2N∼ρs, (27)

where the mismatch in units is compensated by a power of the non-singular spin-wave velocity rs () (Sec. VIC), which can be neglected here. Our basic hypothesis is that this proportionality, which is the singular part of a relationship based on matching scaling dimensions, applies in all respects at the upper spatial critical dimension ( for ), such that logarithmic corrections to arise solely due to the logarithmic corrections intrinsic to .

Fisher et al. Fisher1989 () have shown that the critical spin stiffness can be expressed as , where is the correlation length and is the free-energy density. The logarithmic corrections to both and , presented by Kenna in Ref. rkbc (), are

 ξ∼|g−gc|−νln^ν(|g−gc|), (28)

with for the relevant universality class, and

 f∼|g−gc|4νln^α(|g−gc|), (29)

with . The logarithmic correction to is therefore given by

 ρs∼|g−gc|2^νln2^ν+^α(|g−gc|), (30)

and by combining these results with Eq. (27) we obtain

 TN∼|g−gc|νln^τ(|g−gc|), (31)

where . From the values of and given above rkbc (), we obtain the prediction , which is remarkable in that the exponent in the logarithmic correction to should be the same as the one in the sublattice magnetization, (24).

Because the zero-temperature order parameter is a consequence purely of quantum fluctuations, whereas the classical ordering temperature is a consequence primarily of thermal fluctuations, there is a priori no reason to expect that the two should have the same form. Exact numerical calculations are therefore uniquely positioned to provide qualitatively new information in this case. We note that this equality applies to the phase transitions of O() models for all values of ; because and rkbc (), we obtain , the same value as the exponent in Eq. (24). Thus we predict that, in the neighborhood of , will be proportional to , with no multiplicative logarithmic factors, for all values of ; this result was reported for the case in a previous QMC study Jin2012 (), which we now extend sufficiently close to the QCP to observe the cancellation of logarithmic terms.

### vi.2 QMC calculations

Calculating within our QMC simulations is similar to obtaining in Sec. III, but with some important differences of detail. The calculations of Sec. III were performed for a genuinely 4D system, with the imaginary-time axis treated on the same footing as the spatial dimensions. At finite temperatures, this symmetry is broken and the system is 3D with a separate temperature variable, which determines the finite thickness of the time dimension even when . Both the Binder ratio [Eq. (9)] and the spin stiffness [Eq. (11)] are size-independent quantities at the thermal phase transition and hence remain valuable indicators, although the appropriately scaled spin stiffness for the 3D transition is now (instead of , used for analyzing the 4D transition in Sec. III) Jin2012 ().

For each value of the coupling ratio within the Néel phase (), we compute and for a range of system sizes and perform finite-size-scaling extrapolations to deduce the Néel temperature, , in the thermodynamic limit. Similar to Sec. III, we first obtain the crossings of the and data for different system sizes using polynomial fits, as shown in Figs. 8(a) and 8(b) for . The crossing points of both quantities for each successive pair of system sizes, and , are used to extrapolate towards the value from above and below, using power-law forms analogous to Eq. (16). We note that, as in the analysis leading to (Sec. III.2), the data points obtained for and using systems of all sizes () fall within a criterion analogous to Eq. (22) for this type of fit. The extrapolation of for is shown in Fig. 8(c).

We comment here that our determination of for close to is rather less precise than our determination of . The fundamental difference in character of the two quantities, and hence of their calculation, causes the estimators for [the approximate crossings in Figs. 8(a) and 8(b)] to have larger error bars and finite-size effects. Further, the error bars of the crossing points grow rapidly as , while the decrease in leads to longer simulation times (because the space-time volume is proportional to ). After detailed error control, the closest reliable data point to the QCP is , for which is twice as large as for the closest point (Sec. V). We have nevertheless obtained 19 reliable data points, within the QC regime determined from (Fig. 6) and down to unprecedentedly low temperatures, which are fully sufficient to test for evidence of logarithmic corrections to .

The Néel temperature has units of energy and clearly depends on the overall energy scale of the system. Ideally, it should be normalized by an intrinsic energy scale of the system to give a dimensionless quantity. In Ref. Jin2012 () it was shown that normalized to the microscopic energy scale , given by the sum of all couplings of a spin to its neighbors, yields a remarkably system-independent result for Heisenberg antiferromagnets with three different dimerization patterns; for the double cubic lattice, . Other authors Oitmaa2012 (); Kao2013 () have suggested that the appropriate normalization is given by a macroscopic quantity, the spatially averaged spin-wave velocity, (which, it should be noted, does not have units of energy and requires an unknown dimensionful constant). We begin by taking the former approach and return below to address the latter.

The relationship between and is presented in Fig. 9. Once again we show a mean-field scaling line for comparison and once again it cannot provide an adequate fit, suggesting that logarithmic corrections are indeed present. However, a fit to our predicted form, given by Eq. (31) with , describes the data very well, even at the limits of the region classified as QC based on the fit in Fig. 5.

For a fully quantitative test of the exponent we predict for the multiplicative logarithmic term in Eq. (31), we substitute a free exponent for the fixed value and optimize it using fits with different windows of -values. This analysis is precisely analogous to that performed for in Fig. 7. The behavior of and of the optimized exponent, with error bars again estimated using the method of numerical Gaussian noise propagation, is presented in Fig. 10. By taking the inverse-variance-weighted average over all results for which is acceptable, we obtain , in excellent agreement with the prediction . We conclude that the multiplicative logarithmic correction to is, to within our error bars and in agreement with a straightforward scaling argument based on the spin stiffness (Sec. VIA), identical to the correction in Eq. (24).

### vi.3 Spin-wave velocity

The spin-wave velocity, , is uniform in the primary axial directions on the double cubic lattice. As discussed in Sec. II, it can be calculated most straightforwardly and most accurately in the SSE framework from the spatial and temporal winding-number fluctuations in Eq. (13) to define the space-time-isotropic criterion of Eq. (14), which contains the velocity Kaul2008 (); Jiang2011 (); Kao2013 (); Sen2015 (). This technique remains well-defined throughout the critical regime and is expected not to be affected by logarithmic corrections; it was shown in Ref. Sen2015 (), which we follow for technical details, that the winding-number approach produces the correct result for in the Heisenberg chain, a system known to have strong logarithmic corrections to scaling.

In Fig. 11(a) we show the results for of calculations on finite systems of even sizes up to , which we extrapolate to the thermodynamic limit using the relation

 c(g,L)=c(g)+a(g)/L2+b(g)/L3. (32)

This form is found empirically Sen2015 () to provide a very good reproduction of the data and numerical errors due to finite-size effects in the critical regime are clearly small. Figure 11(b) shows the results for the extrapolated spin-wave velocities of the infinite system, by comparing the macroscopic scale with the microscopic quantity discussed above. The almost perfect linearity demonstrates that the two effective energy scales are very closely related, which can be expected from the fact that the velocity of a spin excitation depends directly on the net interaction of a single spin, and both are perfectly valid choices for the normalization of . A graph completely analogous to Fig. 9, showing a logarithmic correction with the same exponent , is obtained if is used to normalize .

### vi.4 Relation between Tn and ms

In experiment it is often difficult to relate an external control parameter to the microscopic coupling constants of a model Hamiltonian. In a quantum antiferromagnet, some aspects of this problem can be circumvented by studying the relationship between and directly, without reference to the control parameter, . A universal relationship between these macroscopic and measurable quantites would be of considerable experimental utility in characterizing the nature of critical phenomena without recourse to detailed microscopic knowledge of the system parameters (such as the pressure dependence of the exchange couplings in TlCuCl). Although an experimental test Merchant2014 () of the linear relationship between and Jin2012 () indicated satisfactory agreement close to the QCP, the issue of how best to normalize was not addressed. We use our systematic data spanning the entire QC regime to test the limits of linear proportionality and discuss the normalization of .

In Ref. Jin2012 (), where a universal linear relation was found in three different models, the authors articulate a mean-field argument based on semiclassical considerations for a direct proportionality of to the effective spin order gauged by at . In Ref. Oitmaa2012 (), these arguments were elucidated in a field-theory context, where it was stated that logarithmic corrections should be negligible for linear proportionality to emerge. In fact these arguments can be reduced to the statement that it should be possible to treat quantum and thermal fluctuations independently, with no mutual interference of their effects Jin2012 (). If one considers that mean-field exponents are valid in high-dimensional systems because thermal fluctuations become independent of quantum fluctuations when the phase space is sufficiently large, then it appears that weak logarithmic corrections could enter the relationship of to at . This possibility, also motivated by the (then) unknown form of the logarithmic corrections to , was investigated directly by QMC simulations for the cubic lattice Kao2013 (), but the results were not conclusive (claims concerning the observation of logarithmic corrections are not justified by the available data range). Here we have presented scaling arguments (Sec. VIA) and numerical data (Sec. VIB) demonstrating that the logarithmic corrections to and have precisely the same form, setting their linear relationship in this class of system on a far firmer foundation.

Our data from Secs. V and VI can be used to probe the relation in detail and confirm that linearity extends much closer to the QCP than previous studies could show. Because our results are for a single type of dimerized model, we are not able to address the question of a universal prefactor Jin2012 (). However, we are able to make a definite statement regarding logarithmic corrections in the relationship between and . Figure 12 shows our data, taken from Figs. 6 and 9, in the form , with the implicit control parameter effectively eliminated. In panel (a), is simply normalized by the energy scale, the interdimer coupling ; in panel (b), we have normalized by the composite scale (where can be considered a function of ), as in Fig. 9; in panel (c), we have normalized using the correctly scaled spin-wave velocity, . The shaded regions again signify our definition of the critical region, based on the strict critical scaling of in Fig. 5.

As shown in Ref. Jin2012 (), we find in Fig. 12(b) that the linearity of and extends well beyond the QC regime; although our data were not selected to focus on this region, we do find complete agreement with previous calculations where the data overlap. Here we demonstrate that the essentially perfect linearity also extends much closer to (), and in particular that it remains valid throughout a regime with explicit logarithmic corrections in the individual quantities. Although we cannot show that the linear relationship extends all the way to the QCP, our data certainly suggest that this is the case, i.e. that it can also be considered a universal property of the QC regime. To the extent that linearity of and is a consequence of the decoupling of thermal and quantum fluctuations, this independence appears to extend from strongly ordered systems, where no logarithmic corrections are expected, to the most strongly fluctuating QC systems. The linear relationship we have demonstrated verifies in full our scaling prediction that the logarithmic corrections to have the same exponent, , as . We also note that the normalization of has no effect (beyond the prefactor) on the linear relationship in the QC regime, but that deviations from linearity clearly differ at higher .

We close by considering in more detail the case where is normalized by [Fig. 12(c)], with a view to making quantitative comparisons with field-theory predictions Oitmaa2012 (). The explicit relationship is

 TN=γc3/2√125ms, (33)

where is the dimensionful prefactor relating the expectation value of the un-normalized field in the action to the order parameter of the lattice model. In the dimensionless units of our work (, …), we obtain , thereby providing a bridge between the quantum field theory and the microscopic lattice Hamiltonian. We suggest that this calculation should be repeated for other dimerized geometries to test the universality of . It is worth repeating in this context the advantages of the double cubic lattice in making the spin-wave velocity equal in all three primary axial directions. The winding-number method used to extract in Sec. VI.3 can be generalized to anisotropic systems jiang11b (), but incurs the significant complication of altering the aspect ratio of the spatial lattice. Different techniques for computing the velocities, such as those based on the hydrodynamic relationship among , the spin stiffness, and the magnetic susceptibility Sen2015 (), may then be more convenient in practice.

### vi.5 Width of the Classical Critical Region

A key question raised by the experiments on TlCuCl Merchant2014 () concerns the width of the region close to where classical critical scaling applies. It was found that this width, , is essentially constant when normalized by . We employ scaling arguments to show that the normalized width, , should indeed be a constant with only a weak logarithmic correction in 3+1 dimensions.

For fixed temperature , we consider the correlation length, defined in terms of the approach of the spin-spin correlation function to its asymptotic long-range value, , when approaching the critical coupling ratio from the ordered side. This quantity has an initial divergence governed by the 4D QCP,

 ξ(T)=ξ4(T)∼[gc(T=0)−g(T)]−ν4, (34)

with mean-field exponent , because the temporal thickness far exceeds and the system cannot sense its finite temporal extent, behaving as at . At the point where reaches , the behavior crosses over to a 3D scaling form,

 ξ(T)=ξ3(T)∼[gc(T)−g(T)]−ν3, (35)

where is the 3D O() exponent. Without logarithmic corrections, the temporal length is simply (more precisely, ), but at the upper critical dimension this relationship is modified by a logarithmic factor,

 Lτ∼|ln(T)|^q/T, (36)

which is obtained by generalizing the classical result of Kenna rkbc (). For the 4D O() universality class, rkbc (); rkb (), and the correlation length itself also has a logarithmic correction,

 ξ4∼(gc−g)−ν4|ln(gc−g)|^ν, (37)

with . The quantum-classical crossover taking place when therefore corresponds to

 |ln(T)|^q/T∼[gc−g)]−ν4|ln(gc−g)|^ν, (38)

which, by converting to a temperature-dependence and keeping only the leading logarithm, yields the crossover temperature

 T∗(g)∼(gc−g)1/2|ln(gc−g)|^q−^ν. (39)

Using our result for [Eq. (31)], the width of the classical critical region on the ordered side of the transition is therefore

 ~W(g)=TN(g)−T∗(g)TN(g)∼1−a|ln(gc−g)|(^q−^ν)/^τ, (40)

with a constant , whose calculation requires further considerations, and a small exponent

 ^q−^ν^τ=1/4−5/223/11=1/12, (41)

on the logarithm. This very weak dependence explains the near-constant behavior found for TlCuCl Merchant2014 (). On general grounds we expect the width of the classical critical regime on the other side of the transition to scale in the same way.

## Vii Summary

We have provided a direct and non-perturbative verification of the existence and nature of multiplicative logarithmic corrections to scaling at the quantum phase transition for three-dimensional dimerized quantum Heisenberg antiferromagnets. These systems correspond to the field theory of an O() quantum field in dimensions, which is the upper critical dimension () for all models with O() universality. With the exception of the Ising model () Kenna1994 (), no such demonstration exists to date, despite a significant body of analytical and numerical work on quantum criticality in dimerized quantum antiferromagnets.

Our results are obtained from large-scale quantum Monte Carlo calculations based on state-of-the-art simulation techniques and detailed finite-size-scaling analysis. These enabled us to extract the precise logarithmic corrections to the leading critical properties at the quantum phase transition from a non-magnetic state of dominant dimer correlations to a Néel-ordered antiferromagnetic state. Specifically, we have obtained the multiplicative logarithmic corrections to the mean-field behavior of the order parameter, the zero-temperature staggered magnetization (), on the control parameter, the coupling ratio . We have verified that these are governed by an exponent , a value we specify with numerical (statistical) precision under , matching precisely the prediction of perturbative renormalization-group calculations Zinn-Justin2002 (); rkbc ().

No prediction was previously available for the analogous logarithmic correction to the Néel temperature, . We have implemented a scaling Ansatz exploiting the known logarithmic corrections of other physical quantities to obtain its form. Our prediction is that has exactly the same exponent in its logarithmic correction, , as the order parameter, and our numerical results for are in excellent agreement. We have thereby established an exact linearity between and throughout the quantum critical regime. We have also demonstrated a different kind of logarithmic correction, in the size-dependence of the staggered magnetic susceptibility at the four-dimensional quantum critical point, where we verify the predicted -independent scaling form Kenna2004 ().

The numerical task of finding logarithmic corrections is not a straightforward one. We have established that the appropriate scaling regime is . Within this region, obtaining reliable evidence for logarithmic corrections is critically dependent on having many high-precision data points at very small values of , which mandates accurate calculations at large system sizes. After establishing the location of the critical point to approximately one part in , , we were able to obtain highly accurate extrapolations of the physical observables for coupling ratios as close to as 4.834, i.e. with . This required working with linear system sizes as large as , meaning a system containing interacting spins, and at temperatures as low as . From this perspective, it becomes obvious why previous studies Jin2012 (); tbk07 (); Kao2013 (), with only a handful of data points in the quantum critical regime (none closer than ) were not able to find any meaningful evidence for logarithmic corrections.

Our results are directly relevant to the pressure-induced quantum phase transition in TlCuCl rrs2004 (); Ruegg2008 (); Merchant2014 (). Detailed experiments on this material by elastic and inelastic neutron scattering have measured the staggered magnetization, the Néel temperature, the gap of the quantum disordered phase, and the magnetic excitation spectrum on both sides of the transition. On the assumption that the leading dependence of the control parameter (the ratio of antiferromagnetic superexchange parameters) is linear in the applied pressure, both and show good mean-field exponents and a close linear relation over much of the accessible pressure range Merchant2014 (). On the grounds that the available data follow mean-field scaling around the quantum critical point, it cannot be argued that they provide any evidence for logarithmic corrections, although the size of the experimental errors and the shortage of data very close to the QCP certainly mean they cannot be excluded.

It has been argued very recently Scammell2015 (), based on a field-theoretic treatment, that the apparent suppression of and visible in the experimental data for TlCuCl rather far from the QCP (at pressures 2-4 times the critical pressure) arises due to logarithmic scaling of the coupling constant. In Ref. Merchant2014 () it was assumed that these effects are in fact a consequence of departures from the quantum critical scaling regime, evident also in the violation of linearity between and beyond the point where the order parameter is 60% of the classical moment. Although a direct comparison with our results is not possible without a microscopic treatment of the relationship between the applied pressure and the control parameter, a similar downturn is visible, and better described by including the multiplicative logarithmic corrections, beyond