Integral equation and simulation studies of the Heisenberg spin fluid in an external magnetic field F. Lado* and E. Lomba Instituto de Quı ´mica Fı ´sica Rocasolano, CSIC, Serrano 119, 28006 Madrid, Spain J. J. Weis Laboratoire de Physique The ´orique et Hautes Energies, Ba ˆtiment 211, Universite ´ de Paris-Sud, 91405 Orsay Cedex, France Received 26 February 1998 We develop a general method to study inhomogeneous liquids in an external field using orthogonal poly- nomials tailored to the one-body density. The procedure makes integral equation calculations of these systems no more difficult than those of ordinary homogeneous molecular fluids. We apply this method to the ferro- magnetic Heisenberg spin fluid in an external magnetic field using both the reference-hypernetted chain closure and a reference version of the Zerah-Hansen closure, with no further approximations. The calculation includes a mapping of the two-phase region for several values of the external magnetic field. Comparison with Monte Carlo simulation data shows the integral equation procedure yielding nearly exact results, in particular for nonzero external fields. S1063-651X9815109-7 PACS numbers: 61.20.Gy, 64.60.-i, 75.10.-b, 75.30.-m I. INTRODUCTION The Gibbsian N -body density function of a Hamiltonian H N that is rotationally and translationally invariant must it- self be rotationally and translationally invariant, as must then also be all reduced n -body density functions of this Hamil- tonian. In particular, the one-body density is a constant and the two-body density depends only on relative coordinates. An external field destroys this homogeneity, producing an- isotropy or nonuniformity 1,2, and so makes necessary the joint calculation of the coupled one-body and two-body den- sity functions. A striking if familiar example of the response of a bulk system to an external field is ferromagnetism. In this paper, we shall use this particular case to present a gen- eral procedure 3to compute the coupled one-body and two- body density functions of an inhomogeneous classical fluid in an external field. Remarkably, the procedure is no more difficult to carry through than similar calculations for ordi- nary homogeneous systems. Perhaps the simplest model of a disordered continuum system exhibiting ferromagnetic behavior is a fluid of hard spheres with embedded Heisenberg spins described using classical statistical mechanics 4–8. On the one hand, the simplifications of this model compared to magnetic dipole- dipole interactions are significant for simulation. The dipole- dipole model presents substantial conceptual difficulties in connection with the nonexistence of the thermodynamic limit for orientationally ordered phases. Moreover, its simulation results are strongly dependent on the boundary conditions. Simulations for the Heisenberg spin model are free of these problems. The integral equation formalism developed below serves equally well for either potential.On the other hand, disordered Heisenberg spin systems are of interest in their own right. It has been known for some time that Co/P alloys 9and Co/Au melts 10have a tendency to form amor- phous ferromagnets, although their ‘‘true liquid’’ nature has been questioned because of the technical difficulties posed by undercooling below the Curie temperature. While these systems might be more representative of quenched spin flu- ids, very recently Albrecht and co-workers 11have finally managed to undercool a Co 80 Pd 20 melt below its Curie tem- perature at zero field, thus obtaining the first evidence of ferromagnetic behavior in a liquid metal in conditions where the Heisenberg exchange interaction absolutely dominates the magnetic dipole-dipole contribution. The exchange inter- action is the crucial term in the present paper and it is our main aim to analyze how it models the phase behavior with and without the presence of external fields. The Heisenberg spin fluid in an external magnetic field B 0 is defined by the canonical partition function Z = 1 N ! 3 N j =1 N d r j d j exp j j B 0 - i j u 0 r ij - i j u ss r ij , i , j . 1 This can be factored into an ideal part and the excess, Z =Z id Z ex , where Z id = 1 N ! 3 N j =1 N d r j d j exp j j B 0 = 1 N ! V 3 N 4 sinh  B 0  B 0 N , 2 *On leave from Department of Physics, North Carolina State Uni- versity, Raleigh, NC 27695-8202. PHYSICAL REVIEW E SEPTEMBER 1998 VOLUME 58, NUMBER 3 PRE 58 1063-651X/98/583/347812/$15.00 3478 © 1998 The American Physical Society