# Uniform susceptibility of classical antiferromagnets in one and two dimensions in a magnetic field

In recent years, investigations of two-dimensional antiferromagnets concentrated primarily on the quantum model with $`S=1/2`$. A practical reason for that is its possible relevance for the high-temperature superconductivity. On the other hand, the identification with the quantum nonlinear sigma model (QNL$`\sigma `$M) in the low-energy sector allowed using field-theory methods . Although the QNL$`\sigma `$M results for the $`S=1/2`$ model proved to be in a good argeement with quantum Monte Carlo (QMC) simulations (see, e.g., Ref. ), the requirement of low energies confines the validity region of the QNL$`\sigma `$M to rather low temperatures already for $`S1`$. High-temperature series expansions (HTSE) for $`S1`$  and QMC simulations  for $`S=1`$ in the experimentally relevant temperature range, as well as experiments on model substances with $`1S5/2`$, showed much better accord with the pure-quantum self-consistent harmonic approximation (PQSCHA) , than with the field-theoretical QNL$`\sigma `$M predictions. In contrast to the QNL$`\sigma `$M, the PQSCHA maps a quantum system on the corresponding classical system on the lattice, which, in turn, can be studied by classical MC simulations or other methods. The parameters of these classical Hamiltonians are renormalized by quantum fluctuations and given by explicit analytical expressions.

The above arguments show that in most cases the classical model can be used as a good starting point for studying quantum systems. In fact, most of nontrivial features of two-dimensional antiferromagnets, such as impossibility of ordering at nonzero temperatures in the isotropic case, are universal and appear already at the classical level. The main theoretical problem is that due to Goldstone modes, a simple spin-wave theory at $`TJS^2`$ is inapplicable to two-dimensional magnets.

Despite their importance, classical antiferromagnets received much less attention than the quantum $`S=1/2`$ model. In particular, the initial uniform susceptibility $`\chi (T)`$ for the square lattice having a flat maximum at $`TJ`$ has been simulated for $`S=1/2`$ in Refs.  and for $`S=1`$ in Ref. , but there are no results for the classical model yet! For the latter, only the old MC data for the energy  are available up to now.

On the other hand, classical magnets can be theoretically studied with the help of the $`1/D`$ expansion, where $`D`$ is the number of spin components . In Ref. , $`\chi (T)`$ has been calculated for the square lattice and linear chain to first order in $`1/D`$ for all temperatures, the solution interpolating between the exact result at $`T=0`$ and the leading terms of the HTSE at high temperatures. In contrast, the low-energy approaches such as “Schwinger-boson mean-field theory”  or “modified spin-wave theory”  break down at $`TJ`$ and fail to reproduce the maximum of $`\chi (T)`$. It should be noted that for quantum magnets there is a method consisting in the expansion in powers of $`1/N`$ where $`N`$ is the number of flavors in the Schwinger-boson technique . This method, which is nonequivalent to the $`1/D`$ expansion in the limit $`S\mathrm{}`$, is supposed to work for all $`T`$, in contrast to the low-energy QNL$`\sigma `$M. Unfortunately, only the results for $`m(T,H)`$ of ferromagnets  are available.

The $`1/D`$ expansion also works in the situations with nonzero magnetic field, which are not amenable to the methods of Refs.  imposing an external condition $`m=0`$. An especially interesting issue is the singular behavior of $`\chi (H,T)`$ for $`H,T0`$ for the square-lattice and linear-chain models. For any $`H0`$, the spins with lowering temperature come into a position nearly perpendicular to the field, thus $`lim_{H0}lim_{T0}\chi (H,T)=1/(2J_0)`$, where $`J_0`$ is the zero Fourier transform of the exchange interaction, $`J_0=zJ`$, $`z`$ is the number of nearest neighbors. This value coincides with the susceptibility of the three-dimensional classical antiferromagnets on bipartite lattices in the direction transverse to the spontaneous magnetization. For $`H=0`$, the spins assume all directions, including that along the infinitesimal field, for which the susceptibility tends to zero at $`T0`$. Thus $`lim_{T0}lim_{H0}\chi (H,T)=1/(2J_0)(11/D)`$. One can see that the difference between these two results is captured exactly in the first order in $`1/D`$. According to Ref. , for any $`H0`$ with lowering temperature $`\chi (H,T)`$ increases, goes through the flat maximum, decreases, attains a minimum and then goes up to the limiting value $`1/(2J_0)`$.

The existence of the interesting features described above, which should be also pertinent to quantum antiferromagnets, have never been checked numerically. That is why we have undertaken MC simulations for classical AFMs in square-lattice and linear-chain geometries.

Our systems are defined by a classical Heisenberg Hamiltonian

$$=𝐇\underset{i}{}𝐒_i+\frac{1}{2}\underset{ij}{}J_{ij}𝐒_i𝐒_j$$
(1)   
where $`𝐒`$ is a $`D`$-component normalized vector of unit length ($`|𝐒|=1`$), $`𝐇`$ is a magnetic field and the exchange coupling $`J_{ij}`$ is $`J>0`$ for nearest neighbors and zero otherwise. The mean-field transition temperature is given by $`T_c^{\mathrm{MFA}}=J_0/D=zJ/D`$. Although there is no phase transition in our model, it is convenient to choose $`T_c^{\mathrm{MFA}}`$ as the energy scale and to introduce dimensionless temperature, magnetic field, and susceptibilities

$$\theta T/T_c^{\mathrm{MFA}},hH/J_0,\stackrel{~}{\chi }_\alpha J_0\chi _\alpha ,$$
(2)   
where $`\chi _\alpha S_\alpha /H_\alpha `$ and $`\alpha =x,y,z`$.

In the limit $`D\mathrm{}`$, the model Eq. (1) is exactly solvable and equivalent to the spherical model. The solution includes an integral over the Brillouin zone taking into account spin-wave effects in a nonperturbative way. The latter leads to the absence of the phase transition for the spatial dimensionalities $`d2`$.

The $`1/D`$ corrections to the spherical-model solution have been obtained in Refs. . They include double integrals over the Brillouin zone and are responsible for the maximum of the antiferromagnetic susceptibility at $`\theta 1`$ . For small fields and temperatures, $`h,\theta 1`$, the field-induced magnetization $`m`$ for the square-lattice model simplifies to

$$m\frac{h}{2}\left[1\frac{1}{D}+\frac{\theta }{\pi D}\mathrm{ln}\left(1+\frac{h^2}{16}e^{\pi /\theta }\right)+\frac{\theta }{D}\right],$$
(3)   
which follows from Eqs. (4.9) and (2.23) of Ref. . The log term of the above expression is responsible for the singularity of both transverse and longitudinal (with respect to the field) susceptibilities,

$$\stackrel{~}{\chi }_{}m/h,\stackrel{~}{\chi }_{}m/h,$$
(4)   
which was mentioned above. For $`h=0`$ they have the form $`\stackrel{~}{\chi }[11/D+\theta /D]/2`$, whereas for $`h0`$ the limiting value at $`\theta =0`$ and the slope with respect to $`\theta `$ are different: $`\stackrel{~}{\chi }\{1[\theta /(\pi D)]\mathrm{ln}[16/(e^\pi h^2)]\}/2`$. In the latter case, $`\chi `$ has a minimum at $`\theta \theta ^{}=\pi /\mathrm{ln}(16/h^2)`$. There are corrections of order $`\theta ^2`$ and $`1/D^2`$ to Eq. (3). The latter renormalize the last, regular term in Eq. (3) (see Eq. (8.2) of Ref. ). The $`1/D^2`$ corrections cannot, however, appear in the log term of Eq. (3), because this would violate the general properties of $`\chi (H,T)`$ discussed above.

For the linear chain, the magnetization in the region $`h,\theta 1`$ to first order in $`1/D`$ is given by 

$$m\frac{h}{2}\left[1\frac{\theta }{D\sqrt{h^2+\theta ^2}}+\frac{\theta }{D}+O(\theta ^2)\right].$$
(5)   
The transverse susceptibility of the linear chain behaves qualitatively similarly to that of the square lattice. The minimum of $`\chi _{}`$ is attained at $`\theta =h^{2/3}`$ which is smaller than in two dimensions. The longitudinal susceptibility $`\chi _{}`$ corresponding to Eq. (5) has a minimum at $`\theta 3^{1/3}h^{2/3}h`$ and a maximum at $`\theta 3^{1/2}h^{3/2}h`$.

For comparison, the zero-field Takahashi’s results  for the Heisenberg model on the linear chain and square lattice can for $`\theta 1`$ be rewritten in the form 

$$\stackrel{~}{\chi }\frac{1}{3}\{\begin{array}{cc}[1\theta /3]^1,\hfill & d=1\hfill \\ 2\left[1+\sqrt{1(4/3)\theta }\right]^1,\hfill & d=2,\hfill \end{array}$$
(6)   
where the exponentially small terms are neglected. For both lattices the low-temperature expansion is the same to order $`\theta `$: $`\stackrel{~}{\chi }=(1/3)+(1/9)\theta +\mathrm{}`$, and the results diverge at $`\theta 1`$. The coefficient in front of $`\theta `$ here is at variance with the $`1/D`$-expansion results above for $`D=3`$. It was argued in Ref.  that the correct general-$`D`$ form of the low-temperature expansion of the zero-field susceptibility for both square lattice and the linear chain reads

$$\stackrel{~}{\chi }=\frac{1}{2}\left(1\frac{1}{D}\right)+\frac{1}{2D}\left(1\frac{1}{D}\right)\theta +O(\theta ^2),$$
(7)   
i.e., it is reproduced to order $`\theta `$ at the second order of the $`1/D`$ expansion. This formula is in accord with Takahashi’s theory.

In order to check the validity of the analytic results from the $`1/D`$ expansion above for the most realistic case of $`D=3`$, we performed Monte Carlo simulation for three-component classical spins on a chain with length $`N`$ as well as on a square lattice of size $`N=L\times L`$, both with periodic boundary conditions. In our Monte Carlo procedure, a spin is chosen randomly and a trial step is made where the new spin direction is taken randomly with equal distribution on the unit sphere. This trial step does not depend on the initial spin direction. The energy change of the system is computed according Eq. (1) and is accepted with the heat-bath probability. One sweep through the lattice and performing the procedure described above once per spin (on average) is called one Monte Carlo step (MCS). We start our simulation at high temperature and cool the system stepwise. For each temperature we wait 6000MCS (chain) and 4000MCS (square lattice), respectively, in order to reach equilibrium. After thermalization we compute thermal averages $`\mathrm{}`$ for the next 8000MCS (chain) and 6000MCS (square lattice), respectively.

The relevant quantities we are interested in are the magnetization $`mm_z=M_z`$ and the components of the susceptibility $`\chi _\alpha =\frac{N}{T}(M_\alpha ^2M_\alpha ^2)`$, where the $`z`$ axis is directed along H, $`\alpha =x,y,z`$, and $`M_\alpha \frac{1}{N}_iS_i^\alpha `$. We have used the formula above for $`\chi _\alpha `$ to simulate the zero-field and longitudinal susceptibility, $`\chi _{}\chi _z`$. For the transverse susceptibility, $`\chi _{}\chi _x=\chi _y`$, at nonzero field it is more convenient to use Eq. (4). For $`h=0`$ the transverse and longitudinal susceptibilities are identical and calculated as $`\chi _{}=\chi _{}=(\chi _x+\chi _y+\chi _z)/3`$.

With intent to minimize the statistical error and to be able to compute error bars we take averages over $`N_r=100`$ independent Monte Carlo runs. The error bars we show are the mean errors of the averages $`\sigma /\sqrt{N_r}`$, where $`\sigma `$ is the standard deviation of the distribution of thermal averages following from the independent runs.

We start the comparison of theoretical and numerical results with the square lattice. Fig. 1 shows the temperature dependence of the reduced longitudinal susceptibility $`\stackrel{~}{\chi }_{}`$ and reduced transverse susceptibility $`\stackrel{~}{\chi }_{}`$ for different values of the magnetic field, both for the system size $`L=64`$. The corresponding results for the spin chain with system size $`L=100`$ are presented in Fig. 2.

We investigated possible finite-size effects by varying the lattice size. However, we did not find any significant change of our data for lattice sizes in the range $`L=16\mathrm{}64`$ (square lattice) and $`L=40\mathrm{}100`$ (linear chain). Also, we did not find any systematic change of our results for longer Monte Carlo runs so that we believe to present data corresponding to thermal equilibrium.

Note, that for all Monte Carlo data shown the error bars of the transverse susceptibility are smaller than those of the longitudinal one since the transverse susceptibility follows directly from the $`z`$ component of the magnetization while the longitudinal susceptibility is calculated from the fluctuations of the $`z`$ component of the magnetization. In the case $`h=0`$ the transverse and longitudinal susceptibility are identical and follow from fluctuations of the magnetization so that the error bars are larger.

For the square lattice as well as for the chain the numerical data confirm the non-analytic behavior of $`\chi `$ in the limit of temperature $`T0`$, i. e. the limiting values $`\stackrel{~}{\chi }_{}=\stackrel{~}{\chi }_{}=1/2`$ for $`h0`$ and $`\stackrel{~}{\chi }_{}=\stackrel{~}{\chi }_{}=1/3`$ for $`h=0`$.

Especially for the square lattice, the Monte Carlo data agree reasonable with the first-order $`1/D`$ expansion in the whole range of temperatures. On the other hand, at low temperatures the agreement with Takahashi’s theory within error bars is achieved. Our numerical data thus confirm that the coefficient in the linear-$`\theta `$ term in $`\chi `$ in Takahashi’s theory is accurate. For $`h=1`$ and $`\theta 1`$, the MC data fall slightly below the $`1/D`$-expansion curve. Both are again in accord with each other for $`\theta 3`$ (not shown).

The maximum of the longitudinal susceptibility of the square-lattice model for $`h=1`$ looks much sharper than that of the theoretical curve. This feature, as well as the hump on the $`h=0.1`$ curve at slightly lower temperature, are possible indications of the Berezinsky-Kosterlitz-Thouless (BKT) transition. The reason for that is an effective reduction of the number of spin components by one at sufficiently low temperatures in the magnetic field (the effect mentioned in the introduction), so that the Heisenberg model becomes effectively $`D=2`$ and it can undergo a BKT transition in two dimensions. We have not, however, studied this point in detail in this work.

For the antiferromagnetic chain our MC simulation data are in a qualitative agreement with the $`1/D`$ expansion, although the discrepancies are stronger.

Unfortunately, we could also not perform simulations for even lower values of the field $`h`$ for the following reason: The singular behavior of $`\chi `$ stems from the fact that for $`h>0`$ the spins tend to come into a position perpendicular to the field. For fields as small as $`h=0.01`$ (curve 4 in Figures 1 and 2) the amount of energy related to this ordering field is 100 times smaller than the exchange interaction energy. Therefore the corresponding relaxation for this energetically favorable state takes very long in a Monte Carlo simulation, especially for these low temperatures, where this effect occurs for low fields.

Our MC simulations showed for the first time the singular behavior of the susceptibility of classical antiferromagnets at low temperatured and magnetic fields. The results are in accord with predictions based on the first-order $`1/D`$ expansion . It would be interesting to try deriving the corresponding low-temperature results \[cf. Eqs. (3) and (5)\] without using the $`1/D`$ expansion. One of the formulas of this type already exists: It is Eq. (7). A candidate among theoretical approaches is the chiral perturbation theory of Ref. , which is applicable to quantum models, as well.

The features manifested here by classical antiferromagnets should be pertinent to quantum models, as well. The effects observed here could be checked with the help of the QMC simulations which achieved recently a substantial accuracy (see, e.g., Refs.  and ). Another possibility is to map the quantum model on the classical one  and to perform classical MC simulations. One should also mention an alternative way of mapping of quantum magnetic Hamiltonians on classical ones with the help of the coherent-state cumulant expansion , which is a rigorous expansion in powers of $`1/S`$.