Finite size scaling for a first order transition where a continuous symmetry is broken: The spin-flop transition in the 3D XXZ Heisenberg antiferromagnet

TL;DR

This paper develops a finite size scaling theory for first order phase transitions involving continuous symmetry breaking, validated by Monte Carlo simulations of the 3D XXZ Heisenberg antiferromagnet's spin-flop transition.

Contribution

It introduces a Gaussian approximation with a degeneracy factor to predict universal finite size scaling behavior at first order transitions with symmetry breaking.

Findings

Universal intersection points in moments and cumulants at the transition.

Predicted degeneracy factor q=π matches numerical evidence.

Finite size corrections scale as inverse volume.

Abstract

Finite size scaling for a first order phase transition where a continuous symmetry is broken is developed using an approximation of Gaussian probability distributions with a phenomenological "degeneracy" factor included. Predictions are compared with data from Monte Carlo simulations of the three-dimensional, XXZ Heisenberg antiferromagnet in a field in order to study the finite size behavior on a $L \times L \times L$ simple cubic lattice for the first order "spin-flop" transition between the Ising-like antiferromagnetic state and the canted, XY-like state. Our theory predicts that for large linear dimension $L$ the field dependence of all moments of the order parameters as well as the fourth-order cumulants exhibit universal intersections. Corrections to leading order should scale as the inverse volume. The values of these intersections at the spin-flop transition point can be expressed in terms of a factor $q$ that characterizes the relative degeneracy of the ordered phases. Our theory yields $q=\pi$, and we present numerical evidence that is compatible with this prediction. The agreement between the theory and simulation implies a heretofore unknown universality can be invoked for first order phase transitions.

Tables

Table 1. Table 1: Estimates for q e ​ f ​ f subscript 𝑞 𝑒 𝑓 𝑓 q_{eff} from the ratio of probability distributions of the weights of the peaks for different values of L 𝐿 L .

$L$	$H^t/J=3.83830$	$H^t/J=3.838305$	$H^t/J=3.83831$
60	$3.31\left(23\right)$	$3.42\left(30\right)$	$3.53\left(28\right)$
80	$3.15\left(24\right)$	$3.36\left(26\right)$	$3.75\left(29\right)$
100	$2.82\left(30\right)$	$3.19\left(29\right)$	$3.82\left(30\right)$

Equations

$$\mathcal{H}=J\sum _{⟨i,j⟩}[\Delta (S_{ix}{S}_{jx}+S_{iy}{S}_{jy})+S_{iz}{S}_{jz}]-H\sum _i{S}_{iz}$$

$$\tilde{m}_{z}=\frac{1}{L^{3}}\Big{(}\sum\limits_{i\in 1}S_{iz}-\sum\limits_{i\in 2}S_{iz}\Big{)}\enskip.$$

$$m_{z}=\frac{1}{L^{3}}\Big{(}\sum\limits_{i\in 1}S_{iz}+\sum\limits_{i\in 2}S_{iz}\Big{)}\enskip.$$

$$\psi_{\alpha}=\frac{1}{L^{3}}\Big{(}\sum\limits_{i\in 1}S_{i\alpha}-\sum\limits_{i\in 2}S_{i\alpha}\Big{)}\enskip,\quad\alpha=(x,y)\enskip.$$

$$m_{z,\infty }=-(\partial F/\partial H{)}_T.$$

$$\Delta F_{AF}=(H^t-H)m_{z,\infty }^{AF},H\le H^t,$$

$$\Delta F_{SF}=(H^t-H)m_{z,\infty }^{SF},H\ge H^t.$$

$$\Delta F\equiv \Delta F_{SF}-\Delta F_{AF}=(H-H^t)(m_{z,\infty }^{AF}-m_{z,\infty }^{SF})$$

$$=-(H-H^t)\Delta m$$

$$a_{AF}=\exp (\Delta F{L}^3/k_{B}T)/[q+\exp (\Delta F{L}^3/k_{B}T)],$$

$$a_{SF}=q/[q+\exp (\Delta F{L}^3/k_{B}T)],$$

$$P(\vec{\psi })=\delta (\vec{\psi })+\sum _{k=1}^q\delta (\vec{\psi }-{\vec{\psi }}_k)$$

$$P({\tilde{m}}_z,\psi )=[\delta ({\tilde{m}}_z-{\tilde{m}}_{z,\infty })+\delta ({\tilde{m}}_z+{\tilde{m}}_{z,\infty })]\delta (\psi )+{\int }_0^{2\pi }\delta ({\tilde{m}}_z)\delta (\psi -{\psi }_{\infty })d\varphi$$

$$P(\psi )=2\delta (\psi )+2\pi \delta (\psi -{\psi }_{\infty }),$$

$$P_{\infty }({\tilde{m}}_z,\psi )=[\delta ({\tilde{m}}_z-{\tilde{m}}_{z,\infty })+\delta ({\tilde{m}}_z+{\tilde{m}}_{z,\infty })]\delta (\psi )+2q\delta ({\tilde{m}}_z)\delta (\psi -{\psi }_{\infty }).$$

$${\mathcal{N}}_{\infty }=\underset{-1}{\overset{+1}{\int }}d{\tilde{m}}_z\underset{0^{-}}{\overset{1}{\int }}d\psi P_{\infty }({\tilde{m}}_z,\psi )=2+2q.$$

$$⟨|\vec{\psi }|{⟩}_{\infty }={\psi }_{\infty }q/(1+q),$$

$$⟨{\psi }^2{⟩}_{\infty }={\psi }_{\infty }^{2}q/(1+q),$$

$$⟨{\psi }^4{⟩}_{\infty }={\psi }_{\infty }^{4}q/(1+q),$$

$$⟨|{\tilde{m}}_z|{⟩}_{\infty }={\tilde{m}}_{z,\infty }/(1+q),$$

$$⟨{\tilde{m}}_z^2{⟩}_{\infty }={\tilde{m}}_{z,\infty }^2/(1+q),$$

$$⟨{\tilde{m}}_z^4{⟩}_{\infty }={\tilde{m}}_{z,\infty }^4/(1+q).$$

$$U_{\infty }^{xy}\equiv 1-⟨{\psi }^4{⟩}_{\infty }/[3⟨{\psi }^2{⟩}_{\infty }^2]=1-\frac{1+q}{3q},$$

$$U_{\infty }^z\equiv 1-⟨{\tilde{m}}_z^4{⟩}_{\infty }/[3⟨{\tilde{m}}_z^2{⟩}_{\infty }^2]=(2-q)/3.$$

$$\displaystyle P_{L}(m_{z})\propto\frac{a_{AF}}{\sqrt{\chi^{AF}_{zz}}}\exp\Big{\{}-\frac{[m_{z}-(m^{AF}_{z,{\infty}}+\chi^{AF}_{zz}\Delta H)]^{2}}{2k_{B}T\chi^{AF}_{zz}/L^{3}}\Big{\}}$$

$$\displaystyle\quad+\frac{a_{SF}}{\sqrt{\chi^{SF}_{zz}}}\exp\Big{\{}-\frac{[m_{z}-(m^{SF}_{z,{\infty}}+\chi^{SF}_{zz}\Delta H)]^{2}}{2k_{B}T\chi^{SF}_{zz}/L^{3}}\Big{\}}\enskip,$$

$$\Delta H\equiv H-H^t$$

$${\chi }_{zz}^{\max }\approx \frac{{\chi }_{zz}^{AF}+{\chi }_{zz}^{SF}}{2}+\frac{(\Delta m{)}^2{L}^3}{4k_{B}T}$$

$$(H^{\max }-H^t)/k_{B}T=-[\Delta m{L}^3{]}^{-1}\ln q.$$

$$P^{AF}_{L}(\vec{\psi})=\mathcal{N}\exp\Big{(}-\frac{{\vec{\psi}}{{\>}^{2}}}{2k_{B}T\tilde{\chi}^{AF}_{xy}/L^{3}}\Big{)}\enskip,$$

$$P^{SF}_{L}(\vec{\psi})=\mathcal{N}\exp\Big{[}-\frac{({\vec{\psi}}{{\>}^{2}}-\psi^{2}_{\infty})^{2}}{8{\psi^{2}_{\infty}}k_{B}T\tilde{\chi}^{SF}_{xy}/L^{3}}\Big{]}\enskip,$$

$$⟨{\psi }^2{⟩}_L=a_{AF}⟨{\psi }^2{⟩}_{AF}+(1-a_{AF})⟨{\psi }^2{⟩}_{SF}$$

$$⟨{\psi }^4{⟩}_L=a_{AF}⟨{\psi }^4{⟩}_{AF}+(1-a_{AF})⟨{\psi }^4{⟩}_{SF}$$

$$a_{AF}=1/[1+q\exp (\mathcal{Z})].$$

Full text

I Introduction

Finite size scaling at both first order and second order transitions between phases with discrete numbers of states is now relatively well established and extremely successful at describing phase transition behavior in the thermodynamic limit from Monte Carlo data produced for finite size systems  MEF1 ; VP ; DPL ; BL1 ; BK ; CLB ; MCbook . The simplest case of a first order transition, namely the first order, field driven transition in the 2-dimensional Ising model below its critical point, was studied using Monte Carlo simulations by Binder and Landau BL1 ; and to a good approximation the probability distribution of the order parameter at the first order transition could be described by the sum of the Gaussians representing the two coexisting states in the finite system. For the temperature driven first-order transition in the $q$-state Potts model a similar theoretical development could be used, but a factor of “$q$” needed to be included in the Gaussian representing the $q$-fold degenerate ordered state BK . Monte Carlo simulations of the $q=10$ Potts model on $L\times L$ lattices verified this finite size behavior CLB .

In the case of a first order phase transition involving the breaking of a continuous symmetry, however, there are neither good data from simulations nor theoretical predictions regarding the finite size behavior. A good “testing ground” system is thus needed to help provide an understanding of this case, and we believe that the uniaxially anisotropic, three-dimensional ($3D$) Heisenberg antiferromagnet in an external field, $H$, is exactly such a candidate model. It has attracted substantial interest for a number of decades, largely to clarify the phase diagram and ordered structures for this model and to identify the nature of the multicritical point MEF ; DRJMME ; KNF ; DLKB ; MHJ ; PCAPEV ; RFYH ; GBWS ; WS ; HTL . Although different scenarios have been proposed, we now believe that the phase diagram of this model contains a low temperature, low field antiferromagnetic (AF) phase in which the spins point in opposite directions along the axis given by the anisotropy as shown in Fig. 1(a) and a spin-flop (SF) phase in which the spins are tilted with continuous rotational symmetry about the field direction (see Fig. 1(b)) at low $T$ and higher $H$. A paramagnetic (P) phase with no long range order exists at high $T$ and/or at high $H$.

The SF to P and AF to P phase transition lines are of second order and belong to the XY and the Ising universality classes, respectively. In contrast, a line of first-order transitions separates the AF and SF phases. The point $T=T_{b}$ where the three phases meet was determined to be a bicritical point in the three-dimensional (3D) Heisenberg universality class. The resultant phase diagram in the vicinity of the bicritical point is shown in Fig. 2. In earlier work the spin-flop boundary for the anisotropic Heisenberg antiferromagnet in a field was located rather precisely  WS ; HTL , so this model is indeed a fertile testing ground for the study of finite size effects at a first order transition where a continuous symmetry is broken.

In the next section we will clearly define the model to be studied, describe the Monte Carlo methodology used to generate data for different lattice sizes, and provide a theoretical formulation for the understanding of the finite size behavior of this model (mathematical details will be assigned to Appendix A). In Sec.III we will present the results of careful Monte Carlo simulations which are then used to test the theory. Our conclusions are summarized in Sec.IV.

II Model and Methods

II.1 Model

The XXZ antiferromagnetic model studied here is described by the Hamiltonian in Eq. (1).

[TABLE]

where the classical spins Si are unit vectors with components $(S_{ix},S_{iy},S_{iz})$ on sites $i$ of a simple cubic lattice with linear size $L$, and $J>0$ is the exchange coupling between nearest-neighbor pairs of spins. The first summation is over all $\langle i,j\rangle$ pairs of nearest-neighbor sites and the second summation is over all $N=L^{3}$ spins on the lattice. $\Delta$ is the uniaxial exchange anisotropy, which we set to $\Delta=0.8$ in this work. An external magnetic field $H$ is applied along the $z-$axis, which is the easy axis of the model.

The phase diagram in the temperature $T$ and field $H$ plane is shown in Fig. 2, where the meeting point of the three phase transition lines is estimated to be WS ; HTL at $k_{B}T_{b}/J=1.025\pm 0.0025$ and $H_{b}/J=3.89\pm 0.01$. Note that here we are not concerned with the phase boundaries near the bicritical point but rather shall study the finite size effects associated with the first order transition from the AF phase to the SF phase at $T=0.95J/k_{B}$ at a transition field $H^{t}$.

II.2 Monte Carlo Methods

To carry out Monte Carlo simulations we used two different Monte Carlo methods and considered simple cubic lattices with even values of $L$ and periodic boundary conditions. Simulations for $L\leq 60$ were first performed using a standard Metropolis algorithm metropolis with the R1279 shift register random number generator MCbook . Runs of length $3\times 10^{7}$ MCS were performed for all lattice sizes and the number of independent runs ranges from $10$ for $L=30$ to $1035$ for $L=60$. For $L=60$ Metropolis sampling had great difficulty tunneling between the two states on opposite sides of the spin-flop transitions. Therefore, to insure that the sampling was truly ergodic for $L=60$, $L=80$ and $L=100$, we implemented multicanonical sampling BergNeuhaus . Multicanonical simulations were performed for the entire range of sizes so that results could be compared with those from Metropolis sampling. The multicanonical sampling probability was determined iteratively for each $L$ and then runs of length $10^{7}$ MCS were carried out. To determine averages and error bars a total of $100$ independent runs were made for $L=30$ and the number increased with increasing size until $900$ independent runs were used for $L=100$. For smaller lattices, there was agreement between the data generated using the two different sampling methods and the results could be combined for the analysis. For the multicanonical runs the Mersenne Twister random number generator was used MCbook .

We employed histogram reweighting techniques reweighting to extract thermodynamic quantities for fields near the values used in the simulations. For the largest lattices, comparisons were made between runs made at adjacent field values and reweighted results to insure that we were not reweighting beyond the reliable range of fields.

As we shall see shortly only data for $L\geq 40$ were in the asymptotic finite size scaling regime, so smaller lattice data were omitted from plots of raw data that will be shown in the results section of this manuscript.

II.3 Theory

II.3.1 General relations

For small enough fields, $H<H^{t}(T)$, and low temperatures, $T$, the anisotropic Heisenberg antiferromagnet exhibits Neél-type two-sublattice order on the simple cubic (or other bipartite three-dimensional) lattices. This order is described by the staggered magnetization (with the two interpenetrating sublattices of the $L\times L\times L$ lattice denoted by indices $1$ and $2$ )

[TABLE]

For $H>0$ we also expect to have a non-zero uniform magnetization $m_{z}$ (per spin)

[TABLE]

At the transition field, $H=H^{t}(T)$, there is a first-order phase transition to the “spin-flop” phase described by a two-component order involving the transverse spin components

[TABLE]

Both $\psi_{x}$, $\psi_{y}$ are equivalent and form the components of a vector order parameter $\vec{\psi}$ with XY symmetry. Fig. 2 shows the phase diagram, and Fig. 3 describes the schematic variation of the free energy with the magnetic field $H$ (for simplicity, the temperature dependence of $H^{t}(T)$ is suppressed in Fig. 3) along with its derivative $m_{z,\infty}=\langle m_{z}\rangle_{T,L\rightarrow\infty}$ as well as the order parameter $\psi_{\infty}=\sqrt{\langle\psi^{2}_{x}+\psi^{2}_{y}\rangle_{T,L\rightarrow\infty}}$. Henceforth, a subscript “$\infty$” means that we are referring to properties in the thermodynamic limit. Note that the variable $m_{z,\infty}$ is the thermodynamically conjugate variable to the magnetic field, with $F$ being the Gibbs free energy per spin free_energy ,

[TABLE]

Defining the transition field in the thermodynamic limit as $H^{t}$, we note that $m_{z,\infty}$ at $H^{t}$ must jump from $m^{AF}_{z,\infty}$ in the antiferromagnetic phase to a larger value $m^{SF}_{z,\infty}$ in the spin-flop phase. At the transition field the free energies of the two phases are equal, $F=F_{t}$, and the free energy differences (per spin) in the two phases relative to this value are, to leading order

[TABLE]

[TABLE]

The total free energy difference between phases then becomes (continuing both phases into their metastable region and ignoring the inequalities in Eqs. (6), (7))

[TABLE]

where we have introduced $\Delta{m}=(m^{SF}_{z,\infty}-m^{AF}_{z,\infty})$ to represent the jump in the magnetization $m_{z}$ in the thermodynamic limit at the transition (see Fig. 3). Note that $\Delta F>0$ for $H<H^{t}$ since $m^{AF}_{z,\infty}<m^{SF}_{z,\infty}$.

We can use Eqs. (6) - (II.3.1) to construct the statistical weights of the two phases in a large, but finite, system in a field $H\approx H^{t}$. The most naive assumption would be $a_{AF}\propto\exp(-\Delta F_{AF}L^{3}/k_{B}T)$, $a_{SF}\propto\exp(-\Delta F_{SF}L^{3}/k_{B}T)$. However, this assumption, using a common normalization factor, disregards the difference in degeneracy of the two phases. While the $AF$ phase is two-fold degenerate ($n=1$, one-component order parameter), in the SF phase a continuous (XY-model like) symmetry is broken ($n=2$, two-component order parameter). How this difference enters in the weights is not obvious, unlike in the simpler case of the thermally driven $q$-state Potts model where the high-temperature phase is non-degenerate and the low-temperature phase is simply $q$-fold degenerate. There, an extra factor $q$ appears in the weight of the low-temperature phase multiplying the Boltzmann factor.

It is unclear (at least to us) what this factor $q$ must be when dealing with a continuous symmetry. Thus, we introduce an analogous factor $q$ $\it phenomenologically$ to account for the difference in degeneracy between the phases at the coexistence point in the thermodynamic limit $(H=H^{t})$. As we shall see, in full analogy with the Potts model, this factor leads to shifts of characteristic finite size induced features (e.g. position of the maximum of the magnetic susceptibility or specific heat, minimum of the various cumulants, etc.). Observing these features in the simulations should provide numerical estimates for this degeneracy factor $q$ for the present problem. In the next subsection we shall advance a hypothesis for the value of $q$.

We now postulate the “equal weight rule” BK for the statistical weights of the two phases, i.e. $a_{AF}=\mathcal{N}\exp(-\Delta F_{AF}L^{3}/k_{B}T)$ and $a_{SF}=\mathcal{N}q\exp(-\Delta F_{SF}L^{3}/k_{B}T)$, with a normalization factor $\mathcal{N}$. Requiring $a_{AF}+a_{SF}=1$ yields

[TABLE]

[TABLE]

as expected.

Formally, the factor $q$ can be absorbed by redefining the weights as $a_{AF}=\exp(\Delta F^{\prime}L^{3}/k_{B}T)/[1+\exp(\Delta F^{\prime}L^{3}/k_{B}T)]$, $a_{SF}=1/[1+\exp(\Delta F^{\prime}L^{3}/k_{B}T)]$, with $\Delta F^{\prime}=\Delta F-(k_{B}T/L^{3})\ln q$. This shows that finite size induced shifts of characteristic features scaling as ($k_{B}T/L^{3})\ln q$ will occur.

II.3.2 Order parameter distribution at the transition in the thermodynamic limit

Before going further, we use the transition between the disordered and ordered states for the $q$-state Potts model to provide some insight about the effective degeneracy factor for the AF to SF transition in the anisotropic Heisenberg model. For the $q$-state Potts model in the thermodynamic limit exactly at the transition temperature, the probability distribution of the order parameter $P(\vec{\psi})$ is simply the sum of $q+1$ weighted delta functions,

[TABLE]

where the $\vec{\psi}_{k}$ are the discrete values of the order parameter in the ordered phase Wu . The first term on the right hand side of Eq. (11) represents the disordered phase and the second term represents the (degenerate) ordered phase. In the current case a similar expression holds except that the $\vec{\psi}$ are continuous. We, therefore, conjecture that

[TABLE]

where the order parameter $\vec{\psi}$ is written in terms of the magnitude $\psi$ and angle $\phi$ in the $(\psi_{x},\psi_{y})$ plane and we have integrated over $\phi$. The index “$\infty$” indicates that the thermodynamic limit was taken first and then $H\to H^{t}$. Since there is no dependence upon $\phi$, the integral gives $2\pi$. Integration over $\tilde{m}_{z}$ then yields

[TABLE]

or the relative weight of the two phases is simply $\pi$!

Eq. (12) merely indicates that in the thermodynamic limit and for $H=H^{t}$, we have phase coexistence between pure AF phases ($\tilde{m}_{z}=\pm\tilde{m}_{z,\infty}$, $\vec{\psi}=0$) and pure SF phases ($\vec{\psi}=(\psi_{\infty},\phi)$, in polar coordinates in the ($\psi_{x},\psi_{y})$-plane, and $\tilde{m}_{z}=0$). The distribution of the order parameters is simply characterized by the appropriate Dirac delta-functions. Making contact with formulation of Eqs. (9), (10), where the relative weights of the two phases at $H^{t}$ was denoted by the phenomenological parameter $q$, we find that the joint (unnormalized) distribution of the order parameters $\tilde{m}_{z}$, $\psi=|\vec{\psi}|$ becomes

[TABLE]

The normalization constant for this distribution is

[TABLE]

From Eqs. (14), (15) we can easily obtain the moments and cumulants of both order parameters (the notation $\langle\cdots\rangle_{\infty}$ means that an average over both phases at the transition point in the thermodynamic limit is taken)

[TABLE]

[TABLE]

Hence, the cumulants simply become

[TABLE]

and

[TABLE]

Of course, these results do not invoke the assumption of Gaussian distributions of the order parameters for finite $L$ {see Eqs. (II.3.3), (28), (29) or a similar assumption for $P_{L}(\tilde{m}_{z})$ that will be used below}. Thus, Eqs. (16)-(23) are not affected in any way by deviations from Gaussian distributions in the wings of the actual distributions for finite $L$.

Eqs. (16)-(23) permit stringent tests of this theory by simulations using the following recipe: Suppose an accurate estimate of $H^{t}$ is known from suitable finite size extrapolation (e.g. using Eq. (27); other choices will be given below). Then, a very large system can be simulated (for which no transitions between the pure phases occur for very long runs) right at $H=H^{t}$, once starting in the AF phase and once starting in the SF phase, to obtain very accurate estimates of $\tilde{m}_{z,\infty}$ and of $\psi_{\infty}$. From the distributions of the order parameters in these pure phases, accurate estimates of the staggered susceptibilities $\tilde{\chi}^{AF}_{xy}$, $\tilde{\chi}^{SF}_{xy}$, $\tilde{\chi}^{AF}_{zz}$, and $\tilde{\chi}^{SF}_{zz}$, can also be extracted. Extrapolations of the estimates for $\langle|\vec{\psi}|\rangle_{L}$, $\langle|\tilde{m_{z}}|\rangle_{L}$, $\langle\psi^{2}\rangle_{L}$, $\langle\tilde{m}^{2}_{z}\rangle_{L}$, $U^{xy}_{L}$ and $U^{z}_{L}$ at $H=H^{t}$ towards $L=\infty$ should provide estimates for the factors $q/(1+q)$ and $1/(1+q)$ in Eqs. (16), (19) as well as the cumulants, Eqs. (22), (23).

Eqs. (16) - (23) also define (almost) universal intersection points when we study analogous averages for finite $L$ as a function of the field $H$. $\langle|\vec{\psi}|\rangle_{L}$, $\langle|\tilde{m}_{z}|\rangle_{L}$, $\langle\psi^{2}\rangle_{L}$, $\langle\tilde{m}^{2}_{z}\rangle_{L}$, $U^{xy}_{L}$ and $U^{z}_{L}$ are all analytic functions of $H$, saturating for $|\Delta F|L^{3}/k_{B}T>>1$. For the cumulants these saturation values are trivial, e.g. $U^{z}(H<H^{t})=1$, $U^{z}(H>H^{t})=0$. But we will show later that they will agree with Eqs. (16) - (23), up to corrections of order $L^{-\frac{3}{2}}$ or $L^{-3}$ for $|\Delta F|L^{3}/k_{B}T=0$. Unlike second-order transitions for which only cumulants, e.g. $U_{L}^{xy}$ and $U_{L}^{z}$, have unique intersection points at the transition, for this first order transition both the individual moments and the cumulants exhibit this feature of common intersection points at $H=H^{t}$.

II.3.3 Two-Gaussian approximation for the magnetization distribution

While Eqs. (16)-(23) describe the behavior of the system when we first set $H=H^{t}$ and then take the limit $L\rightarrow\infty$, it is also of great interest to explore the leading corrections to the limiting behavior when $L$ is large but finite.

Following general considerations of statistical physics L&L , in pure phases, for large but finite size systems, we expect Gaussian distributions for the densities of extensive thermodynamic variables rather than $\delta$-functions. A simple case to consider is the uniform magnetization for which Gaussian distributions for the (scalar) quantity $m_{z}$ in the two phases would give a distribution,

[TABLE]

where

[TABLE]

and $\chi^{AF}_{zz},$ $\chi^{SF}_{zz}$ are the susceptibilities at $(H=H^{t})$ in the two phases.

Invoking the analogy of these equations to the case of the Potts model energy distribution [cf. Eq. (IV.21) in Ref. MCbook ], we conclude that the susceptibility peak should scale as

[TABLE]

in analogy to the specific heat of the Potts model. Note that the location of this maximal magnetization fluctuation (i.e. susceptibility of the z-component of the uniform magnetization) occurs when the two weights are equal, i.e. $a_{AF}=a_{SF}=1/2$. This condition readily yields $q\exp(-\Delta FL^{3}/k_{B}T)=1$, i.e. $\Delta F/k_{B}T=\ln q/L^{3}$, or

[TABLE]

Since $m^{AF}_{z,\infty}<m^{SF}_{z,\infty}$, the position of the susceptibility maximum relative to the transition point must shift to smaller fields, $H^{\rm max}<H^{t}$, and scale with size like $L^{-3}$.

The susceptibility at the transition point $H^{t}$ is smaller by a factor $4q/(1+q)^{2}$ than $\chi_{zz}^{\max}$ for $L\to\infty$.

One important caveat, however, is that for the weight $a_{-}$ for the low temperature phase of the Potts model, the factor $q$ reflecting the degeneracy of the ordered phase is known. In contrast, here the value of a similar factor representing the difference in degeneracies of the spin-flop and antiferromagnetic phases is $\it unknown$ unless we rely on the hypothesis of Eq. (13) that $q=\pi$.

II.3.4 The SF phase order parameter distribution

We next consider the distribution of the SF order parameter $\vec{\psi}$. For $H<H^{t}$, i.e. in the AF phase, there is simply a Gaussian distribution about zero since the transverse spin component is disordered,

[TABLE]

where we introduced the notation $\tilde{\chi}_{\alpha\beta}$ for the tensor of staggered susceptibilities, and $\tilde{\chi}^{AF}_{xy}$ stands for the $xy$-components of the staggered susceptibility in the AF phase.

The order parameter distribution in the SF phase is more interesting,

[TABLE]

where now $\tilde{\chi}^{SF}_{xy}$ denotes the $xy$-component of the staggered susceptibility in the spin-flop phase. Note that a 4th order polynomial in $\vec{\psi}$ is needed in the argument of the exponential function in Eq. (29) to bring out the spherical symmetry in the $(\psi_{x},\psi_{y})$-plane correctly. Near the peak $(|\vec{\psi}|\approx\psi_{\infty})$ the argument of the exponential reduces to the expected quadratic form, i.e. $-(\psi-\psi_{\infty})^{2}/(2k_{B}T\tilde{\chi}_{xy}^{SF}/L^{3})$.

In the vicinity of the transition field, $H^{t}$, we now make the standard superposition approximation,

[TABLE]

where

[TABLE]

We have written $\mathcal{Z}$ =$-\Delta{F}L^{3}/k_{B}T$, and the moments $\langle\cdots\rangle_{AF},$ $\langle\cdots\rangle_{SF}$ refer to the order parameter distributions in the “pure” AF and SF phases respectively.

Of course, we could repeat the calculation of Eqs. (14)-(23) for finite $L$, replacing the delta functions by the appropriate Gaussian distributions, e.g. Eqs. (28) and  (29). These calculations are straightforward, but clumsy, so we confine the details to Appendix A and only give a few final results here.

From the general result for the fourth order cumulant (see Appendix A) we immediately conclude that for $H^{t}$, $U^{xy}_{L}$ differs from $U^{*}$=$U^{xy}_{\infty}$ only by corrections of order $L^{-3}$ and

[TABLE]

If we use the value $q=\pi$, we find that $U_{\infty}^{xy}\approx 0.56056$.

Taking the derivative of the expression for the cumulant to find the minimum, and writing $Y$=$q\exp(\mathcal{Z})$, we find that the location of the minimum is given by $Y_{\min}\approx 2k_{B}T\tilde{\chi}^{AF}_{xy}/((\psi_{\infty})^{2}L^{3})$. Hence

[TABLE]

Thus, the shift in the location of the minimum scales as $L^{-3}$, but the leading term is actually ${\ln}L/L^{3}$. At the minimum, the value of $U_{L}^{xy}$ is

[TABLE]

which means that $U_{L,\min}^{xy}$ approaches $-\infty$ proportional to $-L^{3}$ as $L$ approaches infinity.

II.3.5 AF phase order parameter distribution

The AF order parameter distribution for this one-component order parameter in the AF phase is a double Gaussian, analogous to Eq. (II.3.3),

[TABLE]

whereas in the SF phase it is given by a single Gaussian

[TABLE]

From these distribution and superposition approximations analogous to Eqs. (30) - (32), it is straightforward to calculate the various moments of the staggered magnetization and the cumulant (see Appendix A). Then, at the transition, the fixed point value of the cumulant is

[TABLE]

This means that, as expected, the asymptotic value of the cumulant coincides with the result from the treatment using delta function distributions. Again the correction to this result is small, of order $L^{-3}$; and choosing $q=\pi$, we find the fixed point value of the cumulant is $U_{*}^{z}\approx-0.38$.

II.3.6 Maximum slope of the cumulant and cumulant crossings

We have shown that the cumulant at $H=H^{t}$ is of order unity but reaches a deep minimum (of order $-L^{3}$ for $L\rightarrow\infty$) at a value $H^{t}-H\propto L^{-3}$ (on the AF side of the spin-flop transition when the cumulant of the SF-order parameter is considered). In order to achieve a variation of order $L^{+3}$ in an interval of order $L^{-3}$, the maximum slope of the cumulant in the interval between the minimum and the crossing point must then be of order $L^{6}$. On the other hand, we can easily show that the slope of the cumulant right at $H=H^{t}$ still is only of order $L^{3}$, namely

[TABLE]

Thus, significant curvature should appear in the plot of $U_{L}^{xy}$ vs. $H$ near $H^{t}$, since the slope first increases from zero at the cumulant minimum to a value of order $L^{6}$ and then decreases to a value of order $L^{3}$ at $H=H^{t}$. The location of the intersection point, which corresponds to a cumulant value of order unity, hence must be very close to $H^{t}$ ($|H^{\rm cross}-H^{t}|\propto L^{-6}$) also. On a scale of $H-H^{t}$ of order $L^{-3}$ (the regime over which the variable $\mathcal{Z}$ exhibits a significant variation) differences in cumulant intersections of order $L^{-6}$ are completely negligible. Consequently, the spread of the cumulant intersections is very small, although in a strict sense there is no unique cumulant crossing point. Two cumulants for linear dimensions $L$ and $L^{\prime}=L+\delta$ that differ slightly can be shown to cross for $H^{\rm cross}$ given by

[TABLE]

This equation is derived by a Taylor expansion of $U_{L+\delta}$ simultaneously in the small variables $H-H^{t}$ and $\delta/L$. This result verifies the above argument that the scale for the shift of the cumulant crossing is negligibly small in comparison with that for the shift of the cumulant minimum. The shift of the intersection point thus scales as the square of the inverse volume, while the regime over which the transition is spread out is given by

[TABLE]

The situation is analogous to the case of the (temperature driven) transition in the Potts model, cf. Vollmayr et al.  Vollmayr . As in the latter case, and unlike the simple, field driven first order transition of an Ising ferromagnet below $T_{c}$, there is no “equal height rule” of the order parameter distribution at the transition. Eqs. (28), (A) show that the peak height at $\vec{\psi}=0$ scales proportional to the volume, $L^{3}$, while the height of the “rim” $|\vec{\psi}|=\psi_{\infty}$ only scales like the square root of the volume, $L^{3/2}$, cf. Eqs. (13), (14), and (23) of Vollmayr at al. Vollmayr

Useful information can be readily extracted about the mean square order parameter $\langle\psi^{2}\rangle$ as a function of the field $H$ near the transition field $H^{t}$.

The superposition approximation (Eqs. (30), (31)) readily shows that $\langle\psi^{2}\rangle$ is described by a simple scaling function of $\mathcal{Z}$, namely

[TABLE]

where we have introduced the abbreviation

[TABLE]

Note that the quantities $\psi_{\infty}$, $\tilde{\chi}^{AF}_{xy}$ (and $\Delta{m}$ which is needed to convert the scale of $H-H^{t}$ to $\mathcal{Z}$) can be estimated directly from simulations. For runs of modest length for very large systems precisely at $H=H^{t}$, transitions between the phases can be avoided due to metastability. Starting in the AF and SF states respectively will permit measurements in the pure phases. ($H^{t}$ is already known with high precision). For the correct choice of $q$, all choices of (sufficiently large) $L$ should then lead to perfect collapse on a master curve, Eq. (42), that is explicitly predicted.

The value of the order parameter at the location of the maximum slope is

[TABLE]

which is independent of the value of “$q$”. Since $\psi_{\infty}$ and $\tilde{\chi}_{xy}^{AF}$ can be measured independently, we have a non-trivial test of the double Gaussian approximation. The value of the slope at the transition field $H=H^{t}$ is

[TABLE]

We also note that for $H=H^{t}$, i.e. $\mathcal{Z}=0$, Eq. (42) predicts $\langle\psi^{2}\rangle/{\psi_{\infty}^{2}}=(b+q)/(1+q)\approx{q/(1+q)}$, as expected from Eq. (17). When $\langle\psi^{2}\rangle/{\psi_{\infty}^{2}}$ is plotted vs $H$, all curves for large $L$ will intersect for $H=H^{t}$ at this value.

II.3.7 Limitations of the phenomenological theory

At this point, we comment on an important distinction between the order parameter cumulant intersection for first order transitions and for second order transitions. At second order transitions, for $L\rightarrow\infty$ corrections to scaling can be ignored and the cumulant is a regular function of the variable $\mathcal{Z}=tL^{1/\nu}$, $\nu$ being the correlation length critical exponent. The slope of the cumulant $dU_{L}/dt$, $t$ being the reduced distance from the transition point, at the transition point hence is of the same order as the inverse shift $L^{-1/\nu}$ of the susceptibility maximum.

The $L^{-3}$ correction to Eqs. (33), (38) is simply a “correction to finite size scaling” for first order transitions, analogous to those that appear in finite size scaling at second order transitions. However, in the latter case these corrections involve a second, non-trivial exponent (different from $1/\nu$), while here the inverse volume $L^{-3}$ is the only variable that leads to correction terms in the finite size scaling description. Since the resolution of the numerical data is insufficient to quantitatively resolve any of these correction term effects, they shall not be discussed further. We also note that further terms are expected due to corrections to the Gaussian approximation for the probability distribution, e.g. the result $U^{AF}_{L}=1/3$ for the cumulant of the XY-order parameter in the AF phase is also expected to have a $1/L^{3}$ correction (related to higher order correlation functions) which has been ignored here but seems rather relevant numerically. In addition, corrections of order $\exp(-L/\xi)$, where $\xi$ is the appropriate correlation length of the AF and/or SF order, have been ignored. These become important, however, if $H^{t}$ is close to the bicritical point.

There is another, rather different limitation to our phenomenological treatment: The superposition approximation, Eqs. (30) and (31), which assumes that the total order parameter distribution $P_{L}(\psi_{x},\psi_{y}$) is a superposition of Gaussians for the AF phase (centered at $\psi_{x}=0$, $\psi_{y}=0$) and of the ordered phase (centered at $\psi^{2}_{x}+\psi^{2}_{y}=\psi^{2}_{\infty}$) with appropriate weights, is not accurate away from the peaks of the distribution. The same applies for the distribution of the order parameter $P(\tilde{m}_{z})$ for the AF phase. In the latter case, the problem is well understood: far from the transition $(H\ll H^{t})$ the distribution $P(\tilde{m}_{z})$ near $\tilde{m}_{z}=0$ is dominated by “slab configurations” where a domain with $\tilde{m}_{z}=-\tilde{m}_{z}^{\rm spont}$ is separated by two domain walls (“antiphase domain boundaries”) from domains with $\tilde{m}_{z}=+\tilde{m}_{z}^{\rm spont}$ (where $\tilde{m}^{\rm spont}_{z}$ is the value of the order parameter where $P(\tilde{m}_{z})$ has its peak). Thus, $P(\tilde{m}_{z})$ is not controlled by a Boltzmann factor containing the volume $L^{3}$, but rather by the surface area $L^{2}$,

[TABLE]

with $f_{\rm int}$ the interfacial excess free energy per unit area (interfacial tension). This “slab configuration” with two planar interfaces is compatible with the periodic boundary conditions, of course.

For a two-component order parameter, the “antiphase domain boundaries” are spread out over the entire volume, and the “phase” of the order parameter gradually rotates from zero to $2\pi$ as the system is traversed. Thus, for $H\gg H^{t}$, we have

[TABLE]

where $\Gamma$ is essentially the “helicity modulus”  helicity . For H near $H^{t}$, however, mixed phase configurations will occur with states $m_{z}\approx{m_{z}^{AF}}$, $\tilde{m}_{z}\approx\pm{\tilde{m}_{z,\infty}}$, $\psi\approx 0$ coexisting with states $m_{z}\approx{m_{z}^{SF}}$, $\tilde{m}_{z}\approx 0$, $|\vec{\psi}|\approx\psi_{\infty}$ with comparable weights. Such mixed phases require more complex “interfaces” in which the order parameters “interpolate” between their coexisting phase values. The generalization of Eqs. (46) and (47) is unknown.

III Monte Carlo Results

The simulations were performed below the bicritical point $T_{b}$ at fixed temperature $T=0.95J/k_{B}$, and we varied the external field $H$ in order to determine the phase transition from AF to SF. When not shown, error bars in the figures showing our results are smaller than the size of the symbols.

The probability distributions of the energy $E$ per site at $T=0.95J/k_{B}$ are shown at the transition field $H^{t}_{L}$ for different lattice sizes $L$ in Fig. 4. For each size we chose the finite size transition field to be located at the point at which the (symmetric) peaks in the probability distribution for the energy were of equal heights. While the dip between the peaks was rather shallow for $L=40$, it rapidly became quite deep for increasing values of $L$ in agreement with the predictions of Binder Binder1982 and Lee and Kosterlitz LeeKosterlitz for a first order transition.

For smaller values of $L$ the probability distribution was almost flat, and the resultant thermodynamic properties showed such substantial finite size rounding that it was not possible to extract useful information about the asymptotic behavior. For this reason, we shall not show raw data for $L<40$ in the figures that follow.

The probability distributions of the z-component for the magnetization $m_{z}$ per site at $T=0.95J/k_{B}$ are shown at the transition field for different lattice sizes $L$ in Fig. 5. For each size we chose the finite size transition field to be located at the same point at which the (symmetric) peaks in the probability distribution for the energy were of equal heights. These data show that the distributions of $m_{z}$ contain two clear, asymmetric peaks of different heights; moreover, these data cannot be described solely by the sum of two Gaussians. Recall that the minimum between the two peaks represents phase coexistence inside the simulation box with one slab in a state $|\tilde{m}_{z}|\approx{\tilde{m}_{z,{\infty}}},\psi=0$ and the other slab having $|\tilde{m}_{z}|\approx 0,|\vec{\psi}|=\psi_{\infty}$ with the two slabs separated by a complex interface connected by the periodic boundary conditions. If we “separate” the distributions into two peaks by choosing the minimum probability as the separation point, we can measure the “weight” of each peak by numerically integrating the probability under each peak. To a high degree of precision, the peaks for each of the lattices sizes, $L$, then have equal weights.

The values of the transition field for each lattice size, $L$, as determined by the equal heights of the two peaks in the probability distributions for the energy, are plotted in Fig. 6. This figure shows very nicely that the variation is linear with $L^{-3}$ for $L\geq 40$. The estimated transition field in the thermodynamic limit is $H^{t}/J=3.83830(5)$. Also shown are the positions of the minima of the 4th order cumulant of the energy. These agree almost perfectly with the values extracted from the locations of the peaks in the probability distributions and will be discussed in more detail later.

From the extrapolation of the peak positions of $P_{L}(m_{z})$ (and of the values of the first moments $\int dm_{z}m_{z}P_{L}(m_{z})$ for each peak) we can estimate $m^{AF}_{z,\infty}$ and $m^{SF}_{z,\infty}$ and hence obtain the difference $\Delta m=m^{SF}_{z,\infty}-m^{AF}_{z,\infty}\approx 0.0352$. Since the factor $\Delta mL^{3}$ exceeds $k_{B}T/J$ by a factor from about $10^{3}$ to $1.8\times 10^{4}$ when $L$ varies from $L=30$ to $L=80$, it is plausible that $H^{t}$ can be located with excellent precision, as shown in Figs. 6 and 7. The fact that all extrapolations, using both quantities from the AF phase and from the SF phase, yield the same transition field $H^{t}$ to very good precision, reinforces the conclusion of Ref. HTL that there is a direct 1st order transition between the two phases with no intervening biconical phase.

As a check on the assumptions about the degeneracy of the SF order parameter, in Fig. 8 we show the order parameter distribution in the xy-plane at the transition field $H^{t}_{L}$ for an $L=60$ lattice. The contours of constant absolute value are almost perfectly circular and show a clear jump from the non-zero value in the SF-phase to a small value in the AF-phase which differs from zero only because of finite size effects.

The variation of the positions of the peaks in the susceptibilities of the uniform magnetization as well as both the z-component of the staggered magnetization and of the SF-order parameter are shown in Fig. 7. Excluding the values for $L=30,40$ as probably being outside the asymptotic region, we fitted the remaining values to obtain an asymptotic value. The positions of all three susceptibilities extrapolate with $L^{-3}$ to a value of $H^{t}/J=3.83830(8)$. This is in perfect agreement with the result of the extrapolations presented in Fig. 6. Using the slope of the susceptibility for $m_{z}$ and Eq. (27) we estimate an effective value $q\approx 3.7$ which is slightly larger than the estimate of $\pi$ as suggested in Sec.II.C.2.

The behavior of the 4th order cumulant of the z-component of the magnetization, $m_{z}$, seen in Fig. 9, shows a clear minimum for each lattice size which sharpens and moves towards larger fields as the lattice size increases. With increasing $L$ the value of the minimum decreases. For $L\geq 50$ the positions of these minima extrapolate linearly with $L^{-3}$ to a value of $H^{t}/J=3.83831(5)$ in excellent agreement with the extrapolated value obtained from the location of the transition from the “equal weight” rule.

The 4th order cumulant of the energy, seen in Fig. 10 also shows a single minimum that sharpens and moves slowly towards higher fields and becomes deeper as the lattice size increases. For $L\geq 40$ the positions of these minima extrapolate linearly with $L^{-3}$ to a value of $H^{t}/J=3.83831(5)$, in excellent agreement with the behavior of the z-component of the magnetization (both the “equal weight” rule and the cumulants).

The probability distributions of the SF order parameter $\psi$ at $T=0.95J/k_{B}$ are shown at the transition field for different lattice sizes $L$ in Fig. 11. For each size we chose the finite size transition field to be located at the same point at which the (symmetric) peaks in the probability distribution for the energy were of equal heights. These data show that the distributions contain two clear peaks at $\pm{\psi_{\infty}}$ corresponding to the SF order and a peak centered about zero corresponding to the AF phase. These peaks cannot be described solely by Gaussians since the states describing phase coexistence (from about $|\vec{\psi}|\approx 0.2$ to $|\vec{\psi}|\approx 0.36$) are not yet strongly suppressed. If we “separate” the distributions into two peaks by choosing the minimum probability as the separation point, we can measure the “weight” of each peak by numerically integrating the probability under each peak. As shown in Table 1 the relative weight of the sum of the “ordering” peaks and the disordered peak depends upon the exact choice of $H^{t}$ and is also slightly dependent upon the choice of $L$. For our best estimate of $H^{t}/J=3.838305$ the value appears to be converging for large $L$ at the estimate of $q=\pi$ that was obtained earlier using the two Gaussian approximation. Note that the result $q=\pi$ was already predicted from the two delta-function distribution appropriate to the thermodynamic limit in Eq. (13).

For completeness we show the probability of the AF order parameter vs the applied field in Fig. 12. For small systems two peaks are seen at $\pm\tilde{m}_{\infty}$ with a broad plateau in between, but as $L$ increases three distinct peaks develop. One peak, centered about $\tilde{m}_{z}=0$ is for the SF phase with no AF order and the peaks at $\pm\tilde{m}_{\infty}$ are for the AF ordered phase. Again, the minima between the peaks for $L\geq 50$ can be interpreted in terms of phase coexistence.

In Fig. 13 we show the variation of the 4th order cumulant of the z-component of the order parameter with field. As the lattice size increases the crossing points move systematically towards slightly larger fields and the values of the cumulant at the crossing points decrease. In fact, the reduction in the value of ${U}_{L}^{z}$ at the crossing seems to accelerate and there is no indication of convergence for the range of lattice sizes studied so from these data alone we cannot tell if the prediction from our simple phenomenological theory is verified. However, in Appendix B we give a tentative interpretation of this behavior in terms of crossover behavior between critical behavior dominated by the bicritical point (prevailing for small $L$) towards first order finite size behavior.

Data for the 4th order cumulant of the xy-component of the order parameter, shown in Fig. 14, reveals similar behavior except that the minimum occurs at fields lower than the field of the crossing point and moves towards higher fields as the size increases. The crossing point values are just below $U_{\infty}^{xy}=0.6$ and are higher than the predicted value of $U_{\infty}^{xy}=0.56056$ assuming an effective value of $q=\pi$. However, the bottom portion of this figure shows a clear tendency for the crossing point values to decrease slightly with increasing size, and the predicted value is not inconsistent with an extrapolation to $L=\infty$. (Note: Use of $U_{\infty}^{xy}=0.6$ in Eq. (22) would yield $q=5$).

The locations of the minima for the cumulants of both order parameters as well as those of the energy and the uniform magnetization $m_{z}$ extrapolated to the thermodynamic limit are also shown in Fig. 15 vs the inverse volume, $L^{-3}$. The observed variation with lattice size agrees with the predictions of our simple double Gaussian theory and, again, we find a common intersection point of $H^{t}/J=3.83830(7)$. The asymptotic size regime, however, appears to begin only for $L\geq 40$. The observation that the depth of the minima decreases very strongly with increasing $L$ is compatible with a scaling as $\approx-L^{3}$ as observed in the Potts model Vollmayr .

Perhaps the most striking results of our study emanate from Eqs. (16), (17), (19) and (20). Plots of ${\langle|\vec{\psi}|\rangle}_{L}$, ${\langle\psi^{2}\rangle}_{L}$, ${\langle|\tilde{m}_{z}|\rangle}_{L}$, and ${\langle\tilde{m}_{z}^{2}\rangle}_{L}$, vs $H$ should show common intersection points for different $L$ at $H^{t}$. These features would $\bf not$ occur at a second order transition, but are nicely consistent with our phenomenological theory. In Figs. 16, 17, and 18 we show data for the first and second moments of $\psi$, $\tilde{m}_{z}$ and, for completeness, $m_{z}$ near the transition. From Eqs. (16) and (17) we conclude that

[TABLE]

and from Eqs. (19) and (20)

[TABLE]

Here too, the data show small but systematic shifts with increasing system size and we can only say that the predictions are consistent with the current data which yield

[TABLE]

These values are in quite reasonable agreement with predictions although more precise values would be needed to draw strong conclusions. However, the discrepancies between the measured and predicted values noted above can probably be attributed to the difference in the location of the intersections and our best estimate for $H^{t}/J=3.83830$. Using more precise data on still larger systems to extrapolate the small finite size variations to $L\to\infty$ could give slightly different estimates than quoted above but would require prohibitively large resources at the present time.

An alternative approach is to base our analysis on Eq. (45) which describes the slope of $\langle\psi^{2}\rangle_{L}$ at the transition field $H^{t}$ as a function of $H/k_{B}T$. From the corresponding plot of $d\langle\psi^{2}\rangle_{L}/d(H/J)|_{H^{t}}$, normalized by $J\Delta{m}L^{3}\psi_{\infty}^{2}/k_{B}T$, vs $L^{-3}$ we expect linear behavior whose intercept is $q/(1+q)^{2}\approx 0.182$ if $q=\pi$. Such a plot, shown in Fig. 19 is, indeed, reasonably compatible with this conjecture.

Near the transition, where $b$ can be neglected, Eq. (42) leads to a simple finite size scaling expression for $\psi^{2}$ that is tested in Fig. 20 in which $\psi^{2}$ is plotted vs. $(H-H^{t})L^{3}$ where $H^{t}/J=3.83830$. For large enough values of $L$ and $H-H^{t}$, the values for $\psi^{2}$ collapse onto a single curve representing a simple analytic scaling function describing the behavior near the transition between these two different ordered phases. For large but negative values of $H-H^{t}$, curves should instead approach the small constant $b$ in Eq. (42) which goes to zero as $L^{-3}$. The data bear out these predictions.

IV Conclusion

In order to provide an understanding of the finite size behavior of a first order transition from a state with simple, discrete degeneracy to a state with an infinite degeneracy, we have performed both theoretical and simulational studies of a uniaxially anisotropic Heisenberg antiferromagnet on finite, simple cubic lattices in an external field $H$ applied along the easy axis.

We first presented a phenomenological theory based upon phase coexistence in the thermodynamic limit with probability distributions of the system in each phase described by delta functions. We hypothesized that the relative weights of the AF and SF phases are $2$ and $2\pi$ by integrating over the angle $\phi$ of the two-component SF order parameter. For finite volume, this description was generalized in terms of suitable Gaussian distributions, which led to the prediction that the moments and cumulants of the order parameters of both phases show common intersection points for large $L$ at the transition field $H^{t}$ (apart from corrections of order $L^{-3}$). The values predicted for these intersections depend upon the effective, relative degeneracy $q=\pi$. We then tested these somewhat speculative predictions via large scale Monte Carlo simulations

We determined the finite size behavior of the model by performing high resolution Monte Carlo simulations at $T=0.95J/k_{B}$. The phase transition can be located quite precisely by using an equal height rule for the probability distribution for the internal energy, and we find that an ”equal weight” rule applies for the order parameter at the transition. The asymptotic finite size behavior of the transition field as well as other quantities does not appear until lattice sizes $L\geq 40$ are reached. In retrospect this is not surprising since the discontinuity in the z-component of the magnetization is quite small, i.e. $\Delta{m_{z}}\leq 0.04$. The locations of the minima in the cumulants for the antiferromagnetic order, spin-flop order, and internal energy extrapolate to the same transition field in the thermodynamic limit as does the “equal weight” rule for the magnetization as predicted by the double Gaussian approximation. Different predictions from the theory yield consistent values for the effective value of the degeneracy $q$ but $\bf only$ for quite large values of $L$. Therefore, we conclude that the simple theory based upon a double Gaussian distribution provides a complete picture of the finite size effects at first order transitions between phases with different symmetries of the order parameters. Since the underlying theory does not depend on the fine characteristics of the model, this means that a heretofore unknown kind of universality at a first order transition has been identified.

Acknowledgements.

DPL thanks the Graduiertenschule MAterials science IN mainZ (MAINZ) and the Alexander von Humboldt Foundation for support. KB thanks the University of Georgia ’s Office of the Vice President for Research at the U. of Georgia for support. This work was supported in part by resources from the Georgia Advanced Computing Resource Center, a partnership between the University of Georgia’s Office of the Vice President for Research and Office of the Vice President for Information Technology.

Appendix A

In this Appendix we want to provide some of the details of the calculations behind the phenomenological theory presented in Sec. II.

The distribution of the order parameter $\vec{\psi}$ in the antiferromagnetic phase, given by Eq. (28) in the main text,

[TABLE]

can be used to evaluate the 4th order cumulant in the following way. The normalization factor $\mathcal{N}$ is given by,

[TABLE]

This factor $\mathcal{N}$ is only applicable for the order parameter moments of a “pure” AF phase and has nothing to do with the normalization factor used to derive Eq. (9). The second moment is simply

[TABLE]

and the fourth moment

[TABLE]

Hence

[TABLE]

Note that for a two-component order parameter, a different normalization of $U_{L}$ would be required for a Gaussian distribution to yield zero in the disordered phase.

In order to obtain the moments of the order parameter of the SF phase using Eq. (29), we have to compute the integrals

[TABLE]

[TABLE]

where the normalization factor $\mathcal{N}$ is given by

[TABLE]

Similarly, writing $\psi^{2}=x$, we conclude

[TABLE]

with negligibly small correction. However, a non-trivial correction term arises in the fourth moment,

[TABLE]

The ratio $U_{L}$ in the ordered SF phase hence becomes

[TABLE]

Invoking the superposition approximation of Eq. (30)- (32) we find that the 4th order cumulant for the SF order then becomes

[TABLE]

Using Eqs. (II.3.5) and  (37) we can straightforwardly obtain the second and fourth moments $\langle\tilde{m}_{z}^{2}\rangle,\langle\tilde{m}_{z}^{4}\rangle$ in both phases and apply a superposition approximation. The resulting cumulant is

[TABLE]

The weighted averages of the order parameter moments in the AF and SF phases become

[TABLE]

and

[TABLE]

A.0.1 Staggered susceptibility maxima

In the absence of symmetry breaking staggered fields, the staggered susceptibility components $\tilde{\chi}^{\prime}_{zz}$ and $\tilde{\chi}^{\prime}_{xy}$, referring to the z-component $\tilde{m}_{z}$ and the xy-components $\vec{\psi}$ of the staggered magnetization, are defined as follows

[TABLE]

[TABLE]

Note that Eq.(65) yields the usual staggered susceptibility in the AF phase, and Eq. (66) in the SF phase. Of course, in the phases where no spontaneous order exists, we simply have from the standard fluctuation relations

[TABLE]

[TABLE]

It is well known that a unique expression for the staggered susceptibility in both phases would require consideration of the limit of an applied staggered field (conjugate to the appropriate order parameter) approaching zero after the thermodynamic limit had been taken. This, however, would be quite inconvenient within the context of simulations.

As already discussed above, the moments $\langle\cdots\rangle_{L}$ are computed from the corresponding moments $\langle\cdots\rangle_{AF},$ $\langle\cdots\rangle_{SF}$ in the pure phases, taking an average with the corresponding weights $a_{AF}$, $1-a_{AF}$. For the second moments, this has already been considered above; here we supply the results for the first moments. For the AF phase we find

[TABLE]

which together with Eq. (65) yields the correct result $\tilde{\chi}^{\prime}_{zz}=\tilde{\chi}^{AF}_{zz}$ in the AF phase. For the SF phase we find

[TABLE]

From Eqs. (30), (II.3.5), (37), (52), (58) (63), and (64) it is straightforward to obtain the expressions for the staggered susceptibilities $k_{B}T\tilde{\chi}^{\prime}_{zz}$, $k_{B}T\tilde{\chi}^{\prime}_{xy}$ defined in Eqs. (65), (66). These expressions show nicely that in the AF phase, i.e. for $\mathcal{Z}\rightarrow-\infty$

[TABLE]

Recall that the AF phase plays the role of the disordered phase for the SF order. Likewise,

[TABLE]

since here the SF phase plays the role of the disordered phase for AF order.

Locating the staggered susceptibility maximum and estimating its height is also an interesting task. Using the abbreviation $q\exp(\mathcal{Z})=Y$, we find from $d(\langle\tilde{m}^{2}_{z}\rangle-\langle|\tilde{m}_{z}|\rangle^{2})/dY=0$ that the maximum occurs for

[TABLE]

and its height is

[TABLE]

This means that the leading correction is reduced by a factor of order $L^{-3/2}$. We also note that exactly at $H=H^{t}$ (i.e., $Y=q$), $\tilde{\chi}^{\prime}_{zz}/\tilde{\chi}^{{}^{\prime},{\rm max}}_{zz}=4q/(1+q)^{2}\approx 0.733$ (if we assume $q=\pi$). Estimation of this ratio offers yet another route to test the value of $q$.

The scaling function for $\tilde{\chi}_{xy}^{{}^{\prime}}$ is given by

[TABLE]

where the expressions $A,B,C$ are given (to the necessary order in inverse powers of $L$) by

[TABLE]

[TABLE]

[TABLE]

So from this expression one can see clearly the “switching” between the two susceptiblities of the “background” phases, which result from this expression when $Y=0$ or $Y=\infty$, respectively, while

[TABLE]

yields the susceptibility maximum.

Thus, to leading order all susceptibility maxima occur at the same location, namely for

[TABLE]

but higher order corrections (of order $L^{-6})$ differ. In both cases, the maximum staggered susceptibility also varies proportional to $L^{3}$ but has a $L^{-3/2}$ correction,

[TABLE]

and $\tilde{\chi}^{{}^{\prime},{\rm max}}_{zz}$ is given by Eq. (74). Here, too, $\tilde{\chi}^{\prime}_{xy}$ for $H=H^{t}$ is smaller than $\tilde{\chi}^{{}^{\prime},{\rm max}}_{xy}$ by the same factor $4q/(1+q)^{2}$ as quoted above.

We now turn to the divergence of the susceptibility $\tilde{\chi}^{SF}_{xy}$ as $L\rightarrow\infty$: Fisher and Privman  FisherPrivman predicted for isotropic, $n$-component magnets that $k_{B}T\tilde{\chi}_{xy}=\langle\psi^{2}\rangle L^{3}=L^{3}\psi^{2}_{\infty}/n+{\rm const}\;L^{2}(n-1)/(n+2)$ while Chen and Landau  ChenLandau predicted $\langle\psi^{2}\rangle L^{3}=L^{3}\psi^{2}_{\infty}/n+{\rm const}\;L^{2}(n-1)/n$ (for $n=3$). Thus, the leading correction is of order $1/L$ (rather than $L^{-3}$, as found in Eq. (70)). Fisher and Privman predict that $L^{3}(\langle\psi^{2}\rangle-\langle|\vec{\psi}|\rangle^{2})$ varies proportional to $(n-1)L$ rather than being the finite constant obtained here. It would be interesting to test these predictions (based on spin wave theory) using suitable numerical results for the present model, but this is a task that must be left for future work.

Appendix B

Here we want to consider the effect of the nearby bicritical point on the finite size behavior of the 4th order cumulant of $\tilde{m}_{z}$ leading to the dramatic variation in crossing points in Fig. 13.

When the normal distance $d_{SF}(T,H)$ from the SF-P phase boundary inside the region of SF order is rather small, the correlation length $\xi_{SF}$ of order parameter fluctuations is very large. As is standard for critical phenomena, this is described by the power law $\xi_{SF}\propto d^{-\nu_{xy}}_{SF}$, where $\nu_{xy}$ is the (universal) critical exponent of the XY model. Likewise, in the region of the AF phase close to the AF-P phase boundary, the correlation length $\xi_{AF}$ of fluctuations of the AF order parameter is very large, $\xi_{AF}\propto d^{-\nu_{I}}_{AF}$, $d_{AF}(T,H)$ being the normal distance from the AF-P phase boundary, and $\nu_{I}$ the critical exponent of the Ising universality class. These power laws, however, apply only for a state point $(T,H)$ that is not too close to the bicritical point $(T_{b},H_{b})$. Close to the bicritical point, all components of the staggered magnetization are simultaneously critical, and the fluctuations are characterized by Heisenberg criticality, with $\xi\propto d^{-\nu_{H}}$, where $d$ is the distance of the state point from the bicritical point. This latter relation applies when $d$ is sufficiently small and outside the two shaded regions in the schematic sketch, Fig. 21. Within the shaded regions, a smooth crossover to the first-order behavior at the spin-flop boundary occurs.

We conclude that near the AF-SF phase boundary for $T<T_{b}(H)$ there is a region where order parameter fluctuations of Heisenberg model type occur, and the correlation length $\xi$ of these fluctuations only gradually diminishes with the distance from the bicritical point. Even at the chosen temperature $T=0.95J/k_{B}$, this correlation length must still be fairly large; and while we have not determined it directly, we can safely conclude this since probability distributions of energy (Fig. 4), magnetization (Fig. 5), and order parameter components are still very broad for $L=40$. Moreover, for smaller $L$ the two coexisting phases can hardly be recognized from these distributions. The broadness of the peaks for $L=40$ in these figures is evidence that the corresponding (staggered) susceptibilities are still very large as well.

From these observations we can also conclude that the behavior of various cumulants at $T=0.95J/k_{B}$ are still Heisenberg-like when $L\ll 40$, and when $L\approx 40$ a gradual crossover from this critical behavior to the behavior characteristic for the first order transition begins to set in.

To provide quantitative evidence for this scenario, we recall that for a Heisenberg antiferromagnet the order parameter distribution $P_{L}({\tilde{\vec{m}}})$ exhibits full rotational symmetry in order parameter space. Using polar coordinates

[TABLE]

we can write P_{L}\Big{(}{\tilde{\vec{m}}}\Big{)}d{\tilde{\vec{m}}}=P_{L}(\tilde{m})\tilde{m}^{2}d\tilde{m}\sin\theta d\theta d\varphi, where only a distribution $P_{L}(\tilde{m})$ of the magnitude $\tilde{m}$ of the order parameter is needed.

The order parameter cumulant of the Heisenberg model

[TABLE]

at criticality is well known 35 ; 36

[TABLE]

But Eq. (83) is not what has been computed in the main text of this paper, where rather cumulants of ${\tilde{\vec{m}}}_{z}$ {Eq. (23)} or the transverse order parameter $\vec{\psi}=({\tilde{\vec{m}}}_{x},{\tilde{\vec{m}}}_{y})$ {Eq. (22)} were considered. However, it turns out that it is straightforward to consider these quantities $U^{z}$, $U^{xy}$ for Heisenberg criticality as well, and actually both of them can be expressed in terms of $U^{H}$.

In order to see this, we first note that

[TABLE]

[TABLE]

the angular part simply cancels out in both expressions. Now from symmetry it is trivial to conclude that

[TABLE]

while in $\langle\tilde{m}^{4}_{z}\rangle$, $\langle\psi^{4}\rangle$ the angular parts contribute, but are straightforward to compute, e.g.

[TABLE]

and similarly

[TABLE]

Thus, we find

[TABLE]

and

[TABLE]

Using Eq. (84), we hence predict that in the bicritical region,

[TABLE]

Looking at Fig. 13, we see that for $L=40$, $L=50$, we still have spurious cumulant crossing somewhere in the region of $U^{z}_{L}\approx 0.34(3)$, while for larger $L$ the crossings move towards significantly smaller values. For the pair $L=(80,100)$ the crossing point is negative. We recall that a degeneracy constant $q=\pi$ (see Eq. (23)) implies $U^{z}_{\infty}=(2-q)/3=-0.3805$. Thus a tentative interpretation of the behavior seen in Fig. 13 is a slow crossover from bicritical behavior to first-order-scaling. With respect to $U^{xy}_{*}$, Fig. 14 did not indicate massive crossover behavior; but this can be understood since $U^{xy}_{*}$ does not differ much from the first-order scaling prediction $U^{xy}_{\infty}\approx 0.56$. Thus, the crossover between these values is “masked” by the standard corrections to finite size scaling.