text:
# Thermodynamics “beyond” local equilibrium

## Abstract

Nonequilibrium thermodynamics has shown its applicability in a wide variety of different situations pertaining to fields such as physics, chemistry, biology, and engineering. As successful as it is, however, its current formulation considers only systems close to equilibrium —those satisfying the so-called local equilibrium hypothesis. Here we show that diffusion processes that occur far away from equilibrium can be viewed as at local equilibrium in a space that includes all the relevant variables in addition to the spatial coordinate. In this way, nonequilibrium thermodynamics can be used and the difficulties and ambiguities associated with the lack of a thermodynamic description disappear. We analyze explicitly the inertial effects in diffusion and outline how the main ideas can be applied to other situations.

Howard Hughes Medical Institute, Department of Molecular Biology, Princeton University, Princeton, New Jersey 08544, USA

Departament de Física Fonamental, Facultat de Física, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Spain

Concepts of everyday use like energy, heat, and temperature have acquired a precise meaning after the development of thermodynamics. Thermodynamics provides us with the basis for understanding how heat and work are related and with the rules that the macroscopic properties of systems at equilibrium follow \[References\]. Outside equilibrium, most of those rules do not apply and the aforementioned quantities cannot unambiguously be defined. There is, however, a natural extension of thermodynamics to systems away from but close to equilibrium. It is based on the local equilibrium hypothesis, which assumes that a system can be viewed as formed of subsystems where the rules of equilibrium thermodynamics apply. Due to the usual disparity between macroscopic and microscopic scales, most systems fall into this category. This is the case of, for instance, the heat transfer from a flame, the flow through a pipe, or the electrical conduction in a wire. Nonequilibrium thermodynamics then extracts the general features, providing laws such as Fourier´s, Fick´s, and Ohm´s, which do not depend on the detailed microscopic nature of the system \[References\].

In contrast, there are other situations where the local equilibrium hypothesis does not hold. Many examples are present in the relaxation of glasses and polymers \[References,References,References\], in the flow of granular media \[References\], and in the dynamics of colloids \[References\]. The main characteristic of such systems is the similarity between microscopic and macroscopic scales, which usually involve internal variables with “slow” relaxation times. The so-called inertial effects in diffusion processes are perhaps the simplest and most illustrative example. In this case, the relaxation of the velocity distribution and changes in density occur at the same time scale. Therefore, local equilibrium is never reached. Here we show how nonequilibrium thermodynamics, as already established in the sixties \[References,References\], can be applied to this situation.

Nonequilibrium thermodynamics \[References\] assumes that the definition of entropy $`S`$ can be extended to systems close to equilibrium. Therefore, entropy changes are given by the Gibbs equation:

$$TdS=dE+pdV\mu dN,$$
(1)
where the thermodynamic extensive variables are the internal energy $`E`$, the volume $`V`$, and the number of particles $`N`$ of the system. The intensive variables (temperature $`T`$, pressure $`p`$, and chemical potential $`\mu `$) are functions of the extensive variables. Local equilibrium means that the Gibbs equation holds for a small region of the space and for changes in the variables that are actually not infinitely slow. Therefore, the internal state of the system has to relax to equilibrium faster than variables change. In this way, all variables retain their usual meanings and the functional dependence between intensive and extensive variables is the same as in equilibrium.

Following this approach, nonequilibrium thermodynamics has been applied to study diffusion processes. The simplest case takes place in one dimension at constant temperature, internal energy, and volume. In this case, from Eq. (1) we obtain a Gibbs equation that depends only on the density and the spatial coordinate $`x`$:

$$Tds(x)=\mu (n,x)dn(x).$$
(2)
Here $`s`$ is the entropy per unit volume and $`n`$ is the density. The chemical potential has the same form as in equilibrium. For instance, for an ideal system —formed of non-interacting particles— it is proportional to the logarithm of the density plus terms that do not depend on the density \[References\]. Notice that these terms can include thermodynamic variables such as temperature and internal energy, and also the spatial coordinate. In the case of non-interacting Brownian particles, its explicit expression is

$$\mu =\frac{k_BT}{m}\mathrm{ln}n+C(x),$$
(3)
where $`m`$ is the mass of the particles, $`k_B`$ the Boltzmann constant, and $`C(x)`$ a function that takes into account possible spatial inhomogeneities. The dynamics of $`n`$ is restricted by the mass conservation law and therefore follows

$$\frac{n}{t}=\frac{J}{x},$$
(4)
with $`J`$ being the flux of mass. An additional assumption of nonequilibrium thermodynamics is that this flux is given by

$$J=L\frac{\mu }{x},$$
(5)
where $`L`$ is the phenomenological coefficient. From this, we obtain the usual diffusion equation

$$\frac{n}{t}=\frac{}{x}\left(D\frac{n}{x}\right),$$
(6)
with the diffusion coefficient $`DL(\mu /n)`$.

When inertial effects are present, changes in density occur at a time scale comparable with the time the velocities of the particles need to relax to equilibrium. The Gibbs equation as stated in Eq. (2) is no longer valid because local equilibrium is never reached. The entropy production depends also on the particular form of the velocity distribution. Both the spatial coordinate, $`x`$, and velocity coordinate, $`v`$, are needed to completely specify the state of the system. Therefore, we consider that local quantities are function of both coordinates. If the system is coupled to other degrees of freedom that relax faster than the velocity and density, a thermodynamic description is still possible. For instance, this is the case of Brownian particles, where the host fluid provides these thermodynamic degrees of freedom. Thus, we consider that diffusion takes place in a two-dimensional space ($`x,v`$) instead of in the original one-dimensional space ($`x`$). In this case, the chemical potential for an ideal system (e.g., non-interacting Brownian particles) is given by

$$\mu =\frac{k_BT}{m}\mathrm{ln}n(x,v)+C(x,v),$$
(7)
where $`C(x,v)`$ is a function that does not depend on the density \[References\]. The form of this function can be obtained by realizing that at equilibrium the chemical potential is equal to an arbitrary constant. If we set this constant equal to zero, we obtain

$$\mu =\frac{k_BT}{m}\mathrm{ln}n+\frac{1}{2}v^2.$$
(8)
Therefore, the Gibbs equation is now

$$Tds(x,v)=\left(\frac{k_BT}{m}\mathrm{ln}n(x,v)+\frac{1}{2}v^2\right)dn(x,v),$$
(9)
The idea of applying the rules of thermodynamics in an internal space was already proposed by Prigogine and Mazur \[References\] and has been used in several situations \[References, References\]. In all of them, however, there was no thermodynamic coupling of these internal degrees of freedom with the spatial coordinate. This is precisely the situation we are considering here.

In the ($`x,v`$)–space the mass conservation law is

$$\frac{n}{t}=\frac{J_x}{x}\frac{J_v}{v}.$$
(10)
Following the standard thermodynamic approach, the flux of mass is given by

$$J_x=L_{xx}\frac{\mu }{x}L_{xv}\frac{\mu }{v},$$
(11)
$$J_v=L_{vx}\frac{\mu }{x}L_{vv}\frac{\mu }{v},$$
(12)
where $`L_{ij}`$, with $`i,j=\{x,v\}`$, are the phenomenological coefficients. There are some restrictions on the values $`L_{ij}`$ can take. Since the system is at local equilibrium in the ($`x,v`$)-space, Onsager relations imply that $`L_{xv}=L_{vx}`$. In addition, the flux of mass in real space, $`\stackrel{~}{J}_x(x)_{\mathrm{}}^{\mathrm{}}vn(x,v)𝑑v`$, has to be recovered from the flux in the ($`x,v`$)-space by contracting the velocity coordinate: $`\stackrel{~}{J}_x(x)=_{\mathrm{}}^{\mathrm{}}J_x(x,v)𝑑v`$. Therefore,

$$_{\mathrm{}}^{\mathrm{}}vn𝑑v=_{\mathrm{}}^{\mathrm{}}\left(L_{xx}\frac{k_BT}{m}\frac{1}{n}\frac{n}{x}+L_{xv}\frac{k_BT}{m}\frac{1}{n}\frac{n}{v}+L_{xv}v\right)𝑑v.$$
(13)
Since $`n(x,v)`$ can take any arbitrary form, the last equality holds if and only if $`L_{xx}=0`$ and $`L_{xv}=n`$. Thus, the only undetermined coefficient is $`L_{vv}`$, which can depend explicitly on $`n`$, $`x`$, and $`v`$.

Previous equations can be rewritten in a more familiar form by identifying the phenomenological coefficients with macroscopic quantities. In this way, with $`L_{vv}=n/\tau `$, the fluxes read

$$J_x=\left(v+\frac{D}{\tau }\frac{}{v}\right)n,$$
(14)
$$J_v=\left(\frac{D}{\tau }\frac{}{x}+\frac{v}{\tau }+\frac{D}{\tau ^2}\frac{}{v}\right)n,$$
(15)
where $`D\frac{k_BT}{m}\tau `$ and $`\tau `$ are the diffusion coefficient and the velocity relaxation time, respectively. The equation for the density is given by

$$\frac{n}{t}=\frac{}{x}vn+\frac{}{v}\left(\frac{v}{\tau }+\frac{D}{\tau ^2}\frac{}{v}\right)n.$$
(16)
This kinetic equation is equivalent to the Fokker-Planck equation for a Brownian particle with inertia since, in an ideal system, the density is proportional to the probability density; i.e. $`n(x,v)=mNP(x,v)`$, where $`P(x,v)`$ is the probability density for a particle to be at $`x`$ with velocity $`v`$, and $`N`$ is the number of particles of the system. The resulting Fokker-Planck equation could have also been derived by following standard techniques of stochastic processes \[References\] or kinetic theory \[References\], which are among the microscopic statistical theories for studying nonequilibrium phenomena.

The approach we have followed, however, explicitly illustrates how thermodynamic concepts can be transferred from equilibrium, through local equilibrium, to far from equilibrium situations. The condition of equilibrium is characterized by the absence of dissipative fluxes ($`J_x=0`$ and $`J_v=0`$). Therefore, from Eq. (14) we obtain that the velocity distribution is Gaussian with variance proportional to the temperature. If deviations from equilibrium are small ($`J_x0`$ and $`J_v=0`$), the local equilibrium hypothesis holds. This is the domain of validity of Fick’s law,

$$J_x=D\frac{n}{x},$$
(17)
which is obtained directly from the equations for the fluxes. In this case, the distribution of velocities is still Gaussian but now centered at a non-zero average velocity and the variance of the distribution is related to the temperature in the same way as in equilibrium. Beyond local equilibrium ($`J_x0`$ and $`J_v0`$), the velocity distribution can take any arbitrary form, from which there is no clear way to assign a temperature. There is, however, a well defined temperature $`T`$: that of local equilibrium in the ($`x,v`$)–space.

In Fig. 1 we illustrate the concepts discussed previously. We show the velocity profiles obtained from Eq. (16) for two representative situations. For fast relaxation of the velocity coordinate, the velocity distribution is Gaussian and centered slightly away from zero, in accordance with local equilibrium concepts. For slow relaxation, however, the velocity distribution loses its symmetry (and its Gaussian form). In this case, the temperature does not give directly the form of the distribution and one has to resort to local equilibrium in the ($`x,v`$)–space to describe the system.

It is important to emphasize that the temperature $`T`$ is the one that enters the total entropy changes and therefore the one related to the second principle of thermodynamics. Other definitions of temperature are possible though. To illustrate this point, let us compute the entropy production $`\sigma `$. This quantity is obtained from local changes in entropy, which are given not only by the production but also by the flow:

$$T\frac{s}{t}=\mu \frac{n}{t}=T\left(\sigma \frac{J_{Sx}}{x}\frac{J_{Sv}}{v}\right),$$
(18)
where $`(J_{Sx},J_{Sv})`$ is the entropy flux. In our case, the expression for the entropy production is

$$\sigma (x,v)=\frac{n(x,v)}{T\tau }\left(v+\frac{k_BT}{m}\frac{\mathrm{ln}n(x,v)}{v}\right)^2.$$
(19)
Now, given a Gaussian velocity distribution $`n(x,v)=n_0(x)e^{mv^2/2k_B\stackrel{~}{T}(x)}`$, we can easily understand the meaning of the temperature $`\stackrel{~}{T}(x)`$ defined through the variance of the distribution: it is the temperature at which the system would be at equilibrium ($`\sigma =0`$). The definition of an effective temperature as that giving zero entropy production can be extended to arbitrary velocity distributions. From Eq. (19), we obtain

$$\frac{1}{\stackrel{~}{T}(x,v)}=\frac{1}{vm}\left(k_B\frac{\mathrm{ln}n(x,v)}{v}\right).$$
(20)
The temperature defined in this way is formally analogous to the equilibrium temperature since the right hand side term of the preceding equation can be rewritten as the derivative of an entropy with respect to an energy:

$$\frac{1}{\stackrel{~}{T}(x,v)}=\frac{s_c(x,v)}{e(v)},$$
(21)
where $`s_c(x,v)=\frac{k_B}{m}\mathrm{ln}n(x,v)`$ and $`e(v)=\frac{1}{2}v^2`$. The term $`s_c(x,v)`$ and $`e(v)`$ can be viewed as the configurational entropy and the kinetic energy per unit of mass, respectively. In general, other definitions of effective temperature are possible. For instance, by considering $`e\left(v\overline{v}(x)\right)`$ instead of $`e(v)`$ in Eq. (21), the resulting temperature would be that of local equilibrium. In this case, however, this temperature does not give zero entropy production but just that of the macroscopic motion. This temperature is then the one at which, once the macroscopic motion is disregarded, the internal configuration of the system would be at equilibrium.

In general, since $`\stackrel{~}{T}(x,v)`$ is not only a function of $`x`$ but also of $`v`$, given a point in space, there is no temperature at which the system would be at equilibrium; i.e. $`\stackrel{~}{T}(x,v)\stackrel{~}{T}(x)`$. If an effective temperature at a point $`x`$ were defined, it would depend on the way the additional coordinate is eliminated. Thus, ambiguities in far from equilibrium quantities arise when considering a lower dimensional space than the one in which the process is actually occurring. This is to some extent similar to what happens with effective temperatures defined through fluctuation-dissipation theorems. In such a case, the effective temperature can depend on the scale of observation \[References\]. It is interesting to point out that all these effective temperatures, despite their possible analogies with the equilibrium temperature, do not have to follow the usual thermodynamic rules since the system is not actually at equilibrium at the temperature $`\stackrel{~}{T}`$.

The idea of increasing the dimensionality of the space were diffusion takes place, so to include as many dimensions as non-equilibrated degrees of freedom the system has, can also be applied to other situations. In a general case, the additional degrees of freedom do not necessarily correspond to the velocity. For instance, let us consider a degree of freedom $`\mathrm{\Theta }(x)`$ that at local equilibrium enters the Gibbs equation in the following way:

$$Tds(x)=\mu dn(x)Bd\mathrm{\Theta }(x),$$
(22)
where $`BB(n,\mathrm{\Theta },T)=T(s/\mathrm{\Theta })_{n,T}`$. In this case, one usually assumes that given $`T`$ , $`n(x)`$, and $`\mathrm{\Theta }(x)`$, the function $`B`$ is completely determined through the equilibrium properties of the system. Far away from equilibrium, we would have to consider explicitly an additional coordinate $`\theta `$, which is related to the degree of freedom by $`\mathrm{\Theta }(x)=\theta n(x,\theta )𝑑\theta `$. The corresponding Gibbs equation

$$Tds(x,\theta )=\mu dn(x,\theta )$$
(23)
would have to take into account the dependence on the coordinate $`\theta `$ through the chemical potential $`\mu `$. Once the Gibbs equation has been obtained, the way to proceed would be analogous to the one we followed for the inertial effects. For instance, some systems with both translational and orientational degrees of freedom can be described by the chemical potential

$$\mu =\frac{k_BT}{m}\mathrm{ln}\frac{n(x,\theta )}{f(\theta )}+U\mathrm{cos}\theta $$
(24)
where $`\theta `$ is now an angular coordinate, $`U\mathrm{cos}\theta `$ is the orientational energy, and $`f(\theta )`$ is a function accounting for the degeneracy of the orientational states (for rotation in three and two dimensions, $`f(\theta )=\mathrm{sin}\theta `$ and $`f(\theta )=1,`$ respectively) \[References\]. This type of systems include, among others, liquid crystals and suspensions of rod-like particles \[References\], field-responsive suspensions \[References\], and polarized systems \[References\]. At local equilibrium some instances of $`B`$ and $`\mathrm{\Theta }`$ are then electric field and polarization; and magnetic field and magnetization. Beyond local equilibrium, by writing the $`(x,\theta )`$ counterpart of Eqs. (10), (11), and (12), one can obtain a kinetic equation that describes the behavior of the system. This equation includes as particular cases the Fokker-Plank equations obtained for those systems by means of microscopic theories \[References,References\].

Along this paper, we have been assuming ideality and locality. The condition of ideality is that the system consists of non-interacting particles. In this case, the chemical potential is proportional to the logarithm of the density plus terms that do not depend on this quantity. Non-ideality can be directly taken into account by considering the right dependence of the thermodynamic quantities on the density and, in general, will give rise to nonlinear partial differential equations. A more difficult aspect to deal with is the presence of non-local effects. In such a case, the interactions between the different parts of the system will need of integro-differential equations to be incorporated in the description.

The main result of our analysis shows that in far from equilibrium diffusion processes, local equilibrium can be recovered when all the relevant degrees of freedom are considered at the same level as the spatial coordinate. In the resulting extended space, thermodynamic quantities, such as temperature and the chemical potential, admit a well defined interpretation. The scheme we have developed may then provide the basis for a consistent formulation of thermodynamics far from equilibrium.

*Acknowledgements*— J.M.R was supported by DGICYT (Spain) Grant No. PB98-1258. J.M.G.V. is an associate of the Howard Hughes Medical Institute.


format: text

text:
# An Unbiased Survey for Outflows in the W3 and W5 Star-Formation Regions

## 1. Introduction

During their birth all stars undergo periods of copious mass loss, frequently characterized by the occurrence of bipolar outflows. These outflows are believed to play a fundamental role in the star formation process. However the exact outflow generating method is obscure at present. To elucidate this problem we are investigating whether the flow properties are correlated over the entire protostellar mass spectrum.

Progress in this area requires that we assemble a statistically valid sample of high-mass outflow systems. This is necessary since existing catalogues of such objects are heterogeneous and statistically incomplete.

## 2. Target Selection and Observations

In order to produce a representative outflow sample a new 12CO survey of the W3 and W5 molecular clouds was initiated at the 14-metre Five College Radio Astronomy Observatory (FCRAO) to look for high velocity gas. W3/5 are intermediate- to high-mass star-forming regions located at a distance of 2 kpc in the Perseus arm.

Survey observations were carried out during April 1999, April 2000, and November 2000. Position-switched 12CO 1–0 emission mapping observations were obtained with beam-width (44′′ at 115 GHz) sampling, using the SEQUOIA spectral plane array.

## 3. Preliminary Results

Since the youngest sources may emit insufficiently in the infrared to be detected by IRAS, the whole data-set was examined for high-velocity gas, indicative of molecular outflows, independent of IRAS data.

Due to the large volume of data an automated potential outflow detection routine was developed. Potential outflow candidates were identified on the basis of line wings present in their spectra. A Gaussian profile was fitted to each spectrum, and then subtracted to leave a residual. The presence of a line wing manifests itself as excess emission (above the noise level) in the residual after subtraction of the Gaussian. (For a more detailed discussion of the potential outflow detection routine see the forthcoming series of papers by Ridge, Bretherton & Moore and Bretherton, Ridge & Moore).

The outflow detection routine successfully identified all known outflows in W3 and W5 (W3-IRS5, W3 (OH), IC1805-W, AFGL 4029 & AFGL 437). In addition it flagged $``$40 regions which might harbour molecular outflows.

Follow-up Nyquist sampled observations of 24 of these outflow candidates were carried out at FCRAO during Spring 2001. Position-switched 13CO observations were also obtained at the central position in order to derive the 12CO optical depth, and to check for multiple-velocity cloud components, which might confuse the outflow detection algorithm. Seven of the potential outflows were discarded on the basis of the 13CO observations.

Figure 1: Contour map of the integrated 12CO 1–0 emission from a new outflow source in W5. (This source was confirmed as a bona fide outflow by observations at the Nobeyama 45m Radio Telescope). This fully-sampled map was made at the FCRAO during the follow-up observations. Spectra corresponding to this region were flagged by the outflow detection routine during reduction of the preliminary programme data. The solid-line contours show the red-lobe emission integrated between -40.8 and -35.8 kms-1. Dashed contours show the blue-lobe emission integrated between -50.8 and -45.8 kms-1. Base contour is 1.8 K kms-1 with contour intervals of 0.5 K kms-1. The triangle indicates the position of an IRAS source, whilst the square indicates the presence of a MSX source.

## 4. Further Work

The next stage of the programme (currently in progress) is the acquisition of high-resolution maps of those outflow candidates not rejected by the FCRAO follow-up observations. High-resolution maps will confirm whether the targets are bona-fide outflows.

Once we have assembled a statistically significant outflow sample we will be able to comment on the applicability of current theories of outflow generating mechanisms. Furthermore the variation of mass and momentum in outflows with source luminosity will be examined, thus enabling stronger constraints to be placed on the actual physics of specific flow models.


text:
# Monte Carlo renormalization group study of the Heisenberg and the XY antiferromagnet on the stacked triangular lattice and the chiral ϕ⁴ model

## I Introduction

The critical behavior of the antiferromagnetic vector spin models on the stacked triangular lattice (STA) is still a controversial issue even after twenty years of extensive studies by means of experimental, field-theoretical, and numerical methods. See  for recent works, and  for a review.

The field-theoretical renormalization group (RG) analysis tells that when the number of spin component $`N`$ is greater than some threshold value $`N_c`$, there are so-called ”chiral fixed points” of the RG which control the critical behavior and are characterized by novel values of critical exponents, while for $`N<N_c`$ such fixed points disappear and the phase transition is of first order (see for a review). Theoretical estimations of $`N_c`$ have been made by various authors by means of $`ϵ`$-expansion , fixed dimensional perturbation , local potential approximation , and effective average action approach , but the estimated values range from negative value to 6.5, depending on the method employed (however, recent results tend to be around 6.0). Thus the critical behaviors of the physical relevant cases $`N=2`$ (XY spin) and $`N=3`$ (Heisenberg spin) are still unclear. A number of experimental studies  and numerical simulations  have yielded results which suggest second-order phase transition for $`N=2`$ and $`3`$. However, recently it has been pointed out that the critical exponent $`\eta `$, calculated from the scaling relation, becomes negative in some of these results, which is unphysical and indicates that the observed critical behavior is in fact pseudo-critical behavior induced by the slow RG flow .

The present work intends to clarify the issue by numerically observing the RG flow and investigating whether there are chiral fixed points or not for several values of $`N`$. The paper is organized as follows: in the next section, we will present the method by which we observe the RG flow in the Monte Carlo simulations. Details of the Monte Carlo simulations are given in Sec. III. In Sec. IV the RG flow diagrams obtained from the simulations for $`N=2,3,8`$ cases are presented. The last section is devoted to the concluding remarks.

## II Numerical observation of the RG flow

### A chiral $`\varphi ^4`$ model

The critical behavior of the $`N`$-component STA model is essentially described by the following $`2N`$-component Ginzburg-Landau-Wilson Hamiltonian :

$$H=𝑑x\left[K\left((\stackrel{}{\varphi }_a)^2+(\stackrel{}{\varphi }_b)^2\right)+r\left(\stackrel{}{\varphi }_a^2+\stackrel{}{\varphi }_b^2\right)+u\left(\stackrel{}{\varphi }_a^2+\stackrel{}{\varphi }_b^2\right)^2+v\left((\stackrel{}{\varphi }_a\stackrel{}{\varphi }_b)^2\stackrel{}{\varphi }_a^2\stackrel{}{\varphi }_b^2\right)\right],$$
(1)
where $`\stackrel{}{\varphi }=(\stackrel{}{\varphi }_a,\stackrel{}{\varphi }_b)`$ is $`N+N`$-component vector field defined on the continuum space. In the field-theoretical studies, $`K`$ is fixed to unity to eliminate ambiguity of coefficients induced by a trivial rescaling of $`\varphi `$, namely, $`\stackrel{}{\varphi }c\stackrel{}{\varphi }`$. In this regularization scheme, $`r`$ plays a role of temperature. In the present work we concentrate on the case $`v0`$. Figure 1 depicts the renormalization group flow of the scaling variables $`r`$, $`u`$, and $`v`$. There are some trivial fixed points, namely:

: high-temperature fixed point $`(r,u,v)=(\mathrm{},0,0),`$

: anisotropic low-temperature fixed point $`(\mathrm{},\mathrm{},\mathrm{}),`$

: isotropic low-temperature fixed point $`(\mathrm{},\mathrm{},0),`$

: Gaussian fixed point $`(0,0,0).`$

The flow on the critical plane projected onto the $`u`$$`v`$ plane is shown in Fig. 2 for (a) $`N<N_c`$ and (b) $`N>N_c`$ cases. When $`N>N_c`$, there are stable and unstable chiral fixed points denoted by $`C_1`$ and $`C_2`$, respectively, beside the $`O(2N)`$ symmetric and the Gaussian fixed points denoted by O and G, respectively. The two fixed points $`C_1`$ and $`C_2`$ approach as $`N`$ decreases, and annihilate each other at $`N=N_c`$.

In the isotropic $`\varphi ^4`$ models, the correction-to-scaling term (distance to the Wilson-Fisher fixed point) rapidly converges to zero as $`b^\omega `$, where $`b`$ is a renormalization factor and $`\omega 0.8`$ is the correction-to-scaling exponent. Therefore the asymptotic critical behavior can be observed in finite systems with moderate size. However, when an anisotropic quartic term is included in the Hamiltonian, it in general induces very slow RG flow along the anisotropy direction. The correction-to-scaling exponent at an anisotropic fixed point (if it exists) is usually of order $`0.1`$ or less, and extremely large system is needed to observe asymptotic critical behavior in the simulation of finite systems. Even when there are no chiral fixed points and the RG flow eventually diverges, the RG flow in the intermediate region is very slow  and it is very hard to observe asymptotic first-order behavior. Thus, in the simulation of finite systems, it is crucial to check the convergence to the fixed point of the renormalized quartic coupling constants $`u`$ and $`v`$. Now let us consider how to observe these renormalized coupling constants in the numerical simulations. In the continuum theory, they describe the behavior of renormalized field variables $`\stackrel{}{\varphi }(𝐱;l)`$ defined as follows:

$$\stackrel{}{\varphi }(𝐱;l)_{|𝐤|<1/l}\mathrm{exp}(i𝐤𝐱)\stackrel{}{\varphi }(𝐤)𝑑𝐤,$$
(2)
where $`l`$ is a cut-off length and $`\stackrel{}{\varphi }(𝐤)`$ denotes Fourier component of the field $`\stackrel{}{\varphi }`$. We denote by $`u_l`$ and $`v_l`$ the renormalized coupling constants of cut-off length $`l`$: the renormalization flow in Fig. 2 is a trajectory of $`u_l`$ and $`v_l`$ on the critical plane when $`l`$ is increased. The renormalization behavior of these quantities may be well observed via their conjugate quantities, which are easier to observe in numerical simulations:

$$(\stackrel{}{\varphi }^2(𝐱;l))^2,$$
(3)
$$\stackrel{}{\varphi }_a^2(𝐱;l)\stackrel{}{\varphi }_b^2(𝐱;l)(\stackrel{}{\varphi }_a(𝐱;l)\stackrel{}{\varphi }_b(𝐱;l))^2.$$
(4)
Of course, these quantities are affected not only by $`u_l`$ and $`v_l`$ but also by other irrelevant scaling fields, such as coefficients of higher order terms like $`|\stackrel{}{\varphi }|^6`$ and $`|\stackrel{}{\varphi }|^8`$. However, these irrelevant scaling fields rapidly vanishes (faster than the leading correction-to-scaling term in the isotropic model, $`l^{0.8}`$) and their effects can be ignored after sufficient renormalization as long as $`u_l`$ and $`v_l`$ converge (or diverge) very slowly.

Note that the definition of the renormalized field $`\stackrel{}{\varphi }(𝐱;l)`$ needs scaling prefactor so that the coefficients such as $`u_l`$ and $`v_l`$ remain finite when $`l`$ goes to infinity. The following scaling scheme is simple and suitable for numerical simulations :

$$\stackrel{}{\varphi }^{}(𝐱;l)\frac{\stackrel{}{\varphi }(𝐱;l)}{\sqrt{\stackrel{}{\varphi }(𝐱;l)^2}}.$$
(5)
Substituting $`\stackrel{}{\varphi }(𝐱;l)`$ by $`\stackrel{}{\varphi }^{}(𝐱;l)`$, equation (3) and (4) lead as follows:

$$\frac{(\stackrel{}{\varphi }^2(𝐱;l))^2}{(\stackrel{}{\varphi }^2(𝐱;l))^2},$$
(6)
$$\frac{\stackrel{}{\varphi }_a^2(𝐱;l)\stackrel{}{\varphi }_b^2(𝐱;l)(\stackrel{}{\varphi }_a(𝐱;l)\stackrel{}{\varphi }_b(𝐱;l))^2}{(\stackrel{}{\varphi }^2(𝐱;l))^2}.$$
(7)
Basically, in the present work, we observe these quantities in the numerical simulations of finite lattice version of the Hamiltonian (1) defined as follows:

$$H=\frac{KN}{2}\underset{<ij>}{}(\stackrel{}{\varphi }(i)\stackrel{}{\varphi }(j))^2+r\underset{i}{}\stackrel{}{\varphi }(i)^2+u\underset{i}{}\left(\stackrel{}{\varphi }(i)^2\right)^2+v\underset{i}{}\left[(\stackrel{}{\varphi }_a(i)\stackrel{}{\varphi }_b(i))^2\stackrel{}{\varphi }_a(i)^2\stackrel{}{\varphi }_b(i)^2\right],$$
(8)
where $`\stackrel{}{\varphi }(i)=(\stackrel{}{\varphi }_a(i),\stackrel{}{\varphi }_b(i))`$ is an $`N+N`$-component vector defined on the lattice site i of an $`L\times L\times L`$ cubic lattice with periodic boundary condition being imposed, the summation $`_{<ij>}`$ runs over all nearest neighbor pairs of the lattice site. We use a restriction $`u+v/42r=0`$ so that the minimum of the Hamiltonian takes place at $`|\stackrel{}{\varphi }(i)|=1`$ to eliminate the ambiguity of the trivial rescaling of $`\varphi `$. Actually we use the following parameterization:

$$r=2\lambda ,u=\lambda (1+A),v=4\lambda A$$
(9)
where $`A`$ controls the strength of the anisotropy and $`\lambda `$ controls the ”hardness” of the spin— the larger $`\lambda `$ is, the smaller the fluctuation of $`|\stackrel{}{\varphi }(i)|`$ is. A line $`\lambda =0`$ corresponds to $`uv/4=0`$, on which instability of the Hamiltonian (8) occurs; thus $`\lambda `$ should be positive.

In the lattice models, renormalized field corresponds to block spin variables. For example, consider a block spin variable $`\stackrel{}{\varphi }_{b,L}`$ defined on a block of $`b\times b\times b`$ spins in an $`L\times L\times L`$ lattice model. One can calculate the quantities (6) and (7) using $`\stackrel{}{\varphi }_{b,L}`$, which we denote by $`U_{b,L}`$ and $`V_{b,L}`$, respectively, and check the convergence or divergence of the renormalization flow. In the conventional MCRG scheme , one calculates (6) and (7) for fixed value of $`L`$ and increasing value of $`b`$. Note that, after the ”renormalization” by the block spin transformation, the number of spins, $`(L/b)^3`$, decreases and the Hamiltonian does not stay in the same Hamiltonian space, unlike the continuum case. Therefore the number of block spins must be sufficiently large, i.e. $`L/b1`$, so that (6) and (7) take the same value at the fixed point for different values of $`b`$.

The present author proposed another scheme , in which the ratio $`L/b`$ is fixed and both $`L`$ and $`b`$ are increased. Let us denote a renormalized Hamiltonian which contains $`l^3`$ block spins by $`H_l(K,r,u,v,\{w_i\})`$, where $`K,r,u,v`$ are the renormalized coefficients of terms which correspond to that in the original Hamiltonian (8), and $`\{w_i\}`$ denote set of coefficients of higher order irrelevant terms such as $`|\stackrel{}{\varphi }|^6`$, which is not included in the original Hamiltonian. Now consider the behavior of $`\stackrel{}{\varphi }_{L,L}`$ and $`\stackrel{}{\varphi }_{2L,2L}`$: the behavior of both quantities are described by a renormalized Hamiltonian which contains only 1 spin, but with different coefficients. Here, $`\stackrel{}{\varphi }_{2L,2L}`$ is a renormalized spin, obtained by applying real-space renormalization of factor-$`2L`$ to the original spins. Let us regard this factor-$`2L`$ renormalization as two successive renormalizations, namely factor-$`2`$ renormalization at first and then factor-$`L`$ renormalization. Accordingly, the original Hamiltonian $`H_{2L}(K,r,u,v,\{0\})`$ is at first renormalized to $`H_L(K^{},r^{},u^{},v^{},\{w_i\})`$, then again renormalized to $`H_1(K^{\prime \prime },r^{\prime \prime },u^{\prime \prime },v^{\prime \prime },\{w_i^{}\})`$. If the starting Hamiltonian is at a fixed point and $`L`$ is sufficiently large, coefficients before and after the factor-$`2`$ renormalization should remain unchanged except for the higher order coefficients $`\{w_i^{}\}`$ . At this stage, $`U_{2L,2L}`$ and $`V_{2L,2L}`$ can be calculated (in other words, factor-$`L`$ renormalization) from $`H_L(K^{},r^{},u^{},v^{},\{w_i\})`$, while $`U_{L,L}`$ and $`V_{L,L}`$ are calculated from $`H_L(K,r,u,v,\{0\})`$. Since the coefficients of these two Hamiltonians are the same (except the irrelevant ones), $`U_{2L,2L}=U_{L,L}+O(L^\omega )`$ and $`V_{2L,2L}=V_{L,L}+O(L^\omega )`$ should be satisfied, where $`\omega `$ denotes correction-to-scaling exponent of the higher order terms such as $`|\stackrel{}{\varphi }|^6`$. The so-called ”phenomenological renormalization group method”  uses this property to locate the critical point.

When the Hamiltonian is not at a fixed point, the coefficients before and after the first factor-$`2`$ renormalization differ: they moves along the renormalization flow. This difference is reflected by $`U_{L,L}`$ and $`V_{L,L}`$; If the renormalization flow slowly converges to some fixed point, $`U_{L,L}`$ and $`V_{L,L}`$ slowly converge to some fixed value, while when the phase transition is of first-order and the renormalization flow diverges, so do $`U_{L,L}`$ and $`V_{L,L}`$ . We will henceforth denote these quantities simply by $`U_L`$ and $`V_L`$, and investigate their behavior at the critical point when $`L`$ is increased. Here their explicit form is written down for the readers’ convenience:

$$U_L=\frac{\left(\stackrel{}{M}^2\right)^2}{\stackrel{}{M}^2^2},$$
(10)
$$V_L=\frac{\stackrel{}{M}_a^2\stackrel{}{M}_b^2(\stackrel{}{M}_a\stackrel{}{M}_b)^2}{\stackrel{}{M}^2^2},$$
(11)
where $`\stackrel{}{M}=(\stackrel{}{M}_a,\stackrel{}{M}_b)=_i\stackrel{}{\varphi }(i)`$. In the past numerical studies in which only a specific model was investigated, the only way to reach the final fixed point (if it exists) or to observe asymptotic first-order behavior was to simulate extremely large systems, and it was impossible in general. In the present work, we investigate various Hamiltonians and scan the Hamiltonian space in search of fixed points. This method allows us to determine the asymptotic critical behavior by simulations of moderately large systems.

Note that, in the theoretical works one can concentrate on the RG flow on the critical plane on which inverse susceptibility vanishes, while in the numerical works one must determine the critical point from the numerical data. For this purpose we use a quantity which corresponds to the following quantity in the continuum theory:

$$\frac{_{|k|<1/l}(𝐤/l)^2\stackrel{}{\varphi }^2(𝐤)d^D𝐤}{(\stackrel{}{\varphi }^2(𝐱;l))}=l^2\frac{(\stackrel{}{\varphi }(𝐱;l))^2}{(\stackrel{}{\varphi }^2(𝐱;l))}.$$
(12)
In the lattice models, the above quantity corresponds to the normalized correlation between two adjacent block spins. Actually we observe the following quantity:

$$C_L=\frac{\stackrel{}{\varphi }(𝐤_1)\stackrel{}{\varphi }(𝐤_1)}{\stackrel{}{M}^2},$$
(13)
where $`\stackrel{}{\varphi }(𝐤_\mathrm{𝟏})=_𝐫\stackrel{}{\varphi }(r)\mathrm{exp}(i𝐤_\mathrm{𝟏}𝐫)`$ with $`𝐤_1=(2\pi /L,0,0)`$.

Figure 3 depicts the RG flow in $`U_L`$$`V_L`$$`C_L`$ space: a number of trivial fixed points, namely, high-temperature fixed point $`(1+1/2N,\frac{(2N+1)(N1)}{8N(N+1)},1)`$, isotropic low-temperature fixed point $`(1,\frac{N1}{4(N+1)},0)`$, and anisotropic low-temperature fixed point $`(1,\frac{1}{4},0)`$ are denoted by H, $`L_0`$, and $`L_1`$, respectively. Beside these points, there are some trivial manifolds: when the probability distribution of the order parameter $`\stackrel{}{M}`$ is isotropic, a relation $`V_L=\frac{N1}{4(N+1)}U_L`$ is satisfied, while in the strong anisotropy limit where $`\stackrel{}{M}_a\stackrel{}{M}_b=0`$ and $`|\stackrel{}{M}_a|=|\stackrel{}{M}_b|`$ hold, $`V_L=U_L/4`$ is satisfied. At the Gaussian fixed point, the behavior of the finite system is governed by the zero-mode , therefore the Gaussian fixed point lies on the $`C_L=0`$ plane. The RG flow on the critical plane, projected onto the $`U_L`$$`V_L`$ plane, is shown in Fig. 4 for (a) $`N<N_c`$ and (b) $`N>N_c`$ cases. In both cases, the RG flow along the (approximately) horizontal direction rapidly converges, while along the (approximately) vertical direction the flow is expected to slowly converge/diverge.

To obtain the RG flow on the critical plane, we investigate the Hamiltonian (8) for various values of $`A`$ and $`\lambda `$, tuning $`K`$ so that $`C_L=C_{2L}`$ is satisfied and observe the difference between $`(U_L,V_L)`$ and $`(U_{2L},V_{2L})`$. It should be noted that the above definition of critical plane induces systematic deviation. For example, consider the RG flow of the isotropic $`\varphi ^4`$ model (set $`v=0`$ in (8) ). The flow of $`(U_L,C_L)`$ near the Wilson-Fisher fixed point is depicted in Fig. 5 (a) . If one define the critical point by $`C_L(K)=C_{2L}(K)`$, the flow of $`U_L`$ becomes the one shown in Fig. 5(b), thus the arrow tends to “overshoot” the RG fixed point. However, the direction of the flow can be correctly estimated, and the systematic error vanishes as one approaches the fixed point.

### B STA models

We also observe the RG flow of the following STA model:

$$H=K\underset{<ij>}{}J_{ij}\stackrel{}{S_i}\stackrel{}{S_j},$$
(14)
where $`\stackrel{}{S_i}`$ denotes a two-component (XY) or three-component (Heisenberg) vector spin with $`|\stackrel{}{S}_i|=1`$ defined on the lattice site $`i`$ of the stacked triangular lattice. We set $`J_{ij}=1`$ (antiferromagnetic) for intra-plane nearest neighbor pairs, $`J_{ij}=3/4`$ (ferromagnetic) for inter-plane nearest neighbor pairs, and $`J_{ij}=0`$ otherwise, so that the Fourier transform of $`J_{ij}`$ near the antiferromagnetic mode  $`\stackrel{}{k}_{AF}=(2\pi /3,0,0)`$ becomes isotropic in the $`\stackrel{}{k}`$-space, i.e. $`J(\stackrel{}{k}_{AF}+\stackrel{}{k})J(\stackrel{}{k}_{AF})+c|\stackrel{}{k}|^2+O(|\stackrel{}{k}|^4)`$. Note that the critical values of quantities such as $`U_L,V_L,`$ and $`C_L`$ do not depend on the microscopic lattice structure, but do depend on the macroscopic lattice structure such as boundary conditions and aspect ratios . Therefore we use rectangular system which contains $`L_x\times L_y\times L_z`$ spins (see Fig. 6), imposing periodic boundary condition for all three directions. The aspect ratio then becomes $`L_x:\sqrt{3}L_y/2:L_z`$. We simulate the system with $`(L_x,L_y,L_z)=(21,24,21)`$, $`(42,48,42)`$, $`(84,96,84)`$, and $`(126,144,126)`$, all of which give aspect ratio $`1:0.99\mathrm{}:1`$.

The order parameter is defined as follows :

$$\stackrel{}{M}_a=\text{Re}[\stackrel{}{S}(\stackrel{}{k}_{AF})]=\stackrel{}{M}_A\frac{1}{2}\stackrel{}{M}_B\frac{1}{2}\stackrel{}{M}_C,$$
(15)
$$\stackrel{}{M}_b=\text{Im}[\stackrel{}{S}(\stackrel{}{k}_{AF})]=\frac{\sqrt{3}}{2}\stackrel{}{M}_B\frac{\sqrt{3}}{2}\stackrel{}{M}_C,$$
(16)
where $`\stackrel{}{S}(\stackrel{}{k}_{AF})`$ denotes Fourier component of $`\stackrel{}{S}_i`$ at $`\stackrel{}{k}_{AF}`$, and $`\stackrel{}{M}_A,\stackrel{}{M}_B,\stackrel{}{M}_C`$ denote the magnetization on the three sublattices. Then $`U_L`$ and $`V_L`$ are calculated from Eq. (10) and Eq. (11). $`C_L`$ is calculated as follows:

$$C_L=\frac{|\stackrel{}{S}(\stackrel{}{k}_{AF}+\stackrel{}{k}_1)|^2+|\stackrel{}{S}(\stackrel{}{k}_{AF}\stackrel{}{k}_1)|^2}{\stackrel{}{S}^2(0)},$$
(17)
where $`\stackrel{}{k}_1`$ is the smallest, non-zero momentum in the $`\stackrel{}{k}`$-space. We observed the values of $`C_L`$ for three kinds of direction, namely, $`\stackrel{}{k}_1=(\frac{2\pi }{L_x},0,0)`$, $`(0,\frac{4\pi }{\sqrt{3}L_y},0)`$, and $`(0,0,\frac{2\pi }{L_z})`$, and confirmed that they coincide each other within statistical errors. This indicates that the simulated system is spatially isotropic.

## III Monte Carlo simulation

### A chiral $`\varphi ^4`$ models

We used the usual single-spin update Metropolis algorithm, since there is no spin reflection axis which preserves the anisotropy term in (8) and the cluster algorithm  cannot be applied. In a single-spin update process, the new spin value was chosen as follows:

$$\stackrel{}{\varphi }_{new}(i)=\stackrel{}{\varphi }_{old}(i)+\frac{R_g}{\sqrt{2N}}\stackrel{}{G}$$
(18)
where $`\stackrel{}{G}`$ is an $`2N`$-component vector whose components are independent Gaussian random variables with average $`0`$ and variance $`1`$. The optimum value of the amplitude of random move, $`R_g`$, was determined by simulating small systems. The optimum value was found to be $`R_g0.3`$ for all $`N`$, $`A,`$ and $`\lambda `$. Each Metropolis sweep was followed by one overrelaxation-type update . Hereafter we refer one Metropolis sweep plus one overrelaxation-type sweep as “MCS”.

Physical quantities were observed for every $`C_{INT}\times L^2`$ MCS, where $`C_{INT}`$ is a constant which is adjusted so that the correlations between successively observed quantities become less than $`0.70`$, namely:

$$\frac{X_tX_{t+1}X_t^2}{X_t^2X_t^2}<0.7,$$
(19)
where $`X_t`$ denotes observed value of a physical quantity $`X`$ at the $`t`$th observation. The value of $`C_{INT}`$ varies from $`0.1`$ to $`2.0`$, depending mainly on $`A`$. Actual values will be summarized in the later sections. For each values of $`L`$, $`N`$, $`A`$, and $`\lambda `$, physical quantities were observed $`10^4`$ to $`10^5`$ times, depending on the accuracy required to observe the RG flow. Statistical errors were estimated by the Jackknife method. Histogram method  was used to calculate thermal averages at $`K`$ slightly away from that where the simulation is actually carried out.

Simulations were mainly carried out on Fujitsu VPP5000 vector processors at JAERI. The checkerboard-wise decomposition of the lattice was used for vectorization. Furthermore, we simultaneously simulated a number of independent systems, whose spins were stored into one long vector. This promoted the vectorization, especially when $`L`$ is small, and accelerated the calculations. For $`N=8`$ and $`L=16`$ system, one MCS took 0.84 ms on the VPP5000.

### B STA models

In the frustrated spin models, the cluster update algorithm tends to flip all the spins and does not accelerate the simulation. Therefore we again used only the single-spin update Metropolis algorithm, followed by an overrelaxation sweep, i.e. 180 degree rotation of a spin $`\stackrel{}{S}_i`$ with respect to its local field $`_jJ_{ij}\stackrel{}{S}_j`$. In the Metropolis update, the new spin direction was chosen as follows:

$$S_x=\sqrt{1z^2}\mathrm{cos}(\theta ),S_y=\sqrt{1z^2}\mathrm{sin}(\theta ),S_z=z,$$
(20)
for Heisenberg case and

$$S_x=\mathrm{cos}(\theta ),S_y=\mathrm{sin}(\theta ),$$
(21)
for XY case, where $`z`$ and $`\theta `$ are uniformly distributed random numbers in the range $`[1,1)`$ and $`[0,2\pi )`$, respectively. Physical quantities were observed for every $`L_x\times L_y/20`$ MCS for the Heisenberg case and $`L_x\times L_y/10`$ MCS for the XY case. This interval was long enough to satisfy the condition (19) as long as the energy histogram has no double-peak. Physical quantities were observed $`10^4`$ times for all system sizes up to $`(L_x,L_y,L_z)=(42,48,42)`$.

To observe the energy histogram of the larger system, we also simulated very large systems $`(L_x,L_y,L_z)=(84,96,84)`$ and $`(126,144,126)`$ for the XY case. Owing to the prominent double-peak of the energy histogram, autocorrelation time was much longer than $`L_x\times L_y/10`$ MCS in these sizes. We used $`4\times 10^6`$ MCS for the $`(84,96,84)`$ case and $`7\times 10^6`$ MCS for the $`(126,144,126)`$ case, which were 1000 times longer than the auto-correlation time. One MCS of $`N=3`$ and $`(L_x,L_y,L_z)=(84,96,84)`$ system took 27 ms on the VPP5000.

### C Stiefel’s $`V_{N,2}`$ models

Stiefel’s $`V_{N,2}`$ model  corresponds to the Hamiltonian (8) with $`\lambda =\mathrm{}`$ and $`A=\mathrm{}`$, in which restrictions $`\varphi _a\varphi _b=0`$, $`\varphi _a^2=\varphi _b^2=1/2`$ are imposed to all the spins. We also simulated this model and investigated the RG flow. In the Monte Carlo simulation, only the single-spin Metropolis update was used. We only simulated the $`V_{3,2}`$ model, since the $`V_{2,2}`$ model is known to exhibit strong first-order behavior . In the Metropolis update, the new spin value was chosen as follows:

$`\stackrel{}{\varphi }_a`$ $`=`$ $`(z,\sqrt{1z^2}\mathrm{cos}\theta _1,\sqrt{1z^2}\mathrm{sin}\theta _1)/\sqrt{2},`$ (22)
$`\stackrel{}{\varphi }_b`$ $`=`$ $`(\sqrt{1z^2}\mathrm{sin}\theta _2,\mathrm{sin}\theta _1\mathrm{cos}\theta _2z\mathrm{cos}\theta _1\mathrm{sin}\theta _2,\mathrm{cos}\theta _1\mathrm{cos}\theta _2z\mathrm{sin}\theta _1\mathrm{sin}\theta _2)/\sqrt{2}`$ (23)
where $`z`$ and $`\theta _{1,2}`$ are uniformly distributed random numbers in the range $`[1,1)`$ and $`[0,2\pi )`$, respectively. Physical quantities were observed for every $`L^2`$ Metropolis sweeps. This interval was long enough to satisfy the condition (19) for all system sizes up to $`L=48`$. Physical quantities were observed $`10^4`$ times for these system sizes.

To observe the energy histogram of the larger system, we also simulated a very large system $`L=64`$ and $`L=80`$. Like the STA-XY case, a prominent double-peak of the energy histogram emerges and autocorrelation time gets much longer than $`L^2`$ MCS in these sizes. We used $`6\times 10^6`$ MCS for the $`L=64`$ case and $`2\times 10^7`$ MCS for the $`L=80`$ case, which were about 1000 times longer than the auto-correlation time. For the largest size $`L=80`$, one Metropolis sweep took 22 ms on the VPP5000.

## IV results

### A XY case

Parameters $`\lambda `$, $`A`$, and $`K`$ at which simulations were carried out are summarized in Table I. Figure 7 shows the histogram of energy per spin, $`E=J_{ij}\stackrel{}{S}_i\stackrel{}{S}_j/L_xL_yL_z`$ of the STA-XY model with $`(L_x,L_y,L_z)=(126,144,126)`$ and $`(84,96,84)`$ at $`K=0.77262`$. It can be seen that the valley between the two peaks deepens as the system size increases. This is a clear evidence of first-order transition. This peak does not appear when the size is smaller, therefore the double peak has not been observed in the past studies of STA-XY model in which smaller systems have been used. Note that the first-order behavior has already been observed in the Monte Carlo simulations of other models which possess the same symmetry as the STA-XY model, such as Stiefel model and restricted STA model , and quasi-one dimensional STA model . These models have stronger anisotropy than the STA-XY model and the first-order behavior can be observed in smaller systems.

The two peaks of the histogram correspond to ordered and disordered phases, and the energy difference between the peaks provides lower bound of the latent heat per spin $`\mathrm{\Delta }`$. On the other hand, experimental measurement of specific heat of the STA-XY materials such as CsMnBr<sub>3</sub> or Holmium  indicates second-order transition or very weak first-order transition whose latent heat is beneath the resolution of measurements, thus the upper bound of the latent heat is given. To compare the strength of the first-order transition, an adimensional parameter $`c=\mathrm{\Delta }/k_BT_c`$ is usually used, where $`k_B`$ denotes Boltzmann constant and $`T_c`$ denotes the transition temperature. From Fig. 7, the upper bound of the latent heat per spin is estimated as $`\mathrm{\Delta }>7.8\times 10^3`$ and $`c>6.0\times 10^3`$, while Ref. provides estimates $`c<2.9\times 10^5`$ for CsMnBr<sub>3</sub> and $`c<2.8\times 10^5`$ for Holmium. These upper bound of $`c`$ obtained from experiments is far smaller than the lower bound of $`c`$ obtained in the present work, thus there is an apparent contradiction between them.

Note that the value of $`c`$ may depends on the strength of the easy-plane anisotropy of the spin and the ratio between inter and intra layer coupling constants. The result of Ref.  indicates that when the inter-layer coupling is much larger than the intra-layer coupling, the first-order transition becomes stronger. On the other hand, a weak easy-plane anisotropy will lead to Heisenberg-to-XY crossover which weaken the transition, since the transition is weakened as the number of spin component is increased and eventually becomes second-order. Actually, we will see in the following section that the STA-Heisenberg model does not exhibit as strong a transition as the STA-XY model. Thus the most plausible explanation of the apparent discrepancy between the experimental results and the present work may be that the finiteness of the easy-plane anisotropy in the experimental materials induces Heisenberg-to-XY crossover and the transition is weakened.

Figure 8 and 9 show the RG flow of $`(U_L,V_L)`$ at the strong and weak anisotropy region, respectively, of the $`N=2`$ chiral $`\varphi ^4`$ model. Figure 8 also shows the same plot for the STA-XY model. One can see that the STA-XY model is on the “runaway” trajectory, which is another evidence of first-order transition. It can also be seen that there are no stable fixed points other than the $`O(4)`$ symmetric one. This means that any models or materials which possess the same symmetry as STA-XY model, namely $`Z_3\times O(2)`$, exhibit first-order transition.

Fig. 10 shows the vertical velocity of the RG flow $`V_{2L}V_L`$ plotted against $`V_L`$. The plots for the two different sizes $`L=6`$ and $`L=8`$ coincide within statistical errors. This indicates that the system size is large enough to eliminate finite-size artifact of the RG flow. Thus we can safely extrapolate the asymptotic critical behavior from the finite-size RG flow. Moreover, the plots of the chiral $`\varphi ^4`$ model and the STA-XY model are essentially the same. This means that the critical behavior of the STA-XY model is well reproduced by the chiral $`\varphi ^4`$ model.

### B Heisenberg case

Parameters $`\lambda `$, $`A`$, and $`K`$ at which simulations were carried out are summarized in Table II. Unlike the STA-XY case, the STA-Heisenberg model does not show the double-peak behavior within the simulated size $`(L_x,L_y,L_z)=(84,96,84)`$. However, the $`V_{3,2}`$ model with $`L=80`$ shows a clear double-peak behavior. Figure 11 shows the histogram of energy per spin, $`E=3L^3J_{ij}(\stackrel{}{S}_i\stackrel{}{S}_j)^2`$ at $`K=0.217685`$. Although it is necessary to simulate larger system and show that the valley between the two peaks deepens to confirm the first-order transition, this double-peak strongly indicate that the transition is of first-order.

Figure 12 and 13 show the RG flow of $`(U_L,V_L)`$ at the strong and weak anisotropy region, respectively, of the $`N=3`$ chiral $`\varphi ^4`$ model. Figure 12 also shows the same plot for the STA-Heisenberg model and the $`V_{3,2}`$ model. It can be seen that the STA-Heisenberg model and the $`V_{3,2}`$ model are on the “runaway” trajectory, therefore undergo first-order transition in the thermodynamic limit. It can also be seen that there are no fixed points other than the $`O(6)`$ symmetric one. This means that any models or materials which possess the same symmetry as STA-Heisenberg model, namely $`Z_3\times O(3)`$, exhibit first-order transition. From Fig. 12, we roughly estimate that $`L=96`$ STA-Heisenberg model needs factor-2 renormalizations three or more times to reach the strong first-order region, which means that the clear first-order behavior will be observed in an $`L800`$ system.

Fig. 14 shows the vertical velocity of the RG flow $`V_{2L}V_L`$ plotted against $`V_L`$. In this case, this quantity depends not only on $`V_L`$ but also on $`U_L`$, although the dependence is not so strong. Therefore the plot can not be fitted to a single line well. However, there is no systematic deviation between the plots of the two different sizes. Thus we can safely extrapolate the asymptotic critical behavior from the finite-size RG flow. As in the XY case, the plots of the chiral $`\varphi ^4`$ model and the STA-Heisenberg model are essentially the same and the critical behavior of the STA-Heisenberg model is well reproduced by the chiral $`\varphi ^4`$ model.

### C $`N=8`$ case

It is necessary to demonstrate that the MCRG scheme presented in Sec. II is able to detect the chiral fixed point if it does exists. For this purpose, we investigate $`N=8`$ case, which is greater than any theoretical estimate of $`N_c`$ in literature . Parameters $`\lambda `$, $`A`$, and $`K`$ at which simulations were carried out are summarized in Table III. Figure 15 shows the RG flow of $`(U_L,V_L)`$ for $`L=816`$. The stable chiral fixed point $`C_1`$ and the unstable chiral fixed point $`C_2`$ are clearly seen.

Fig. 16 shows the vertical velocity of the RG flow $`V_{2L}V_L`$ plotted against $`V_L`$, using only the data near the stable fixed point $`C_1`$. Unlike the $`N=2`$ and the $`N=3`$ case, this plot has a zero. The plots of the two different size coincide within statistical error, which indicates that the simulated system is large enough to observe the RG flow.

## V conclusion

The large scale Monte Carlo simulation of the STA-XY model has revealed that the transition is of first-order. This is the first numerical work to confirm the first-order transition of the simple and isotropic STA-XY model. The result also indicate that the transition is strong enough to be observed in experiments. But the finiteness of the easy-plane spin anisotropy of the actual materials will weaken the transition compared to the STA-XY model simulated in the present work, in which the easy-plane anisotropy is infinitely strong.

We have also investigated the RG flow the $`2N`$–component chiral $`\varphi ^4`$ model numerically and have found that there are no chiral fixed points in $`N=2`$ and $`N=3`$ cases, and we conclude that the STA-Heisenberg model also undergoes first-order transition. By a rough estimate, we expect that the STA-Heisenberg model exhibits strong first-order behavior when the size is larger than $`L=800`$. For $`N=8`$ case, we have found the chiral fixed points, thus the value of $`N_c`$ is estimated as $`3<N_c<8`$. Further study to determine the precise value of $`N_c`$ may be of some interest, because such a study will serve to test the accuracy of field theoretical techniques. However, once the nature of physical relevant case $`N=3`$ and $`N=2`$ is clarified, the precise value of $`N_c`$ is of less importance.

## VI Acknowledgments

The present author is grateful to D. Loison, B. Delamotte, K. B. Varnashev, and H. Kawamura for useful comments.


text:
# Final state interaction contribution to the response of confined relativistic particles

## I Introduction

Scattering of high energy probes from composite systems, such as electron scattering by nuclei Frois and Sick (1991) and nucleons Ellis et al. (1996), or neutron scattering by liquid helium Silver and Sokol (1989), is often used to study the structure of the bound system. The common assumption is that in deep inelastic scattering (DIS) at sufficiently high energy the probe is incoherently scattered by the constituents of the system. In the plane wave impulse approximation (PWIA), which neglects final state interaction (FSI) effects, DIS is directly related to the momentum and energy distribution of the constituents in the target.

The role of FSI effects has been studied extensively in electron scattering from nuclear targets Benhar et al. (1991); Benhar and Pandharipande (1993) and neutron scattering from liquid helium Silver and Sokol (1989). Recently it has been suggested that they may also influence DIS of leptons by hadrons Brodsky et al. (2001). In the present study we focus on scattering from targets with confined constituents. The corresponding physical case concerns DIS from nucleons where, in distinction from the nuclear and liquid helium cases, the constituents are confined in both the initial and final states.

We Paris and Pandharipande (2001) have recently studied the response of a massless particle confined by a linear (flux-tube) potential using the simple model Hamiltonian

$$H=\sqrt{|𝐩|^2}+\sqrt{\sigma }r,$$
(1)
to gain insights into DIS of leptons by hadrons. The response to a hypothetical scalar probe is calculated as a function of the momentum and energy transfer, $`𝐪`$, $`\nu `$ in the lab frame. It is given by

$$R(𝐪,\nu )=\underset{I}{}|I|e^{i𝐪𝐫}|0|^2\delta (E_IE_0\nu ),$$
(2)
where the sum is over all energy eigenstates $`I`$. The natural scaling variable in the many-body theory approach to DIS is $`\stackrel{~}{y}=\nu |𝐪|`$ Benhar et al. (2000). In the large $`|𝐪|`$ or scaling limit the response depends only on $`\stackrel{~}{y}`$, not on $`|𝐪|`$ and $`\nu `$ independently. This scaling of the response is equivalent to the Nachtmann scaling, since the Nachtmann Nachtmann (1973) scaling variable $`\xi =\stackrel{~}{y}/M`$, where $`M`$ is the hadron mass. In Ref. Paris and Pandharipande (2001) the $`R(𝐪,\nu )`$ is calculated for $`|𝐪|10`$ GeV for the typical value of the string tension $`\sqrt{\sigma }=1`$ GeV/fm in QCD, by calculating all the relevant states $`|I`$ contributing to the response.

In this work we study the effects of the FSI of the struck particle on the response; analytic calculations of the width of the response are presented in Sec.II, and the numerical results for a linear confining potential are given in Sec.III. Both indicate that the FSI increase the width of the response beyond that predicted by PWIA. The analytic calculations also consider the nonrelativistic problem, in which $`𝐪`$ is large compared to all the momenta in the target, but smaller than the constituent mass $`m`$. The main differences between the nonrelativistic and the relativistic response are that the former peaks at $`\nu =|𝐪|^2/2m`$ and has a width proportional to $`|𝐪|`$, while the latter peaks at $`\nu |𝐪|`$, and has a constant width in the scaling limit. The folding function Benhar (1999) representing the effects of the FSI on the response is discussed in the last section, Sec.IV, where conclusions are given.

## II Moments of the Response

In the case of a single confined particle, the state of the system after the probe has struck the target is

$$|X=e^{i𝐪𝐫}|0,$$
(3)
where $`|0`$ denotes the ground state of the particle. The state $`|X`$ is not an eigenstate of the Hamiltonian and therefore has a distribution in energy. It has a unit norm, $`X|X=0|e^{i𝐪𝐫}e^{i𝐪𝐫}|0=1`$. The total strength of the response, given by the static structure function

$$S(|𝐪|)=_0^{\mathrm{}}𝑑\nu R(|𝐪|,\nu ),$$
(4)
is therefore unity. In many-body systems $`S(|𝐪|)`$ is not necessarily equal to one. Subsequent formulas pertain to the general case and show factors of $`S(|𝐪|)`$ explicitly.

The mean excitation energy of the state $`|X`$ is given by the first moment of the response:

$$\overline{\nu }(|𝐪|)=\frac{1}{S(|𝐪|)}X|HE_0|X=\frac{1}{S(|𝐪|)}_0^{\mathrm{}}𝑑\nu \nu R(|𝐪|,\nu ).$$
(5)
The width of the distribution in energy is characterized by the second moment of the energy about the mean:

$$\mathrm{\Delta }^2(|𝐪|)=\frac{1}{S(|𝐪|)}X|\left(H\frac{X|H|X}{S(|𝐪|)}\right)^2|X.$$
(6)
For example, if the response is Gaussian, it is completely determined by these three moments:

$$R_G(𝐪,\nu )=\frac{S(|𝐪|)}{\sqrt{\pi \mathrm{\Delta }^2(|𝐪|)}}e^{(\nu \overline{\nu }(|𝐪|))^2/\mathrm{\Delta }^2(|𝐪|)}.$$
(7)
The full-width at half-maximum (FWHM) of a Gaussian response is given by $`2\mathrm{\Delta }(|𝐪|)\sqrt{\mathrm{ln}2}\frac{5}{3}\mathrm{\Delta }(|𝐪|)`$. Obviously many more moments are required to describe a response of more general shape. Nevertheless, a necessary condition for the occurrence of the $`\stackrel{~}{y}`$ or Nachtmann $`\xi `$ scaling is that $`S(|𝐪|)`$, $`\overline{\nu }(|𝐪|)|𝐪|`$, $`\mathrm{\Delta }(|𝐪|)`$, and all the higher moments become independent of $`|𝐪|`$, since the response depends only on $`\nu |𝐪|`$ when $`|𝐪|\mathrm{}`$.

In the remainder of this section we will calculate the average excitation energy, $`\overline{\nu }(|𝐪|)`$ and the width, $`\mathrm{\Delta }(|𝐪|)`$ for the model Hamiltonian of Eq.(1) in the scaling limit, $`|𝐪|\mathrm{}`$. They are calculated analytically for a general potential $`V(r)`$. Section III presents numerical results for the case of a linear confining potential. This analysis allows the comparison of the results for the exact response with the following two common approximations. The PWIA, which assumes that the final state of the struck particle can be approximated by a plane wave in an average potential, is commonly used in nuclear physics. And the on-shell approximation (OSA), in which the struck particle is assumed to be in a plane wave state with the energy of that of a free particle, before and after the interaction with the probe, used in hadron physics.

We proceed with the evaluation of the average excitation energy \[Eq.(5)\] for the case of the exact response. The kinetic energy term of the matrix element $`X|H|X`$ is

$$X|\sqrt{|𝐩|^2}|X=\frac{d^3k}{(2\pi )^3}n(k)|𝐤+𝐪|,$$
(8)
where we have transformed to momentum space and $`n(k)=|\psi _0(k)|^2`$. In the scaling limit, with $`𝐪=|𝐪|\widehat{𝐳}`$, we can expand $`|𝐪+𝐤|`$ and obtain:

$$|𝐤+𝐪|=|𝐪|+k_z+\frac{1}{2|𝐪|}k_{}^2\frac{1}{2|𝐪|^2}k_{}^2k_z+𝒪\left(\frac{1}{|𝐪|^3}\right).$$
(9)
Here $`𝒪\left(\frac{1}{|𝐪|^3}\right)`$ denotes the neglected terms of that and higher order, and $`k_{}^2=k_x^2+k_y^2`$. Substituting this expansion into Eq.(5) yields,

$$\overline{\nu }(|𝐪|)=|𝐪|+V_0E_0+\frac{1}{3|𝐪|}k^2_0+𝒪\left(\frac{1}{|𝐪|^3}\right);$$
(10)
the term of $`𝒪\left(\frac{1}{|𝐪|^2}\right)`$ is zero by symmetry. Note that $`E_0V_0=T_0`$, where $`T`$ denotes the kinetic energy. Thus $`\overline{\nu }(|𝐪|)=|𝐪|T_0`$ in the limit $`|𝐪|\mathrm{}`$. The requirement that $`\overline{\nu }(|𝐪|)|𝐪|`$ becomes constant is naturally satisfied in this limit.

The width \[Eq.(6)\] depends on the matrix element $`X|H^2|X`$; $`H^2=p^2+\sqrt{p^2}V(r)+V(r)\sqrt{p^2}+V^2(r)`$. The first term of $`H^2`$ gives:

$$X|p^2|X=|𝐤+𝐪|^2_0=|𝐪|^2+k^2_0.$$
(11)
The average of the second and third terms in $`H^2`$ may be shown to be equal to each other and

$`X|\sqrt{p^2}V|X`$ $`=`$ $`{\displaystyle \frac{d^3k}{(2\pi )^3}\psi _0^{}(k)|𝐪+𝐤|[V\psi _0](k)}`$ (12)
$``$ $`|𝐪|V_0+{\displaystyle \frac{1}{2|𝐪|}}k_{}^2V_0+𝒪\left({\displaystyle \frac{1}{|𝐪|^3}}\right),`$
where $`[V\psi _0](k)=d^3re^{i𝐤𝐫}V(r)\psi _0(r)`$. Substitution into Eq.(6) with $`S(|𝐪|)=1`$ gives:

$$\mathrm{\Delta }^2(|𝐪|)=\frac{1}{3}k^2_0+V^2_0V_0^2+\frac{2}{3|𝐪|}\left(k^2V_0k^2_0V_0\right)+𝒪\left(\frac{1}{|𝐪|^2}\right).$$
(13)
This expression demonstrates that the width of the exact response is independent of $`|𝐪|`$ in the limit $`|𝐪|\mathrm{}`$ as necessary for $`\stackrel{~}{y}`$ scaling. It also shows that the width has a kinematic contribution dependent upon the target momentum distribution, and an additional interaction contribution.

As mentioned, the PWIA assumes that a constituent of momentum $`𝐤`$, after being struck by the probe, may be described by a plane wave with momentum $`𝐤+𝐪`$ in an assumed average potential chosen to give the exact $`\overline{\nu }`$ of Eq.(10). The PWIA response is

$$R_{PWIA}(|𝐪|,\nu )=\frac{d^3k}{(2\pi )^3}n(k)\delta (e(𝐤+𝐪)E_0\nu )$$
(14)
where $`e(𝐤+𝐪)`$ is the energy of the plane wave, taken to be:

$$e(𝐤+𝐪)=|𝐤+𝐪|+V_0.$$
(15)
The first term is the kinetic energy of the struck particle. The second term is the average of the potential in the final state of the system given by, $`X|V|X=0|e^{i𝐪𝐫}Ve^{i𝐪𝐫}|0=V_0`$. The $`R_{PWIA}(|𝐪|,\nu )`$ exhibits $`\stackrel{~}{y}`$ scaling as can be easily seen by expanding the argument of the $`\delta `$ function in $`R_{PWIA}`$. In the large $`|𝐪|`$ limit $`R_{PWIA}`$ depends only on $`\stackrel{~}{y}=\nu |𝐪|`$.

The average energy and the width of the PWIA response are calculated using:

$$_0^{\mathrm{}}𝑑\nu \nu ^nR_{PWIA}(|𝐪|,\nu )=_0^{\mathrm{}}𝑑\nu \nu ^n\frac{d^3k}{(2\pi )^3}n(k)\delta (e(𝐤+𝐪)E_0\nu ).$$
(16)
The average excitation energy, obtained with $`n=1`$, agrees with the exact result in Eq.(10) by construction.

The width of the PWIA response, however, is not the same as the width of the exact response.

$$\mathrm{\Delta }_{PWIA}^2=\frac{1}{3}k^2_0+𝒪\left(\frac{1}{|𝐪|^2}\right)$$
(17)
contains only the first term of the exact result \[Eq.(13)\] due to the target momentum distribution. The second term, $`V^2_0V_0^2`$, of the exact $`\mathrm{\Delta }^2`$ represents the FSI contribution neglected in the PWIA. It does not vanish in the $`|𝐪|\mathrm{}`$ limit for relativistic kinematics.

In the non-relativistic case, $`H_{NR}=\frac{|𝐩|^2}{2m}+V(r)`$, the exact $`\overline{\nu }`$ is given by:

$$\overline{\nu }_{NR}=\frac{1}{2m}|𝐪|^2.$$
(18)
The nonrelativistic-PWIA (NR-PWIA) response:

$$R_{NRPWIA}(|𝐪|,\nu )=\frac{d^3k}{(2\pi )^3}n(k)\delta \left(\frac{|𝐪+𝐤|^2}{2m}+V_0E_0\nu \right)$$
(19)
gives the exact $`\overline{\nu }`$. For the width of the NR-PWIA response we obtain:

$$\mathrm{\Delta }_{NRPWIA}^2(|𝐪|)=|𝐪|^2\frac{k^2_0}{3m^2}+\frac{1}{4m^2}\left(k^4_0k^2_0^2\right).$$
(20)
Note that in Eqs.(18) and (20) we have not taken the $`|𝐪|\mathrm{}`$ limit.

The width of the exact NR response is:

$$\mathrm{\Delta }_{NR}^2(|𝐪|)=\mathrm{\Delta }_{NRPWIA}^2+\frac{1}{m}\left(Vk^2_0V_0k^2_0\right)+V^2_0V_0^2.$$
(21)
It differs from $`\mathrm{\Delta }_{NRPWIA}`$ in terms of order $`1/|𝐪|`$ which can be neglected in the scaling limit. Thus, in contrast to the relativistic case, the FSI do not increase the width of the NR-PWIA response at large $`|𝐪|`$.

Finally we consider the on-shell approximation (OSA) in which the energy of the struck constituent is that of a free relativistic particle before and after the interaction with probe, as assumed in the quark-parton model. The response in OSA is

$$R_{OSA}(|𝐪|,\nu )=\frac{d^3k}{(2\pi )^3}n(k)\delta \left(|𝐤+𝐪||𝐤|\nu \right);$$
(22)
it depends only on the momentum distribution of target constituents and obeys $`\stackrel{~}{y}`$ scaling. The average excitation in OSA is

$$\overline{\nu }_{OSA}=|𝐪|T_0+\frac{1}{3|𝐪|}k^2_0+𝒪\left(\frac{1}{|𝐪|^3}\right),$$
(23)
and the width is given by:

$$\mathrm{\Delta }_{OSA}^2=\frac{1}{3}k^2_0+k^2_0T_0^2+𝒪\left(\frac{1}{|𝐪|}\right).$$
(24)
The exact value of $`\overline{\nu }`$ \[Eq.(10)\] is reproduced by the OSA for any potential. However, the $`\mathrm{\Delta }_{OSA}^2`$ has $`k^2_0T_0^2`$ in place of the $`V^2_0V_0^2`$ in the leading term of the exact $`\mathrm{\Delta }^2`$ \[Eq.(13)\]. For a massless particle in a linear confining potential, i.e. for the Hamiltonian of Eq.(1), $`T_0=V_0`$, and $`k^2_0=V^2_0`$. Therefore for this particular Hamiltonian the OSA reproduces the exact value of $`\mathrm{\Delta }`$.

The calculated responses shown in the next section however indicate that the shape of the OSA response is incorrect. We have calculated the third moment of the $`H`$ in the state $`|X`$ and observed that the exact and OSA responses indeed give different values. We expect the moments to differ also for orders higher than the third since the shapes of the exact and the OSA response are very different.

## III Numerical results

We first compare the response functions for $`|𝐪|=10`$ GeV before comparing their moments. In Ref. Paris and Pandharipande (2001) it has been shown that the scaling limit is obtained for such values of $`|𝐪|`$. The exact response, Eq.(2), is a sequence of $`\delta `$ functions at $`\nu =E_IE_0`$. In order to obtain a smooth response we assume decay widths $`\mathrm{\Gamma }_0`$ for all the excited states. Note that the energies of the states $`|I`$ that contribute to the response at $`|𝐪|=10`$ GeV are large, therefore their decay widths are not affected by the energy dependent terms assumed in Ref. Paris and Pandharipande (2001). The response including decay widths is given by:

$$R(𝐪,\nu )=\underset{I}{}|I|e^{i𝐪𝐫}|0|^2\left(\frac{\mathrm{\Gamma }_0}{2\pi }\right)\frac{1}{(E_IE_0\nu )^2+\mathrm{\Gamma }_0^2/4}.$$
(25)
The responses obtained with $`\mathrm{\Gamma }_0=100`$ and 50 MeV are shown in Fig. 1, along with the PWIA and OSA responses for $`|𝐪|=10`$ GeV and for $`|𝐪|\mathrm{}`$. The difference between the exact responses for $`\mathrm{\Gamma }_0=100`$ and 50 MeV are much smaller than those between the exact and the approximate.

We note that the shape of the PWIA response is qualitatively similar to that of the exact, however, its width is too small. This is a direct consequence of the neglect of interaction terms in $`\mathrm{\Delta }`$ \[Eq.(13)\] as discussed in the last section. The width $`\mathrm{\Delta }`$ of the response for $`\mathrm{\Gamma }_00`$, is 409 MeV, while the $`\mathrm{\Delta }_{PWIA}=326`$ MeV. Note that for a Gaussian response, Eq.(7) the FWHM is $`5\mathrm{\Delta }/3`$. The FWHM of the exact and the PWIA responses shown in Fig. 1 are larger because they are not exactly Gaussian.

The OSA results in the discontinuous curves shown in Fig. 1. They are discontinuous at the lightline ($`|𝐪|=\nu `$) because the response of free particles is limited to the spacelike region $`\nu <|𝐪|`$. The discontinuity at $`\stackrel{~}{y}=0`$ is in clear conflict with the exact response which is continuous across the lightline and is non-zero in the timelike ($`\stackrel{~}{y}>0`$) region. Therefore the OSA appears to be unsatisfactory even though for the special case of a linear potential it has the exact values of $`S(|𝐪|)`$, $`\overline{\nu }(|𝐪|)`$ and $`\mathrm{\Delta }(|𝐪|)`$.

## IV FSI Folding function

The exact response may be obtained from $`R_{PWIA}(|𝐪|,\nu )`$ via convolution with the folding function Silver and Sokol (1989); Benhar et al. (1991); Benhar (1999). The states $`|𝐤+𝐪`$ used in the PWIA are not eigenstates of the $`H`$. Therefore they must have a distribution in energy, and we may express the exact response as:

$$R(|𝐪|,\nu )=𝑑\nu ^{}F(|𝐪|,\nu \nu ^{})R_{PWIA}(|𝐪|,\nu ^{}).$$
(26)
The folding function $`F`$ defined above, is primarily meant to describe the width of the plane wave states. Scaling occurs even in presence of interactions because this folding function becomes independent of $`|𝐪|`$ at large $`|𝐪|`$ Benhar (1999); Weinstein and Negele (1982). Occurrence of scaling does not imply that either the PWIA or the OSA is valid. Fig. 2 shows the folding function obtained from the PWIA and exact responses for $`|𝐪|=10`$ GeV and $`\mathrm{\Gamma }_0=100`$ and 50 MeV. They are calculated numerically from $`R(|𝐪|,\nu )`$ and $`R_{PWIA}(|𝐪|,\nu )`$ by deconvolution.

In order to extract information on the ground state wave function from the measured response, we must first remove the FSI effects represented by the folding function $`F(|𝐪|,\nu )`$. However, in most cases the folding function is not known. In some nonrelativistic cases it can be calculated with quantum Monte Carlo methods Carraro and Koonin (1991). It can also be estimated using the Glauber approximation Silver and Sokol (1989); Benhar (1999). Here we consider the possibility of approximating the $`F`$ by a Gaussian folded with the Lorentzian distribution used in Eq.(25). In this case the width of the Gaussian is given by the interaction contribution to the exact $`\mathrm{\Delta }`$ Eq.(13):

$`\mathrm{\Delta }_I^2`$ $`=`$ $`V^2_0V_0^2`$ (27)
$`F_G(\nu \nu ^{})`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\pi \mathrm{\Delta }_I^2}}}e^{(\nu \nu ^{})^2/\mathrm{\Delta }_I^2},`$ (28)
and the convolution of $`F_G`$ with the Lorentzian is given by:

$`F_A(\nu \nu ^{})`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\nu ^{\prime \prime }F_G(\nu \nu ^{\prime \prime })F_L(\nu ^{\prime \prime }\nu ^{}),`$ (29)
$`F_L(\nu ^{\prime \prime }\nu ^{})`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Gamma }_0}{2\pi }}\right){\displaystyle \frac{1}{(\nu ^{\prime \prime }\nu ^{})^2+\mathrm{\Gamma }_0^2/4}}.`$ (30)
The $`F_A`$ is shown in Fig. 2. It does not have the finer structures in the exact $`F`$, and it peaks at $`\nu \nu ^{}=0`$, while the $`F`$ peaks at $`0.2`$ GeV. Nevertheless, a much better approximation to the response is obtained by folding the PWIA response with the $`F_A`$ as shown in Fig. 3.

In conclusion, we find that the effects of FSI, neglected in the PWIA, become independent of $`|𝐪|`$ at large $`|𝐪|`$ for a particle confined in a linear potential. They increase the width of the response, and must be removed before extracting momentum distributions or other target structure information from the observed response. In general one can always represent the FSI effects by a folding function; in this simple case a Gaussian folding function provides a reasonable approximation. Finally, the OSA predicts a response of the wrong shape but with correct values for the first three moments in the case of a linear potential.

###### Acknowledgements.

The authors thank Omar Benhar and Ingo Sick for many discussions. This work has been partly supported by the US National Science Foundation via grant PHY 00-98353.


text:
# Magnetic properties of 3d-impurities substituted in GaAs

## Abstract

We have calculated the magnetic properties of substituted 3d-impurities (Cr- Ni) in a GaAs host by means of first principles electronic structure calculations. We provide a novel model explaining the ferromagnetic long rang order of III-V dilute magnetic semiconductors. The origin of the ferromagnetism is shown to be due to delocalized spin-uncompensated As dangling bond electrons. Besides the quantitative prediction of the magnetic moments, our model provides an understanding of the halfmetallicity, and the raise of the critical temperature with the impurity concentration.

As the electronic device sizes are shrinking, physicists are asked to provide new device concepts. Devices that work, not despite quantum mechanics but because of quantum mechanics. One track is the area of magnetoelectronics, where one in addition to the charge of the electron tries to exploit its spin degree of freedom . A central problem is here the spin injection into the semiconductor . The idea of using a ferromagnetic semiconductor, i.e. a diluted magnetic semiconductor (DMS), for the spin injection is now being worked on worldwide . So far, the mechanism for the ferromagnetic long range order in DMSs has been a matter of dispute. Akai  and more recently Dietl et al.  proposed the double-exchange mechanism to explain the ferromagnetic order.

In this letter, we propose a novel mechanism for the ferromagnetism of DMSs. We show that the ferromagnetic order is caused by itinerant spin-uncompensated As dangling bond (DB) states: The 3d-impurity atom interacts with the As DBs. Because of the spin-polarization of the 3d-atom, the hybridization with the spin-up As DBs is larger than with the spin-down ones. This has two consequences: i) Since the DB states are delocalized, the hybridization leads to an itinerant magnetic moment. ii) The filling of the As DBs remains incomplete. An eventual filling of these holes destroys the ferromagnetic order in agreement with experiment .

To support our analysis we perform ab initio calculations of Ga-substituted impurities (Cr, Mn, Fe, Co, and Ni) in a GaAs host. The magnetic properties were calculated using a plane wave pseudopotential code (VASP ) within the local spin density approximation (LSDA) . The atoms are described by ultra soft Vanderbilt type pseudopotentials  as supplied by G. Kresse and J. Hafner . The wave functions are expanded in plane waves with an energy cut-off of 295 eV. The electron density was calculated using special k-point sets  corresponding to a 2 x 2 x 2 folding. The resulting number of k-points was shown to be sufficient for the here studied properties. The calculations are performed within a 64 atom supercell, where the 3d-atom is put on a Ga-site in the center of the cell.

In Fig.1 we show the spin-resolved d-projected density of states (DOS) at the impurity site. For the unrelaxed case (left column), i.e. all atoms are sitting on ideal high symmetry positions, GaAs containing Cr, Mn, or Co (Fe or Ni) impurities has become a half-metal (metal). In the right column of Fig.1 we show the results for the position, volume, and shape relaxed case. For Cr or Mn impurities in GaAs, the effect of the relaxation is negligible, whereas for Fe, Co, and Ni the relaxation almost suppresses the ferromagnetism. Co or Ni impurities in GaAs are non-magnetic, for Fe in GaAs the magnetic moment has drastically decreased (Fig.2).

In order to understand these ab initio results, we apply an ”impurity-molecule model”  which we develop in three steps. i) We consider the ideal, i.e unrelaxed, GaAs host with a Ga vacancy. At the vacancy site the four surrounding As neighbours contribute five electrons to the four DBs. We calculate the Ga-vacancy to be non-magnetic. The site-projected vacancy DOS is shown in Fig.3. Close to the Fermi energy ($`E_F`$) a split peak appears in the DOS of the Ga-vacancy (bottom panel of Fig.3). The two peaks have the following origin: Three As DBs give rise to the peak at $`E_F`$ (three top panels), whereas one As DB gives rise to the peak below $`E_F`$ (panel four of Fig.3). We conclude, that the occupation of the As DBs is as follows: One DB is completely occupied (peak below $`E_F`$) and three DBs are only half-occupied (peak at $`E_F`$).

In total the vacancy has 3 holes, i.e. there are 1.5 holes per spin. In a tetrahedral crystal field the $`sp^3`$ orbitals split into ”s ”-like $`a_1`$ orbitals and ”p ”-like $`t_2`$ orbitals . At the vacancy site the $`a_1`$ state is completely filled and the $`t_2`$ states have 3 holes. The $`t_2`$ state of the Ga vacancy is located $`0.06eV`$ above the valence band edge and is thus almost degenerate with the continuum of extended states . The Ga vacancy $`t_2`$ state can thus be treated as a delocalized state. Only then the nonmagnetic state of the vacancy is understandable.

ii) Next the vacancy is filled with the 3d-impurity. In a tetrahedral crystal field the d-states are split into $`e`$ and $`t_2`$ states, where the $`e`$ states lie lower in energy than the $`t_2`$ states . For all 3d impurities (with the exception of Ni) only the spin-up $`t_2`$ ($`t_2^{}`$) states are occupied, i.e. the $`t_2^{}`$ states are completely empty. We assume that the impurity retains its atomic character, i.e. the occupation of the e and $`t_2`$ levels follows Hund’s rules.

iii) Finally, we consider the interaction of the 3d-impurity with the four As DBs. It is well known that a substitutional impurity in a semiconductor is ionized so that the semiconductor holes become filled. One would thus expect the 3d-atom to contribute 3 electrons to the hole-filling, i.e a 3+-state of the 3d impurity. However, for a spin-polarized impurity there are two possibilities: (1) If the exchange energy is larger than the energy gain (bonding energy) due to the complete hole filling of the As DBs, only the spin-up channel will interact with the As DBs, i.e a spin-flip on the 3d impurity site is prohibited. (2) If, in contrast to the first possibility, the bonding energy exceeds the exchange energy, 3 electrons will be transferred to the As DBs, i.e a spin-flip on the 3d impurity site is allowed. In the following, we will discuss the two cases separately.

First assume case (1), i.e that the exchange energy exceeds the bonding energy. This case agrees with the results for the unrelaxed first-principles calculation. Due to their like symmetry and large overlap, only the $`t_2^{}`$ states of the impurity and the As DB hybridize. Therefore the 3d-impurity contributes only 1.5 spin-up electrons to the hole-filling of the As DBs. Because the 3d impurity states are localized, not 1.5 but an integer number of two electrons are promoted to the Fermi energy, of which 0.5 electrons have a high probability to be found on the impurity site itself. Thus one electron is shared between the impurity and the delocalized As DBs. The occupation of the impurity atom $`t_2^{}`$ level is therefore reduced by two electrons (states), i.e. the 3d atom assumes a 2+ state. In the following we discuss this charge transfer in more detail:

The As DBs receive 1.5 $`t_2^{}`$ electrons. In order to minimize the energy (especially the kinetic energy of the host), it seems at first sight, that the optimal situation would be to equally distribute the 1.5 electrons on both spin-channels of the host, i.e. 0.75 electrons had to spin flip. But, this apparently optimal situation is not allowed, because of the following: On the impurity site only spin-up states are available at the Fermi energy ( due to the exchange splitting of the d-states). Since, as mentioned above, one promoted electron is shared between the impurity atom and the delocalized As DBs, the spin of one of the two promoted electrons is restricted to be spin-up. In order to minimize the energy, only one electron of the 1.5 transferred $`t_2^{}`$ electrons, is allowed to be equally distributed among both spin-channels, i.e. only 0.5 electrons are allowed to spin flip. In total then, the spin-up (spin-down) channel of the host has received 1.0 (0.5) electron.

An alternative discussion of the charge transfer is the following: Let us start from the vacancy dangling bonds. As mentioned earlier, three orbitals of the four delocalized As DBs, are half occupied (0.5 electron per spin) and one orbital is completely occupied (1 electron per spin) , because five electrons have to be redistributed on four DBs). Assume now that we add one electron to the vacancy and distribute it equally among both spin channels. Adding 0.5 electrons to each spin channel of the As DBs, gives rise to the following occupation of the As DBs: three orbitals are completely occupied ( 1 electron per spin) and one orbital is completely empty ( 0 electron per spin), because six (five plus one) electrons have to be redistributed on four DBs. This is an energetically very favorable situation, because the DOS at the Fermi energy is drastically reduced to zero, i.e. a bandgap exists. The 3d-impurity donates to the As DBs 1.5 $`t_2^{}`$ electrons, of which one is used to create the bandgap as described above. Now, because not only one electron, but 1.5 $`t_2^{}`$ electrons are transferred form the impurity to the As DBs, the spin-up channel of the As DBs receives an additional 0.5 electron. Therefore the spin-down channel has a bandgap, whereas the spin-up channel has one partly occupied orbital, i.e a state at the Fermi energy (Fig.1). This additional 0.5 electron corresponds to one half of the shared electron. This discussion also shows, that any other distribution of the 1.5 $`t_2^{}`$ electrons within the As DBs, will increase the total energy, because the bandgap in the spin-down channel would have to disappear.

In total the impurity has lost 1.5 spin-up electrons and the As DBs have gained 0.5 spin-down electrons and 1.0 spin-up electrons. Hence the occupation of the spin-up As DBs exceeds the occupation of the spin-down As DBs by 0.5 electrons and this spin-imbalance is accompanied by an itinerant magnetic moment, $`M_{itinerant}=0.5\mu _B`$. Note, that although this amounts to a relatively small moment per formula unit, it constitutes the origin of the long range order, i.e the ferromagnetism of the DMS. Our model also gives a natural explanation for the observation of so called ”hole mediated ferromagnetism” . Namely, any spin-imbalance (ferromagnetism) will be quenched, if the remaining 1.5 holes of the As DBs are filled (disappear) for instance by additional donor-dopants.

We can now express the (local) magnetic moment at the impurity site simply as follows:

$$M_{local}=d^{}d^{}N,$$
(1)
where N is equal to 1.5 (unrelaxed case), and $`d^{},d^{}`$ correspond to the d-occupation of the impurity atom.

The total magnetic moment per impurity atom, $`M_{total}`$, is the sum of the local and itinerant moment:

$$M_{total}=M_{local}+M_{itinerant}=d^{}d^{}1.$$
(2)
In the case of, for example, (Ga,Mn)As, the local Mn magnetic moment amounts to $`M_{local}=501.5=3.5\mu _B`$ and $`M_{total}`$ amounts to $`4\mu _B`$. In Fig.2 we summarize our results for Cr, Mn, Fe, Co, and Ni impurities in a GaAs host. We compare our model (Eq.1-2) with the self consistently calculated magnetic moments. For the unrelaxed case (left column) discussed so far, one notices an overall agreement with some deviations, which however can be explained easily. For Cr and Mn the total magnetic moment agrees perfectly, whereas the itinerant magnetic moment is smaller by about $`0.2\mu _B`$ and the local magnetic moment is larger by the same amount than predicted by our model. In our model, we assumed one electron to be shared equally between As and the 3d-impurity. Realizing that due to the difference in electronegativity, more than half of the shared electron is located on the impurity site, the deviations are obvious and follow the expected chemical trend.

The other deviation we find for Ni and Fe. This can also be easily understood from the calculated d\- and site- projected densities of states (Fig.1). For Fe and Ni we do find spin-down states at the Fermi energy: in case of Fe these are $`e^{}`$ states and in case of Ni $`t_2^{}`$ states. These additional states contribute to the minimization of the total energy and are for simplicity neglected in our model, though it is straight forward to expand the model.

Now we compare our model with the calculated DOS (Fig.1). In agreement with our model, we find two important features in the spin-up channel: (1) one deep lying broad feature originating from the localized $`t_2`$ and $`e`$-states on the impurity site ( $`E1.7eV`$). For Cr this broad feature does not exist, because the Cr atom only has 2 $`t_2`$ electrons. Thus for Cr there is an additional empty $`t_2`$ peak above the Fermi energy. (2) one double $`t_2`$-peak close to the Fermi energy originating from the shared electron between the Ga-vacancy and the impurity.

So far we have not considered any relaxations that increase the overlap between the As DBs and the impurity $`t_2`$ states. With an increased overlap, the gain in bonding energy, $`E_{bo}`$, might become larger than the energy gain due to the exchange splitting, $`E_x`$. Then 3 electrons of the 3d-impurity hybridize with the As DBs, because the cost to spin-flip 1.5 d-electrons (of the impurity) is compensated by the gain in bonding energy. Therefore the impurity will undergo a transition from a high-spin state to a low-spin state. A similar discussion for 3d transition metal ions in Si has earlier been given by Beeler et al. . For Cr and Mn the exchange energy is large and there are no spin-down states at the Fermi level into which a $`t_2`$ spin up electron could flip. The system therefore would loose energy when increasing the overlap between As and Cr (Mn). That is why the As-Cr and As-Mn bond length is more or less maintained.

For Fe, Co and Ni on the other hand, the system gains energy, when 1.5 $`t_2`$ electrons spin-flip to fill all As DB holes, and consequently the bond length is thereby decreased by about 7 %. (Co,Ga)As and (Ni,Ga)As are non-magnetic , whereas $`M_{total}`$ amounts to $`1\mu _B`$ for (Fe,Ga)As. Our model (N=3 in Eq.1) predicts for Fe a local magnetic moment of $`1\mu _B`$, but a zero itinerant magnetic moment. With our first principles calculations we find an itinerant moment of $`0.2\mu _B`$ and a local moment of $`0.8\mu _B`$. This difference is explained as follows: Due to the difference in electronegativity between Fe and As, some of the delocalized electrons become localized on the Fe site. Since the $`e^{}`$ states at the Fermi energy are completely filled, only spin-down electrons can get localized at the Fe site. It follows, that the local magnetic moment is reduced by the amount, $`A`$, of localized spin-down electrons. From Fig.2 one would expect an $`A`$ between -0.1 and -0.2 $`\mu _B`$. This in turn leaves an itinerant magnetic moment of $`A\mu _B`$, in agreement with our first principles calculations. We thus find the ferromagnetism of (Fe,Ga)As to be independent on the number of holes in the DMS.

In order to understand the increase of the critical temperature, $`T_c`$, with the increase of the Mn impurity concentration , consider the following: The itinerant magnetic moment per host atom increases with Mn concentration, whereas the total magnetic moment per impurity atom is independent on Mn concentration. Because in our model only the itinerant magnetic moment is responsible for the long range order, we assume $`T_c`$ to be proportional to the itinerant magnetic moment. Accordingly, $`T_c`$ increases with Mn concentration.

Recent theoretical calculations on (Ga,Mn)As report the As nearest neighbour (nn) local magnetic moment to be antiparallel to Mn . But this does not contradict with our model as we will illustrate in the following: First, notice that the spin down channel of the As DBs is more localized than the spin up channel, because of the following reason: The 0.5 additional electrons in the spin-up channel shifts the spin-up DB states deeper into the valence band, and hence increases its resonance with the valence band, and in effect the spin up DB states become more delocalized.

Second, it thus follows that on the As nn sites, the local As moment might be determined by the spin down channel and accordingly the local As moment may be antiparallel to the impurity moment. The sign of the local As moment depends thus on the chosen local sphere radius. The smaller the chosen sphere radius, the more negative becomes (also in our calculations) the local As moment. However, most important, these recent calculations confirm the existence of an integer total magnetic moment and a non-integer local moment at the impurity site.

Experimentally, often a magnetic moment of $`5.0\mu _B`$ is measured for (Ga,Mn)As. We did a calculation where we add one electron to the unit cell which becomes charge compensated by a jellium background. In this way we simmulate the electronic effect of additional dopants in the unit cell, whose defect-level(s) determine the Fermi level to lie above the pure (Ga,Mn)As Fermi energy. We find a total moment of $`5.0\mu _B`$ and a local moment of $`4.0\mu _B`$. Within our model this is understandable in the following way: Adding one electron, the vaccancy has instead of 3 only 2 holes, i.e. one hole per spin. The energetically most stable situation is with 3 fully occupied orbitals and one empty orbital. Out of the 5 spin-up d-electrons now only one electron is transferred to the As DBs. This electron will maintain the spin, because the system would loose energy, when both (spin-up and spin-down) empty orbitals would become filled by half an electron. Moreover, this transferred spin-up electron is delocalized, i.e the itinerant moment is $`1.0\mu _B`$. Therefore, the Curie temperature of the negatively charged Mn impurity should be higher than of the neutral Mn impurity.

In summary, the present model explains the ferromagnetic long range order as being caused by itinerant spin uncompensated $`t_2`$ orbitals. The simultaneous occurrence of an itinerant and local magnetic moment is explained. Our model quite naturally explains the halfmetallicity of (Ga,Mn)As. Moreover, because the itinerant moment per host atom increases with the impurity concentration, our model explains the dependence of $`T_c`$ on the impurity concentration.

###### Acknowledgements.

We are grateful to the Swedish Foundation for Strategic Research (SSF) , the Swedish Natural Science Research Council (NFR), and the Göran Gustafsson Foundation for financial support. We thank J. Hafner and G. Kresse for giving us the possibility to use VASP.