text:
# Relativistic Jet Response to Precession and Wave-Wave Interactions

## 1 Introduction

Relativistic jets, particularly in extragalactic sources, can exhibit time-dependent curved structures with superluminally moving components, e.g., 3C 345 (Zensus, Cohen & Unwin 1995) or with both superluminally moving and much more slowly moving or stationary components, e.g., M 87 (Biretta, Zhou & Owen 1995; Biretta, Sparks & Macchetto 1999), 4C 39.25 (Alberdi et al. 2000), 3C 120 (Walker et al. 2001). Superluminal motions along curved trajectories can be explained by helical jet models (Hardee 1987; Steffan et al. 1995), and helical patterns are the expected result for Kelvin-Helmholtz and current driven jet instabilities in relativistic flows (Birkinshaw 1991, and references therein; Appl 1996; Istomin & Pariev 1996). It has been proposed that a combination of moving and stationary components can be the result of enhanced features moving with the jet flow through relatively fixed curved helical structures with Doppler boosting leading to fixed components where the average flow comes most nearly towards the line-of-sight (Alberdi et al. 2000; Walker et al. 2001).

Axisymmetric relativistic jet simulations have investigated component motion by introducing velocity perturbations at the origin and studying the subsequent evolution. With repeated velocity perturbations (Gomez et al. 1997; Midouszewski, Hughes & Duncan 1997) the resulting knotty structure moves with the jet flow. In a more recent axisymmetric simulation performed by Agudo et al. (2001), a single velocity perturbation generates a superluminal component associated with the velocity perturbation and multiple trailing sub-luminally moving accelerating components. These simulations address some of the issues involving components moving at different speeds but do not address the issue of components moving through stationary features or curved trajectories. To this end several 3-D hydrodynamic simulations of relativistic jets have been performed (Aloy et al. 1999; Aloy et al. 2000; Hughes et al. 1999, 2001) in which the effect of perturbation induced asymmetries on jet propagation were investigated. However, these studies did not address the issue of internal jet structure associated with asymmetric fluid motion and the component motions that might be observed inside the jets.

In this paper we present fully 3-D hydrodynamic simulations along with a detailed analysis of the time-dependent structures that develop in the jet for three different jet precession frequencies but with no velocity variation other than the relatively small transverse velocity induced by the jet precession. The numerical techniques are described in §2. Results and analysis of jet structures seen in the simulations are presented in §3 and §4. The potential observable features resulting from these structures are shown in §5, and in §6 we summarize and discuss some of the implications of our results.

## 2 CFD Simulations

### 2.1 Solution of the Euler Equations

We assume an inviscid and compressible gas, and an ideal equation of state with constant adiabatic index. We use a Godunov-type solver which is a relativistic generalization of the method of Harten, Lax, & Van Leer (1983), and Einfeldt (1988), in which the full solution to the Riemann problem is approximated by two waves separated by a piecewise constant state. We evolve the mass density $`R`$, the three components of the momentum density $`M_x`$, $`M_y`$ and $`M_z`$, and the total energy density $`E`$ relative to the laboratory frame.

Defining the vector

$$U=(R,M_x,M_y,M_z,E)^T,$$
(1)
and the three flux vectors

$$F^x=(Rv^x,M_xv^x+p,M_yv^y,M_zv^z,(E+p)v^x)^T,$$
(2)
$$F^y=(Rv^y,M_xv^x,M_yv^y+p,M_zv^z,(E+p)v^y)^T,$$
(3)
$$F^z=(Rv^z,M_xv^x,M_yv^y,M_zv^z+p,(E+p)v^z)^T,$$
(4)
the conservative form of the relativistic Euler equation is

$$\frac{U}{t}+\frac{}{x}(F^x)+\frac{}{y}(F^y)+\frac{}{z}(F^z)=0.$$
(5)
The pressure is given by the ideal gas equation of state $`p=(\mathrm{\Gamma }1)(en)`$ where here and in what immediately follows we have set $`c=1`$. The Godunov-type solvers are well known for their capability as robust, conservative flow solvers with excellent shock capturing features. In this family of solvers one reduces the problem of updating the components of the vector $`U`$, averaged over a cell, to the computation of fluxes at the cell interfaces. In one spatial dimension the part of the update due to advection of the vector $`U`$ may be written as

$$U_{}^{n+1}{}_{i}{}^{}=U_{}^{n}{}_{i}{}^{}\frac{\delta t}{\delta x}(F_{i+\frac{1}{2}}F_{i\frac{1}{2}}).$$
(6)
In the scheme originally devised by Godunov (1959), a fundamental emphasis is placed on the strategy of decomposing the problem into many local Riemann problems, one for each pair of values of $`U_i`$ and $`U_{i+1}`$ to yield values which allow the computation of the local interface fluxes $`F_{i+\frac{1}{2}}`$. In general, an initial discontinuity at $`i+\frac{1}{2}`$ due to $`U_i`$ and $`U_{i+1}`$ will evolve into four piecewise constant states separated by three waves. The left-most and right-most waves may be either shocks or rarefaction waves, while the middle wave is always a contact discontinuity. The determination of these four piecewise constant states can, in general, be achieved only by iteratively solving nonlinear equations. Thus the computation of the fluxes necessitates a step which can be computationally expensive. For this reason much attention has been given to approximate, but sufficiently accurate, techniques. One notable method is that of Harten, Lax, & Van Leer (1983; HLL), in which the middle wave, and the two constant states that it separates, are replaced by a single piecewise constant state. One benefit of this approximation, which smears the contact discontinuity somewhat, is to eliminate the iterative step, thus significantly improving efficiency. However, the HLL method requires accurate estimates of the wave speeds for the left- and right-moving waves. Einfeldt (1988) analyzed the HLL method and found good estimates for the wave speeds. The resulting method combining the original HLL method with Einfeldt’s improvements (the HLLE method), has been taken as a starting point for our simulations. In our implementation we use wave speed estimates based on a simple application of the relativistic addition of velocities formula for the individual components of the velocities, and the relativistic sound speed, $`a`$, assuming that the waves can be decomposed into components moving perpendicular to the three coordinate directions.

In order to compute the pressure $`p`$ and sound speed $`a`$ we need the rest frame mass density $`n`$ and energy density $`e`$. However, these quantities are nonlinearly coupled to the components of the velocity as well as to the laboratory frame variables via the Lorentz transformation:

$$R=\gamma n,$$
(7)
$$M^x=\gamma ^2(e+p)v^x,$$
(8)
$$M^y=\gamma ^2(e+p)v^y,$$
(9)
$$M^z=\gamma ^2(e+p)v^z,$$
(10)
$$E=\gamma ^2(e+p)p,$$
(11)
where $`\gamma =(1v^2)^{1/2}`$ is the Lorentz factor and $`v^2=(v^x)^2+(v^y)^2+(v^z)^2`$. When the adiabatic index is constant it is possible to reduce the computation of $`n`$, $`e`$, $`v^x`$, $`v^y`$ and $`v^z`$ to the solution of the following quartic equation:

$`\left[\mathrm{\Gamma }v\left(EMv\right)M\left(1v^2\right)\right]^2`$
$`\left(1v^2\right)v^2\left(\mathrm{\Gamma }1\right)^2R^2=0,`$ (12)
where $`M^2=(M^x)^2+(M^y)^2+(M^z)^2`$. This quartic is solved at each cell several times during the update of a given mesh using Newton-Raphson iteration.

Our scheme is generally of second order accuracy, which is achieved by taking the state variables as piecewise linear in each cell, and computing fluxes at the half-time step. However, in estimating the laboratory frame values on each cell boundary, it is possible that through discretization, the lab frame quantities are unphysical – they correspond to rest frame values $`v>1`$ or $`p<0`$. At each point where a transformation is needed, we check that certain conditions on $`M/E`$ and $`R/E`$ are satisfied, and if not, recompute cell interface values in the piecewise constant approximation. We find that such ‘fall back to first order’ rarely occurs.

### 2.2 Adaptive Mesh Refinement

The relativistic HLLE (RHLLE) method constitutes the basic flow integration scheme on a single mesh. We use adaptive mesh refinement (AMR) in order to gain spatial and temporal resolution. The AMR algorithm used is a general purpose mesh refinement scheme which is an outgrowth of original work by Berger (1982) and Berger & Colella (1989). The AMR method uses a hierarchical collection of grids consisting of embedded meshes to discretize the flow domain. We have used a scheme which subdivides the domain into logically rectangular meshes with uniform spacing in the three coordinate directions, and a fixed refinement ratio of $`\times 3`$. The AMR algorithm orchestrates i) the flagging of cells which need further refinement, assembling collections of such cells into meshes; ii) the construction of boundary zones so that a given mesh is a self-contained entity consisting of the interior cells and the needed boundary information; iii) mechanisms for sweeping over all levels of refinement and over each mesh in a given level to update the physical variables on each such mesh; and iv) the transfer of data between various meshes in the hierarchy, with the eventual completed update of all variables on all meshes to the same final time level. The adaption process is dynamic so that the AMR algorithm places further resolution where and when it is needed, as well as removing resolution when it is no longer required. Adaption occurs in time, as well as in space: the time step on a refined grid is less than that on the coarser grid, by the refinement factor for the spatial dimension. More time steps are taken on finer grids, and the advance of the flow solution is synchronized by interleaving the integrations at different levels. This helps prevent any interlevel mismatches that could adversely affect the accuracy of the simulation.

In order for the AMR method to sense where further refinement is needed, some monitoring function is required. In general we find that recognizing the presence of significant mass density gradients, contact surfaces, or strong shear is effective. The choice of function is determined by the part of the flow that has the most significance in a given study. For the simulations presented here the first two methods were employed, with the location of contact surfaces being found by comparing cell-to-cell pressure differences with cell-to-cell laboratory frame mass density differences. Since the tracking of shock waves is of paramount importance, a buffer ensures the flagging of extra cells at the edge of meshes, ensuring that important flow structures do not ‘leak’ out of meshes during the update of the hydrodynamic solution. The combined effect of using the RHLLE single mesh solver and the AMR algorithm results in a very efficient scheme. Where the RHLLE method is unable to give adequate resolution on a single coarse mesh the AMR algorithm places more cells, resulting in an excellent overall coverage of the computational domain.

### 2.3 Code Validation and Setup

The code has been validated using 1-D relativistic shock tube, 3-D relativistic shock reflection and 3-D relativistic blastwave trials (Hughes, Miller & Duncan 2001). Furthermore, the solver employed in the current code is a direct extension to 3-D (with a recast from Fortran 77 to Fortran 90) of the solver described by Duncan & Hughes (1994). Evidence for the accuracy and robustness of that code comes from, in addition to its application to test problems, the general agreement between studies performed with that code and with independently constructed codes (e.g., Martí et al. 1997).

In the simulations performed here, a ‘pre-existing’ jet flow is established across the computational volume, to represent the case in which a leading Mach disk and bow shock have passed, leaving a flow in pressure balance with a low density external (cocoon) medium – the shocked jet material. For all simulations we take the ratio of rest frame densities to be $`\rho _\mathrm{j}/\rho _\mathrm{x}=10.0`$ where $`\rho `$ is the rest frame mass density. The jet flow has $`v_\mathrm{j}=0.9165c`$ and $`\gamma \left(1\beta ^2\right)^{1/2}=2.5`$. The value of the pressure (and thus the sound speeds inside and outside the jet) is adjusted to yield a generalized Mach number for the jet of $`\gamma \beta /\gamma _\mathrm{s}\beta _\mathrm{s}=8`$. Here $`\gamma _\mathrm{s}\left(1\beta _\mathrm{s}^2\right)^{1/2}`$ where $`\beta _\mathrm{s}a/c`$ and the sound speed, $`a`$ is given by

$$a\left[\frac{\mathrm{\Gamma }p}{\rho +\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }1}\frac{p}{c^2}}\right]^{1/2}\text{ .}$$
(13)
With $`\mathrm{\Gamma }=5/3`$ as the adiabatic index, the relevant sound speeds are $`a_\mathrm{x}=0.6121c`$ and $`a_\mathrm{j}=0.2753c`$.

The computational domain is $`8R_\mathrm{j}\times 8R_\mathrm{j}\times 41R_\mathrm{j}`$, with outflow boundary conditions imposed on all surfaces except the inflow plane, $`z=`$ constant. The inflow plane involves cells cut by the jet boundary for which state variables must be established through a volume-weighted average of the internal and external values. To avoid a ‘leakage’ of jet momentum into the ambient material, fixed, initial values are used across that entire boundary plane at every time step. Schlieren-like images of the pressure, in which the gradient is rendered on an exponential map, provide a way of exploring flow structures with a large range of linear scale and amplitude; an inspection of these images at the last time step of the computations described here reveals that the computations are free of any effects due to reflection from the domain boundaries. Three levels of refinement were admitted, and as the entire jet is fully refined by the AMR algorithm initially there are 27 of the finest level cells across the jet diameter, $`2R_\mathrm{j}`$; outside the jet courser cells are employed initially.

A precessional perturbation is applied at the inflow by imposing a transverse component of velocity with $`v_{}=0.01v_\mathrm{j}`$. Simulations were performed with precessional perturbations $`\omega R_\mathrm{j}/v_\mathrm{j}=`$ (A) $`0.40`$, (B) $`0.93`$, (C) $`2.69`$. The simulations were halted after $`44`$ light-crossing-times for the jet radius, before development of structure had reached a quasi-steady state across the entire computational grid, when the AMR algorithm required refinement that exceeded the available computer memory.

## 3 Simulation Results

Prior to the data analysis we reduced data from all of the refined AMR “patches” in each relativistic simulation into uniformly spaced data. Each simulation contained roughly 300 patches at termination, with most patches at the finest grid resolution and occurring closer to the jet than the outer boundary. Since we wanted to observe the radial and azimuthal response within the jet at the finest resolution possible, our uniform grid along the transverse axes is at the finest resolution of the AMR simulations (13.5 zones/Rj). To reduce the size of the final data sets we have used a moderate scale grid (4.5 zones/Rj) along the jet axis in our uniformly spaced data. The lower resolution along the jet axis allows about 10 zones across 2Rj features in the axial direction; this length is similar to the wavelengths seen in the highest frequency simulation.

We evaluate simulation (theory) results quantitatively by taking 1-D slices through data cubes parallel to the $`z`$-axis at different radial distances along the transverse $`x`$-axis ($`y`$-axis for the theory) as indicated in Figure 1. For the simulation results shown in Figure 2 the $`v_x`$ and $`v_y`$ velocity components correspond to radial, $`v_r`$, and azimuthal, $`v_\varphi `$, velocity components in cylindrical geometry. Axial velocities near to the jet axis show that the simulations are fully evolved out to about $`30R_\mathrm{j}`$ for the low and moderate frequency simulations and out to about $`20R_\mathrm{j}`$ for the high frequency simulation.

Slowing of jet material resulting from surface effects is readily apparent. Surface effects penetrate more deeply into the jet as the precession frequency decreases, e.g., large dips in $`v_z`$ and $`v_\varphi `$ velocity components are observed closer to the jet axis for lower precession frequency. A dominant wavelength of $`\lambda /R_\mathrm{j}`$ (A) 14, (B) 6, (C) 2 is revealed in the axial velocity plots in the 1-D slices farthest from the jet axis and in the transverse velocity plots. The basic helical nature of the structure induced by precession is graphically illustrated by the out of phase oscillation in the transverse velocity components.

In all simulations the pressure structure shows oscillation clearly related to helical motion in 1-D slices farthest from the jet axis, but the structure near to the jet axis can be considerably more complicated – seen mostly in the low and moderate frequency simulations. The maximum pressure fluctuations around the local mean are less than $`\pm 15\%`$ and are smaller near to the jet axis and in the high frequency simulation. Growth in the transverse velocity components, indicative of unstable helical wave growth, is seen only in the low frequency simulation. In the moderate frequency simulation the transverse velocity component’s amplitude remains approximately constant. In the high frequency simulation the transverse velocity components show an initial rapid decline to a plateau. A subtle change in the radial velocity structure occurs at $`z/R_\mathrm{j}10`$ in the high frequency simulation. The change, seen in theoretical fits (see §4), occurs in the amplitude of the 1-D radial velocity slice near jet center relative to the amplitude of the 1-D slice near the jet surface. At larger distance the amplitude near to the center is less than that at the surface but at smaller distance the amplitude at center is larger. A relatively long length scale amplitude modulation of the plateau in the $`v_\varphi `$ velocity component seen beyond $`z/R_\mathrm{j}12`$ in theoretical fits (see §4) may be present to a small degree in the simulation data. As we shall see in §4, subtleties like this are indicative of the presence of multiple wave modes.

Figure 3 shows a spatial Fourier analysis of the 1-D pressure, axial velocity, and radial velocity slices shown in Figure 2 but windowed from $`1.888z/R_\mathrm{j}30.111`$.

All velocity power amplitudes are similarly normalized with the exception of those off the scale (included so power peaks can be seen). Pressure power amplitudes are also similarly normalized. The relatively short length of our window results in coarse wavelength coverage and power is computed at $`\lambda /R_\mathrm{j}`$ 28.2, 14.0, 9.3, 6.9, 5.5, 4.5, 3.8, 3.3, 2.9, 2.6, 2.4, 2.1, 2.0, 1.8, 1.7, etc. We expect power peaks to fall at the wavelength closest to the true wavelength. However, various tests indicated that shift in the location of power peaks to an adjacent wavelength bin could be caused by amplitude changes in the pressure and velocity oscillations, and also that the location of power peaks and their amplitude are somewhat sensitive to the window location. Still, it is apparent that the power distribution depends on the radial location within the jet. In general, power peaks in the radial velocity most nearly indicate the true wavelength. To a certain extent these power peaks in the radial velocity are accompanied by similar power peaks in the pressure and axial velocity at radial locations $`r/R_\mathrm{j}=0.55`$ & $`0.77`$. At radial locations $`r/R_\mathrm{j}=0.11`$ & $`0.33`$ power peaks in the pressure can be identified with similar peaks found in the axial velocity, but there are significant differences. Study of 1-D cuts through the data cubes and comparison with theoretical predictions (see §4) lead to the following general conclusions:

(1) A small amount of power in the radial velocity at wavelengths below the power peak suggests the presence of a second transverse oscillation. The effect is most evident at radial locations $`r/R_\mathrm{j}=0.11`$ & $`0.33`$. Power at these shorter wavelengths increases from about $`1\%`$ to about $`10\%`$ of the peak power as one goes from low to high frequency simulations. Power can also be seen to be enhanced at shorter wavelengths in pressure and axial velocity.

(2) In all simulations a significant power peak in the pressure and axial velocity, particularly apparent at radial location $`r/R_\mathrm{j}=0.11`$, occurs at $`\lambda /R_\mathrm{j}9`$. This feature can be unambiguously identified with a conical pressure wave at the inlet. This conical pressure wave appears in the 1-D pressure slices at $`r/R_\mathrm{j}=0.11`$ in Figure 2 as the pressure dip at axial distance $`z/R_\mathrm{j}78`$ followed by a pressure rise at $`z/R_\mathrm{j}910`$.

Pressure and transverse velocity vectors in planes transverse to the $`z`$-axis at locations of $`z/R_\mathrm{j}=10`$ & $`15`$ are shown in Figure 4. In these cross sections the jet flow is into the page. These locations were chosen as being far enough outwards to be relatively unaffected by inlet effects (beyond about $`5R_\mathrm{j}`$), not so far outwards as to be beyond the quasi-steady region \[beyond (A & B) $`30R_\mathrm{j}`$ or (C) $`20R_\mathrm{j}`$\] or where the jet surface layers are strongly slowed (beyond $`20R_\mathrm{j}`$), and to be far enough separated to see significant changes in the internal structure of the jet. The counterclockwise precession of the jet is revealed by spiral shock waves outside the jet with footpoints on low and high pressure regions on opposite sides of the jet. In the medium just outside the jet, fluid flows from high to low pressure regions around the jet circumference with less flow at higher precession frequency. The highest and lowest pressures and the largest transverse velocities lie at and outside the jet surface, but in general fluid motions in the external medium remain less than the sound speed. The existence of the spiral shocks suggests the potential for significant energy loss from the jet surface layers, and the development of a significant axial velocity shear layer and azimuthal velocity effects appears consistent with such an energy loss.

In order to see more subtle pressure and transverse velocity structure inside the jet we have applied an axial velocity mask with value $`1`$ if $`v_z0.90c`$ and $`0`$ otherwise to the cross sections.

The internal pressure structure shows high and low pressure regions near to the jet surface corresponding to the high and low pressure regions beyond the jet surface in the unmasked images, but with additional complicated structure that changes dramatically between the two locations. Within the jet, transverse velocities are less than in the higher sound speed external medium. In general, transverse velocity shows jet fluid moving towards the azimuthal location of maximum jet displacement at a location farther down the jet – note that spatial rotation outwards is in the clockwise sense. This direction of transverse flow basically proceeds from the leading edge of the high pressure to leading edge of the low pressure region within the jet cross sections – recall that patterns at fixed distance rotate counterclockwise in a temporal sense. Flow patterns inside the jet show some response to internal pressure structure in the high frequency simulation but are essentially straight across the jet at the lower precession frequencies.

The complexities in internal jet structure seen in 1-D slices, as suggested by the power spectrum and by transverse cuts through the jet indicate that more is going on than can be described by a simple helical twist of the jet at one wavelength. In the next section we use all of the findings above along with a theoretical description of the normal mode structures of a cylindrical jet to identify the modes leading to the structures observed in the simulations.

## 4 Jet Structure

Suppose we analyze the structures arising in the simulations by modeling the jet as a cylinder of radius $`R_\mathrm{j}`$, having a uniform density, $`\rho _\mathrm{j}`$, and a uniform velocity, $`v`$. The external medium (cocoon) is assumed to have a uniform density, $`\rho _x`$, and to be at rest. The jet is assumed to be in static pressure balance with the external medium $`p_\mathrm{j}=p_x`$. A general approach to analyzing the time-dependent structures is to linearize the fluid equations along with an equation of state where the density, velocity, and pressure are written as $`\rho =\rho _0+\rho _1`$, $`𝐯=𝐯_0+𝐯_1`$ and $`p=p_0+p_1`$, and subscript 1 refers to a perturbation to the equilibrium quantity. In cylindrical geometry a random perturbation of $`\rho _1`$, $`𝐯_1`$ and $`p_1`$ can be considered to consist of Fourier components of the form

$$f_n(r,\varphi ,z)=f_n(r)\mathrm{exp}[i(kz\pm n\varphi \omega t)]$$
(14)
where flow is along the $`z`$-axis, and $`r`$ is in the radial direction with the flow bounded by $`r=R_\mathrm{j}`$. In cylindrical geometry $`k`$ is the longitudinal wavenumber, $`n`$, is an integer azimuthal wavenumber, for $`n>0`$ the wave fronts propagate at an angle to the flow direction, the angle of the wavevector relative to the flow direction is $`\theta =tan(n/kR)`$, and $`+n`$ and $`n`$ refer to wave propagation in the clockwise and counterclockwise senses, respectively, when viewed outward along the flow direction. The dispersion relation describing the propagation and growth of the “normal” modes, $`n`$, along with the expressions giving the density, pressure and velocity structure associated with the normal modes can be found in Hardee et al. (1998) and Hardee (2000).

The dispersion relation has been solved for k($`\omega `$) using root finding techniques for parameters appropriate to the numerical simulations. Numerical solutions to the dispersion relation for the first two pinch body modes and for the surface and first two helical body modes appropriate to the simulations are shown in Figure 5. In general, the body modes are weakly damped in a frequency range just below where they are not growing (damping rates not shown in the figure).

The solutions have comparable maximum growth rates for these surface and body modes where the spatial growth rate is given by the imaginary part of the wavenumber, kI. Wavelengths, wave speeds and growth (damping) lengths, $`\mathrm{}|k_I|^1`$, for the helical modes at the precession frequencies used in the simulations are given in Table 1. The wave speed is defined by the real part of the phase velocity as $`v_w=(\omega /k)|_{\mathrm{r}eal}`$ and the wavelength is defined by $`\lambda =2\pi v_w/\omega `$. We note that non-relativistic jet simulations have shown wavelength and wave propagation to be given most accurately by the phase velocity and driving frequency and not, for example, by using the real part of the complex wavenumber or the group velocity to determine the wavelength. The wavelengths given in Table 1 along with mode structure and amplitudes deduced from the simulations are used to produce theoretical data cubes.

Producing theoretical data cubes that yield pressure and velocity structures similar to pressure and velocity structures in the numerical simulations was an iterative process. We concentrated on reproducing the general features seen in the simulation 1-D pressure and velocity slices (Figure 2), and then on the jet’s transverse pressure structure as shown in Figure 4. These features and structures along with the simulation power spectra shown in Figure 3 suggested the need for a pinch mode in addition to helical surface and helical body modes. This pinch contribution has been modeled by the first pinch body mode at the wavelength $`\lambda _\mathrm{p}/R_\mathrm{j}=8`$. This wavelength is approximately the length of the conical inlet pressure wave, is equal to the generalized Mach number, and is consistent with the power spectra shown in Figure 3. We note that the second pinch body mode exists but is not growing at this wavelength and higher order pinch body modes do not exist at this long wavelength. The 1-D radial velocity slices, which in the simulations are clearly least affected by surface effects, proved a good indicator of helical mode amplitudes. Features in the 1-D pressure slices near to the jet axis provided an indication of the relative phasing between helical surface and body modes, and, primarily in the case of the low and moderate frequency simulations, indicated the presence of pinching triggered by the conical inlet pressure wave. Fine tuning of relative phasing between modes was provided by comparing pressure structure between simulation and theoretical cross sections.

Results of a combination of helical surface and first body modes along with (A & B only) the first pinch body mode are shown as 1-D pressure and velocity slices in Figure 6.

Qualitatively the pinch mode amplitude grows rapidly to a maximum at $`z/R_\mathrm{j}=10`$, also the approximate location of the tip of the conical pressure wave on the jet axis, and declines rapidly at larger distance, the helical surface mode amplitude grows (A) slowly with distance, (B) is constant, (C) rapidly declines to a constant value, and the helical first body mode amplitude (A) declines slowly with distance, (B) is constant, (C) grows rapidly to a constant value. Quantitatively amplitudes are input as a maximum jet surface displacement as a function of $`z`$, the inlet phase of the helical surface mode, $`\varphi _s`$, is chosen to match the major radial velocity oscillations seen in the simulations, and the inlet phase of the helical body mode, $`\varphi _1`$, is specified relative to the helical surface mode to produce interference effects and cross section structure at the appropriate locations. The inlet phase of the pinch mode is set to produce a maximum on the jet axis at $`z/R_\mathrm{j}=10`$. The expressions used for the individual modes are:

(A)

A$`{}_{p}{}^{}=\{\begin{array}{cc}0.00075z& z<10\\ 0.0075\times 2^{(z/101)}& z10\end{array}`$

A$`{}_{s}{}^{}=0.130+0.005z`$

A$`{}_{1}{}^{}=\{\begin{array}{c}0.004\times e^{z/120}\\ \varphi _1=\varphi _s+0.26\pi \end{array}`$

(B)

A$`{}_{p}{}^{}=\{\begin{array}{cc}0.0004z& z<10\\ 0.004\times 2^{(z/101)}& z10\end{array}`$

A$`{}_{s}{}^{}=0.090`$

A$`{}_{1}{}^{}=\{\begin{array}{c}0.010\\ \varphi _1=\varphi _s+0.5\pi \end{array}`$

(C)

A$`{}_{p}{}^{}=0.0`$

A$`{}_{s}{}^{}=\{\begin{array}{cc}0.038(10.0770z)& z<10.8\\ 0.0065& z10.8\end{array}`$

A$`{}_{1}{}^{}=\{\begin{array}{cc}0.00060\times 2^{z/4}& z<13.3\\ 0.0060& z13.3\\ \varphi _1=\varphi _s+1.6\pi & \end{array}`$

where the $`A`$ and $`z`$ values above are normalized to the jet radius, $`R_\mathrm{j}`$, and $`\varphi _{(s,1)}`$ are the phase angles at $`z=0`$. At the higher frequencies jet surface displacement decreases and the initial phase difference between helical surface and first body mode increases.

In general, the functional forms chosen for the amplitudes are simply a best fit to what is seen in the simulations. While amplitude growth for the pinch mode is chosen to provide a reasonable emulation of the conical pressure wave at the inlet and has no other physical significance, subsequent damping may reflect poor coupling between the conical inlet perturbation and the required structure of this pinch body mode. A similar result was found by Hardee et al. (1998) in axisymmetric relativistic jet simulations. For the first helical body mode in the low frequency precession simulation we have used the exponential damping rate predicted by the theory. Typically the helical mode wave growth or damping seen in the simulations is not exponential, e.g., slow linear growth or rapid linear damping of the surface mode in the low and high precession frequency simulations, respectively. This fact implies that amplitudes seen in the simulations are in the non-linear regime. Interestingly, the initial rapid linear damping of the helical surface mode transverse velocity components in the high frequency simulation is accompanied by rapid nearly exponential growth of oscillations in the axial velocity component. In the simulation, damping of the helical surface mode ceases when the axial velocity oscillations reach the level appropriate to the transverse velocity oscillations predicted theoretically to accompany the helical surface mode at this rapid precession frequency.

While we are able to match the major features in the 1-D slices shown in Figure 2 on the innermost 1-D slices, significant differences between theoretical and simulation 1-D slices appear in the outer half of the jet where velocity shear appears in the simulations. This velocity shear modifies the radial pressure profile, the axial velocity profile and the azimuthal velocity profile relative to the theoretical model. First, the pressure oscillations in the simulation do not attain the amplitudes seen in the theoretical profiles in the outer half of the jet. Second, the azimuthal velocity component oscillations can be greatly exaggerated in the simulations in the outer half of the jet. This second effect becomes more pronounced at lower precession frequency and at the accompanying larger sideways jet displacement. Third, a larger indication of beating between helical surface and first body modes appears in the 1-D theoretical slices. Nevertheless, results based on the linear theory appear to provide reasonable estimates of the pressure and velocity fluctuations observed in the simulations.

Theoretical transverse velocity and pressure structure of individual helical surface and first body modes at the wavelengths given in Table 1 are shown in the top six panels in Figure 7. In these panels jet flow is into the page and spatial rotation of the patterns down the jet is in a clockwise sense. At fixed position temporal rotation of the patterns is in a counterclockwise sense. For the surface mode at the fixed spatial location of the slices, transverse flow is approximately from the leading edge of the high pressure region to the leading edge of the low pressure region. Flow patterns for the first body mode at the longer wavelengths are identical and show flow towards (away from) the regions that will become high (low) pressure regions farther down the jet.

Results of the combined modes are shown in the bottom six panels in Figure 7. These panels show a theoretical fit to the pressure structure seen in the cross sections at $`z/R_\mathrm{j}=10`$ & $`15`$ in the simulations. Because we matched theoretical 1-D slices at radial positions along the transverse $`y`$-axis (see Figure 1) to simulation 1-D slices along the transverse $`x`$-axis, cross sections through the theoretical data cubes are rotated by $`90^{}`$ relative to the simulation cross sections at comparable positions. The theoretical cross sections are taken from a theoretically generated data cube at axial locations of $`z/R_\mathrm{j}10`$ and $`z/R_\mathrm{j}`$ (A) 16, (B) 15, (C) 14. Small differences in mode wavelengths between the theoretical predictions and the simulations lead to the differences in the outermost location of the best fit. We are able to reproduce the basic pressure structure seen in the simulation cross sections and the results prove that simulation structures are the result of a combination of normal modes although there are obvious differences.

In part, differences between simulation and theory are the result of the velocity shear layer that occurs in the simulations but is not included in the theoretical model. We see that transverse flow vectors are rotated relative to the pressure structure by up to about $`45^{}`$ between theoretical cross sections and simulation cross sections. Also we find that pressure structure in the cross sections is very sensitive to the relative phasing between wave modes. For example, the theoretical cross section for simulation A at $`z/R_\mathrm{j}10`$ requires that the pinch mode provide a central pressure enhancement while the high pressure region associated with the surface and first body modes must be very close to $`180^{}`$ out of phase. Only this phasing along with appropriate amplitudes provides a pressure cross section with characteristics similar to that of the simulation. With the exception of this cross section only the two helical modes appear important to other cross section structures. In general the body mode becomes stronger relative to the surface mode at higher precession frequency. Note that transverse velocity structure appears significantly influenced by the body mode only for the high frequency case, as in the corresponding simulation cross section.

Figure 8 shows a spatial Fourier analysis of the theoretical 1-D pressure, axial velocity, and radial velocity slices shown in Figure 6. A window from $`z/R_\mathrm{j}230`$ was used for low and moderate precession frequencies and $`z/R_\mathrm{j}129`$ was used for the high precession frequency. The length of our theoretical window and the point spacing is identical to that used for the simulations, although the normalization is different. Not surprisingly a Fourier analysis of the radial velocity slices shows the best correspondence between theory and simulation as the simulation radial velocity slices are the least influenced by surface effects. Note that relative peak power levels in the radial velocity in the three simulations are mirrored by the theoretical results. The Fourier analysis of pressure and axial velocity slices shows much more variability although similar power peaks are evident in theory and simulation results. We note that normalization, particularly for the pressure slices, is very different for theory and simulation as a result of the large difference in the amplitude of the pressure fluctuation near to the jet surface. Other differences between theory and simulation are the result of different jet structure both near and far from the inlet in the simulations. If factors such as these are considered, the Fourier analysis indicates good agreement between the theoretical models and the simulations.

## 5 Synthetic Emission Images

The extent to which complexities in jet structure yield observable consequences is indicated by the plane of the sky line-of-sight integrations of p2 shown in Figure 9. Simulation images are constructed by integrating over only those computational zones in which $`v_z0.90c`$. Intensities in the different images can be qualitatively but not quantitatively compared as the scales are not identical. Flow and pressure patterns take about $`5R_\mathrm{j}`$ to develop in the simulations so simulation and theoretical images differ on this scale. Beyond this distance simulation and theoretical images are remarkably similar.

The low frequency simulation shows a long wavelength sinusoidal oscillation at the principal wavelength indicated by the 1-D slices in Figure 2. Growth in the amplitude of oscillation is apparent in the image. Decrease in brightness along the jet in the simulation image is a manifestation of the average $`15\%`$ pressure drop between the inlet and $`z/R_\mathrm{j}=30`$ (see Figure 2). There is a modest enhancement in brightness where peaks appeared in the 1-D pressure slices in Figure 2 and at the locations of maximum jet displacement in the plane of the sky. The corresponding theoretical image is shifted by the expected quarter wavelength relative to the simulation image but otherwise exhibits similar structure. Modest brightness enhancement at maximum displacement is a combination of location and shape of the high pressure region within the jet and line-of-sight effects. The effect of pinch overpressure at $`z/R_\mathrm{j}=10`$ is more apparent

in the theoretical image. The primary difference between the two images is a consequence of the jet’s axial pressure decline in the simulation. The presence of the helical body mode, which is at low amplitude, has no apparent influence on the appearance of the jet in simulation or theoretical images.

The moderate frequency simulation shows a sinusoidal oscillation at the principal wavelength indicated by the 1-D slices, but note the lengthening in the apparent wavelength only between $`z/R_\mathrm{j}=1421`$. Some decrease in brightness results from an average $`10\%`$ pressure drop along the jet in the simulation. In this simulation significant brightness enhancement at locations of the maximum transverse displacement is apparent along much of the jet. As suggested by the 1-D velocity slices there is little growth in the oscillation amplitude. The corresponding theoretical image shows the same brightness enhancements at maximum displacement. There is more enhancement than in the low frequency images. It is readily apparent that the theoretical image reproduces the simulation image structure between $`z/R_\mathrm{j}=1421`$. From the theoretical model we deduce that this structure is the result of interference between helical surface and first body modes.

The high frequency simulation shows two distinct sinusoidal oscillation patterns separated by a transition region centered at $`z/R_\mathrm{j}=1213`$ although the oscillation wavelength is nearly identical on either side of the transition. A change in transverse velocity structure occurred just previous to this transition distance and can be seen in the 1-D slices in Figure 2. Similar structure is also apparent in the theoretical image and in the 1-D slices shown in Figure 6. From the theoretical model we deduce that this transition region and change in apparent structure results from decline in the helical surface mode amplitude and growth in the helical body mode amplitude and interference between the two modes. If the theoretical image is extended to larger distance with constant mode amplitudes a beat pattern emerges as a result of the slight difference in wavelength between the two wave modes. Recall that the beat pattern is apparent in the 1-D theoretical pressure and velocity slices shown in Figure 6 and also is suggested in the 1-D simulation velocity slices shown in Figure 2.

## 6 Summary and Discussion

In the simulations, spiral shock waves are driven into the external medium by jet helicity. While fluid motions in the external medium remain less than the sound speed the waves propagate outwards at or slightly above the sound speed. The footpoints of these waves are tied to helical “ripples” in the jet surface which move axially with speeds $`0.86c>v_z^w>0.74c>a_\mathrm{x}=0.61c`$. The azimuthal motion of the high pressure footpoint at fixed axial position increases from (A) low to (C) high precession frequency with $`v_\varphi ^w`$ (A) $`0.37c`$, (B) $`0.85c`$, and (C) $`2.46c`$. In the high frequency simulation the high pressure footpoint moves around the jet circumference superluminally and would appear to move across the jet diameter at about $`1.6c`$. The existence of these spiral shocks suggests the potential for significant energy loss from the jet surface layers. In fact the simulations show the development of a significant axial velocity shear layer and azimuthal velocity effects consistent with such an energy loss. Energy loss appears greatest (as judged by axial and azimuthal velocity effects) at the low precession frequency where physical displacement of the jet surface is largest and least in the high precession frequency simulation where physical displacement of the jet surface is smallest.

Within the jet, fluid motions with respect to the helical pattern speed are less than the jet sound speed, pressure fluctuates less than $`\pm 15\%`$ around the local mean and velocity fluctuations in the jet fluid are much less than the jet sound speed. Axial velocity fluctuation is a very small fraction of the relativistic jet speed, and significant variation in the Lorentz factor does not occur in these relativistic jet simulations. The small transverse velocity induced by jet helicity allows for some angular variation in the flow direction. Angular variation decreases from $`\pm 1.25^{}`$ in the low frequency simulation to about $`\pm 0.3^{}`$ in the high frequency simulation. The angular flow variations seen in these simulations are less than 10% of the beaming angle given by $`1/\gamma `$ and given the relatively small variation in Lorentz factor we would not expect significant Doppler boosting fluctuations at small angles to the line-of-sight. However, we note that larger jet displacements, and larger pressure and velocity fluctuations would occur at larger distances outside the computational grid in the low frequency simulation. On the other hand, it is likely that the jet displacements, and pressure and velocity fluctuations in the high frequency simulation are close to saturation. Larger jet displacements and larger pressure and velocity fluctuations cannot be ruled out for the moderate frequency simulation beyond the length of the computational grid, although fluctuations remain relatively constant across the computational grid.

Details of internal jet pressure and velocity structure can be understood as arising from a combination of the normal modes predicted by theory. In general it is possible to provide good estimates of the velocity and pressure fluctuations in the jet interior but not near the jet surface where there is significant velocity shear in the simulations. Fitting jet structures requires a combination of helical surface and first body modes with a larger relative amplitude for first body mode as the precession frequency increases and where first body mode is more rapidly growing. Typically the helical mode wave growth or damping seen in the simulations is not exponential and this fact implies that amplitudes seen in the simulations are in the non-linear regime. Comparison with theory shows that the initial precession triggers the first helical body mode in addition to triggering the surface mode even if the body mode is not growing or is weakly damped. In the simulations some pinching is observed associated with a conical pressure wave at the inlet. The very oblique pressure wave induced at the inlet cannot couple strongly to an allowed normal pinch body wave whose structure includes a much less oblique pressure wave. Thus, we would expect damping of the initial perturbation beyond the first maximum as is suggested by theoretical fits to the simulations. A similar result was previously found by Hardee et al. (1998) for axisymmetric jets with much higher Lorentz factors.

In the line-of-sight images the sinusoidal oscillation becomes more confined to the jet interior as the precession frequency increases, and the influence of the body mode is enhanced. The high pressure region is somewhat ribbon like in cross section and this leads to enhancement in line-of-sight images at the maximum transverse displacement of the high pressure region. We note that maximum transverse displacement of the high pressure region is shifted axially relative to the maximum surface displacement. We observe additional structure within the jet in addition to the basic sinusoidal oscillation. The image structure can only be adequately understood with the modeling capability afforded by the theory. In line-of-sight images internal structure arises because of the configuration and location of the high pressure region associated primarily with the interacting helical surface and body wave modes. The modes are triggered with some initial phase difference depending on frequency but the phase difference between the modes at the inlet remains constant. In the low frequency simulation the body mode amplitude is small relative to the surface mode amplitude and there are no observable consequences of wave-wave interaction between the surface and first body mode. In the moderate and high frequency simulations where body mode amplitudes are significant, wave-wave interaction between surface and body modes does have observable consequences. In the moderate frequency simulation, line-of-sight integrations reveal the effects of constructive and destructive interference between these wave modes. In the high frequency simulation the surface mode declines and body mode grows to comparable amplitude. The image reveals interference effects similar to those seen in the moderate frequency simulation and shows that the body mode produces effects more confined to the jet interior than the surface mode. In any event we see that wave-wave interaction has observable consequences that are frequency dependent. Additionally, we find that the differences between simulation and theoretical models, particularly in the pressure fluctuations in the outer portion of the jet, do not result in significant differences in line-of-sight images.

Both simulations and theory suggest potentially interesting wave-wave interactions resulting from beating between wave modes. Since the simulations show a fixed phase difference between modes at the inlet, the different mode wavelengths result in distinct wave interaction regions (regions of constructive and destructive interference) that will remain stationary over time although individual wave patterns will move through these regions at different wave (pattern) speeds. The present simulations show that these pattern speeds can be quite different. Similar wave-wave interactions have been found also in non-relativistic numerical simulations (Xu, Hardee & Stone 2000) in a totally different parameter regime relevant to protostellar jets. Thus, this effect requires no fine tuning of parameters and should produce observable consequences on astrophysical jets in general. In the present context one might ask if the presence of moving and fixed components, or quasi-periodic spatial variation in brightness along jets indicates regions of destructive and constructive interference between normal modes excited close to the central engine. In particular, could wave-wave effects lead to the quasi-periodic spacing of brighter regions in the M 87 jet? For example, wave-wave interference could produce particularly complicated structures in a fixed or slowly moving knot, say in knot D in the M 87 jet, but with motions within the knot representative of a combination of pattern and flow speeds with very different apparent velocities. Future work will be specifically designed to address these issues.

P. Hardee, A. Rosen and E. Gomez acknowledge support from the National Science Foundation through grant AST-9802955 to the University of Alabama. P. Hughes acknowledges support from the National Science Foundation through grant AST-9617032 to the University of Michigan.


format: text

text:
# Computer simulation of fluid phase transitions

## I Introduction

One of the fundamental tasks of statistical mechanics is to forge the link between the microscopic (atomic-scale) description of a particular substance and its equilibrium macroscopic (thermodynamic) properties. Typically the former is prescribed in terms of a model, ie. a specification of the constituent molecules and their mutual interactions. Given this, the challenge is to derive the associated macroscopic properties —quantities such as the compressibility, heat capacity and, above all, the phase behaviour i.e. the conditions under which the substance forms a solid, liquid or gas. Computer simulation allows one to do this.

In many respects, a simulation can be viewed as an experiment performed not on a real substance, but on a model system stored in a computer’s memory. As in real experiments, one has to prepare a properly equilibrated sample under the desired thermodynamic conditions; and, just as in real experiments, one can measure the physical properties of that substance, perhaps leading to new discoveries as a result! But, a simulation can have important advantages over a real experiment. Since one is dealing with a numerical model, substances can be studied that are too expensive, too complicated or too dangerous to be tackled by real experiment. Furthermore, since the simulator has access to full and complete information about the state of the simulated system, there are fewer restrictions on just which properties can measured. Accordingly, information and insight can be gleaned from a simulation which would not readily be obtainable by experimental means .

However simulations do have their limitations. The chief drawbacks are constraints on computer speed and memory; most contemporary computers can deal only with systems comprising a few thousand particles –many orders of magnitude fewer than the $`10^{23}`$ typically found in experimental samples. Such restrictions engender so-called finite-size effects in the results, ie. spurious artifacts and systematic discrepancies (compared to the bulk limit) the magnitude of which need to be quantified. It should also be borne in mind that irrespective of the sophistication of the simulation techniques employed, results for macroscopic equilibrium behaviour will reflect reality only to the extent that the underlying model correctly captures the true microscopic nature of the substance under study. As the old computing adage goes: rubbish in, rubbish out.

Simulation strategies for obtaining the equilibrium phase behaviour of classical fluid systems come in a profusion of different shapes and forms, but broadly speaking fall into two main categories: Molecular Dynamics (MD) and Monte Carlo (MC). Both have been previously discussed extensively in a number of introductory texts and articles, see eg. refs. . The MD approach involves computing the phase space trajectories of a system of mutually interacting particles by integrating Newton’s equations of motions through time. Physical properties are measured in terms of configurational or time averages as the simulation evolves. MD represents an attractive method in situations where one is interested in obtaining dynamical information about a system, and can also be employed to obtain single phase thermodynamical properties. However for studies of phase transitions, ie. the process by which one phase spontaneously transforms into another, it is rarely a suitable approach because of the problem of hysteresis (superheating and supercooling) wherein the temperature and pressure at which the transition occurs in a simulation may differ significantly from that of the bulk system. This matter is described in detail in sec. II.

Phase transitions (it now seems quite generally agreed) are best tackled by Monte Carlo (MC) simulation. Here one employs a stochastic Markov process to generate a sequence of equilibrium configurations of the model system of interest; physical properties are measured as configurational averages over the sequence . The actual mechanisms by which the system evolves from one equilibrium configuration to the next are essentially artificial, so information about physical dynamical processes is strictly limited. Nevertheless this loss is compensated for by potential gains in a variety of other areas. Specifically, freed from the strictures imposed by physical dynamics, the simulator is at liberty to construct any number of elaborate schemes via which the simulation efficiently explores the space of possible model configurations. By so doing it becomes possible to bridge time and length scales which

cannot be probed using MD.

In this article we shall focus on one area in which recently developed MC simulation techniques allow one to circumvent an old problem, namely that of accurately obtaining the phase behaviour of model fluids. The approach we describe is in essence a synthesis of a number of existing simulation techniques (developed originally in the context of lattice spin models) which, when married with a standard grand canonical or constant pressure ensemble fluid simulation, provide a powerful and efficient route to phase coexistence properties. Before embarking on this description, however, it is instructive to review some key aspects of phase transition phenomenology, in particular the origins of the hysteresis effect which for many years plagued simulation studies of phase transition and which necessitated the new methodological advances.

## II Hysteresis, Interfaces and the Free Energy Barrier

The phase diagram of a typical simple substance such as Argon, is depicted in schematic form in figure 1. Depending on the imposed conditions of temperature and pressure, the substance can exist in three phases: solid, liquid or gas. The corresponding regions of the phase diagram are delineated by phase boundaries at which transitions occur from one phase to another. Notable features of this phase diagram are the triple point where three phase boundary lines meet and the liquid-gas critical point which terminates the liquid-gas phase boundary and beyond which all distinction between the liquid and the gas vanishes.

Phase diagrams such as that of fig. 1 can be interpreted within the framework of thermodynamics by appeal to the concept of free energy. A system in equilibrium will always choose its state to be that for which the free energy is least . Phase boundaries then arise naturally as that sets of points in the phase diagram for which two phases have the same free energies, being equally favoured thermodynamically. In experiments, however, a transition from one stable phase to another will not always occur exactly on the phase boundary. Instead one generally encounters ‘overshoot’ or ‘hysteresis’, whereby the actual transition point depends on the thermodynamic history of the sample.

The following gedanken experiment will help to explain the origin of hysteresis. Imagine we take a purified quantity of a fluid such as water, place it in a sealed piston-cylinder arrangement at constant atmospheric pressure, and heat it up very slowly so that it always remains close to equilibrium. During this process, the water temperature will rise until at some point it boils –transforming into steam with a concomitant large and abrupt increase in the system volume. This is an example of what is called a first order phase transition. If, however, we stop adding heat before all the water has vapourised we observe coexistence between the liquid and vapour phases i.e. a portion of the container will be occupied by the liquid phase, separated by an interface from the remainder which contains vapour (cf. fig. 2). In transforming from one phase to another, a system always passes through such mixed-phase configurations and it is these that are responsible for hysteresis.

It transpires that mixed-phase configurations possess a higher free energy than pure phase states and since this additional free energy is wholly associated with the interface itself, it is often referred to as the surface tension. Owing to their surface tension, mixed-phase configurations are thermodynamically less favourable than pure phase states at the phase boundary. This has a bearing on phase transitions such as the vaporisation of our water sample. As the water is slowly heated it eventually attains the phase boundary temperature $`T=373.15`$ K above which steam has a lower free energy than liquid water i.e. becomes the thermodynamically favoured state. Nevertheless it is possible to ’superheat’ liquid water some way beyond this temperature without it transforming to steam because, in order to do so, it must pass through the mixed-phase states of higher free energy. A similar effect occurs when cooling steam –it becomes possible to ‘supercool’ it some way below $`373.15`$ K without it liquefying. Hence the temperature at which the transition occurs is not that at which the free energies of the two phases are equal, but instead depends on the initial state of the sample and the direction and rate in which the boundary is traversed.

Much the same thing occurs in simulations of first order phase transitions. But here the problem is also intimately bound up with issues of finite-size effects and simulation timescales. To appreciate how a simulation behaves near a first order phase transition, it is instructive to consider a concrete example, namely the liquid-gas transition of a prototype model fluid –the famous Lennard-Jones fluid (LJF). Within this model, the interaction potential for two point particles separated by linear distance $`r`$ is given by

$`U(r)=4ϵ\left(\left({\displaystyle \frac{\sigma }{r}}\right)^{12}\left({\displaystyle \frac{\sigma }{r}}\right)^6\right)`$where the parameters $`ϵ`$ and $`\sigma `$ set the strength of the interaction and the length scale respectively .

One way of performing a MC simulation of the LJF is to employ the grand-canonical ensemble (GCE) in which, for a given system volume $`V`$, the total configurational energy $`E`$ and particle number $`N`$ are permitted to fluctuate stochastically, but with average values determined by prescribed values of the temperature $`T`$ and chemical potential $`\mu `$. These latter two variables span the phase diagram and by tuning their values, transitions can be induced between gas, liquid and solid phases (cf. fig. 1).

An outline of the operation of a GCE MC simulation is given in box 1. Fluctuations in the energy and particle number occur by means of particle insertions and deletions. MC updates consist of either an insertion or a deletion attempt, each of which is proposed with equal probability. For an insertion, one chooses a random position in the simulation box and calculates the energy change $`\mathrm{\Delta }E^I`$ associated with placing a new particle at that position. The trial insertion is accepted with a probability given by a so-called Metropolis rule designed to ensure that detailed balance is maintained –a necessary condition for attaining thermodynamic equilibrium :

$`p_{\mathrm{acc}}(NN+1)`$ (1)
$`=\mathrm{min}[1,{\displaystyle \frac{zV}{(N+1)}}\mathrm{exp}(\beta (\mathrm{\Delta }E^I))]`$ (2)
where $`z=\mathrm{exp}(\beta \mu )`$ with $`\beta =1/k_BT`$.

Similarly for a particle deletion, one chooses a particle at random from those currently present in the system, calculates the energy change $`\mathrm{\Delta }E^D`$ associated with its removal and then performs the removal with a probability again given by a Metropolis rule

$`p_{\mathrm{acc}}(NN1)`$ (3)
$`=\mathrm{min}[1,{\displaystyle \frac{N}{zV}}\mathrm{exp}(\beta (\mathrm{\Delta }E^D))]`$ (4)
Box 1: Outline operation of a simple GCE simulation program. For brevity only the operations for the $`x`$ coordinate of particle position vector has been given. Analogous operations apply to the y and z coordinates.

Notwithstanding its simplicity, a GCE simulation of this type can provide a highly efficient computational route to equilibrium fluid phase properties. This is especially so if, as is customary, the interparticle potential $`U(r)`$ is truncated at some cutoff radius $`r_c`$ . Contributions to $`\mathrm{\Delta }E^I`$ and $`\mathrm{\Delta }E^D`$ then arise only via local interactions and, by partitioning the simulation box into cubic cells of side $`r_c`$, the subset of contributing particles can be readily identified. Note also, that for studies of liquid-gas phase transitions one gains nothing by implementing explicit particle displacement moves in a GCE simulation. Instead these are realized implicitly by virtue of repeated particle transfers, sole use of which not only constitutes a valid algorithm (it is clearly ergodic) but, moreover, focuses the computational effort on the bottleneck for phase space evolution namely the fluctuations in the particle number density.

Let us now examine how the GCE simulation algorithm of Box 1 behaves in the vicinity of the liquid-gas phase boundary. The primary observables of interest from such a simulation are the particle number $`N`$ and the configurational energy $`E(\{𝐫\})`$ where $`\{𝐫\}`$ denotes the set of particle position vectors i.e. the configuration. The fluctuations of the particle number density, $`\rho =N/V`$, in particular provide much insight into the nature of phase coexistence in a finite-sized system. In a simulation, one can accumulate the probability density function (pdf) of the number density, $`p(\rho )`$, in the form of a histogram averaged over many independent samples of $`\rho `$. The form of such a distribution close to liquid-gas coexistence is shown in schematic form in fig. 3.

The principal feature of this distribution is its bimodal (double-peaked) character. Each of the two peaks corresponds to one of the pure phases –the low density peak to the gas phase and the high density peak to the liquid phase. The location of the liquid-gas phase boundary line in the $`\mu T`$ plane is prescribed formally by the equal peak weight criteria i.e. by the set of values of $`\mu `$ and $`T`$ for which the integrated weights (areas) under the two peaks are equal. Thus to locate liquid-gas coexistence in a simulation one must tune $`\mu `$ and $`T`$ until the measured form of $`p(\rho )`$ is doubly peaked with equal area under each peak. The problem for simulations is that to obtain accurate estimates for the relative peak weights, the simulation procedure must supply bountiful independent samples from each of the two phases, which in turn necessitates that it pass many times back and forth between them.

Unfortunately, the inter-phase route necessarily traverses the mixed-phase configurations in which regions of both phases coexist within the simulation box. Such configurations have, on account of their surface tension, an a-priori probability that is intrinsically low compared to those of the pure phase states. This, of course, is the physical origin of the trough separating the two peaks of $`p(\rho )`$ which may accordingly be regarded as a “probability (or free energy) barrier” to inter-phase transitions . The height of this barrier is taken to be the ratio of the maximum (peak) value of $`p(\rho )`$ to its minimum value in the trough. For large barrier heights, the free energy penalty associated with forming mixed-phase configurations is so high that transitions between the two pure phases can become very rare on the timescale accessible to simulation. As a consequence, the correlation time becomes large, hindering the accumulation of independent statistics on the relative peak weights and thence the accurate location of the phase boundary.

The barrier height (and hence the scale of the difficult) depends on the temperature. At the critical temperature the surface tension vanishes and with it the probability barrier . If, however, one follows the liquid-gas boundary to progressively lower sub-critical temperatures, thermal fluctuations diminish and the surface tension grows. This is manifest in a growth in the barrier height accompanied by a narrowing of the peaks of $`p(\rho )`$. These features are illustrated in fig. 4 which shows the measured form of $`p(\rho )`$ (obtained using the GCE algorithm of box 1.) for the LJF close to the phase boundary for the two temperatures $`T=0.985T_c`$ and $`T=0.965T_c`$. Also shown in fig. 4(b) is the corresponding evolution of the density expressed as a function of the number of Monte Carlo insertion/deletion attempts. For the higher temperature, the barrier height is approximately $`5`$ and inter-phase transitions are fairly frequent. On decreasing the temperature to $`0.965T_c`$, however, the barrier height increases to a factor $`30`$ and inter-phase crossings become comparatively much rarer. This reluctance to transform between the phases is, of course, none other that the hystereis phenomenon, viewed not from the standpoint of thermodnamics, but in terms of the underlying statistical mechanics .

Qualitatively similar effects occur if the linear dimension of the simulation box $`L`$ is increased. The barrier height grows because the surface tension of the interface increases proportionally to its area ($`L^2`$ for $`d=3`$). Additionally, the peaks of $`p(\rho )`$ narrow due to ‘self-averaging’ of fluctuations. Even for quite modest system sizes the barrier can be prohibitively large for inter-phase crossings to occur on the accessible simulation timescale.

Effectively then, one is caught between the ‘rock’ of wishing to minimise finite-size effects by employing a large simulation box, and the ‘hard place’ of seeking to sample on timescales exceeding the correlation time. Evidently therefore, a superior approach to bare GCE simulation is called for if one is to tackle the problem of first order phase transitions at temperatures significantly below that of criticality. Indeed, over recent years considerable effort has been invested in developing new MC simulation methods for circumventing the problems identified above. In the next sections we describe one solution to the problem which is rapidly becoming the method of choice for high resolution studies of fluid phase equilibria.

## III Beating the barrier

Approaches for dealing with the free energy barrier in simulations of phase transitions fall broadly into two categories

* Simulations without interfaces
* Biased sampling techniques

Foremost in category (i) are methods such as the Gibbs Ensemble Monte Carlo (GEMC)  and Gibbs-Duhem integration , both of which have enjoyed widespread use and popularity in studies of phase transitions in model fluids. Although quite distinct in character, both these methods serve to link the pure phases thermodynamically, without traversing mixed-phase configurations. Both are versatile and fairly easy to use. However, they also have significant drawbacks as discussed in Sec. V. Specifically, it seems difficult using the GEMC method to obtain coexistence data of high statistical quality without a large investment of computational effort. By comparison, Gibbs-Duhem integration is more efficient, but potentially suffers from integration errors rendering it difficult to assess the accuracy of results for phase boundary properties.

More recently, an alternative approach has emerged  –one which although arguably less straightforward to implement, offers the rewards of considerably greater efficiency, precision and flexibility than GEMC or Gibbs-Duhem integration. This scheme is based on a synthesis of two existing simulation techniques: Multicanonical biased sampling and Histogram Reweighting. In the remainder of this article, we shall describe how this combined method operates and set out in a step-by-step fashion how one goes about implementing it in practice.

### A Multicanonical Sampling

Multicanonical MC owes its origin to the biased sampling techniques first introduced in the 1970’s by Torrie and Valleau to calculate free energies . Latterly however, such techniques have gained fresh impetus with the realisation  that they permit the bridging of the free energy barrier at a first order phase boundary. In this context the term “Multicanonical sampling” was coined and the method applied successfully to the study of phase transitions and free energy landscapes in a variety of lattice-based spin systems .

The basic idea underpinning Multicanonical MC is to preweight the sampling of configuration space in such a way as to artificially enhance the occurrence of the mixed-phase configurations of intrinsically low probability. By so doing it is possible to overcome the probability barrier separating the two pure phases, thereby allowing the simulation to pass unhindered between them. The result is a great reduction in the correlation time of the sampling process.

Within the GCE framework, the biasing is achieved by use of a preweighting function incorporated into the Metropolis acceptance criteria for particle insertions and deletions (cf. Box 1). Its purpose is to modify (with respect to standard Boltzmann statistics) the probability with which configurations of the various densities are visited, in such a way that the measured number density pdf is approximately flat over the entire density range separating the two pure phases. Of course the results of such biased simulations deviate from Boltzmann statistics and consequently lack direct physical significance. Nevertheless, it is possible (as we describe below) to unfold from the simulation results the unwanted effects of the imposed bias, thence recovering the physically relevant quantities which would have been obtained in an unbiased simulation had sufficient run-time been available. In general, however, there is a price to be paid for this gain, and that is the expenditure of effort involved in finding a suitable form for the preweighting function. Fortunately it transpires that for the purposes of tracing phase boundaries this is not always necessary, it being possible to obtain a suitable preweighting function for ‘free’.

#### 1 Formalism

Let us begin by considering the GCE form of the particle number pdf, $`p(N)`$, at inverse temperature $`\beta =1/k_BT`$ and chemical potential $`\mu `$. For a system of volume $`V`$ this takes the form of a simple average of the Boltzmann factor over all possible particle positions :

$$p(N|V,\beta ,\mu )=\frac{1}{Z}\underset{i=1}{\overset{N}{}}\left\{_V𝑑r_i\right\}e^\beta ,$$
(5)
where $``$ is the configurational Hamiltonian given by $`(\{𝐫\},N)E(\{𝐫\})\mu N`$, while $`Z=Z(\beta ,\mu )`$ is the partition function which serves to normalise the distribution to unit integrated weight.

The multicanonical method operates by sampling not from a simple Boltzmann distribution with Hamiltonian $`(\{𝐫\},N)`$, but from a modified distribution with effective Hamiltonian

$$\stackrel{~}{}=+\eta (N),$$
(6)
where $`\eta (N)`$ is a preweighting function defined on the set of particle numbers $`N`$ . The associated particle number distribution is given by

$$\stackrel{~}{p}(N|V,\beta ,\mu ,\eta (N))=\frac{1}{\stackrel{~}{Z}}\underset{i=1}{\overset{N}{}}\left\{_V𝑑r_i\right\}e^{\beta \stackrel{~}{}}$$
(7)
Let us now suppose for the sake of argument that we are able to choose the preweighting function such that $`\eta (N)=\mathrm{ln}p(N)`$ where $`p(N)`$ is the desired Boltzmann density distribution. Inspection of eqs. 57 reveals that this implies $`\stackrel{~}{p}(N)=\mathrm{constant}N`$. To the extent that such a choice of weight function can actually be realised, the density then performs a $`1`$-d random walk over its entire domain, thereby allowing extremely efficient accumulation of the preweighted histogram $`\stackrel{~}{p}(N)`$.

Unfortunately, this happy state of affairs cannot in general be immediately achieved because the preweighting function $`\eta (N)=\mathrm{ln}p(N)`$ that serves to flatten $`\stackrel{~}{p}(N)`$ is, of course, just the logarithm of the function we are trying to find! Means must therefore be found to obtain a form for $`\eta (N)`$ that approximates $`\mathrm{ln}p(N)`$ sufficiently well that inter-phase transitions occur with an acceptably high frequency. More refined forms for $`\eta (N)`$ can thereafter be obtained in the further course of the study.

Let us assume that a suitable preweighting function has been found, and a simulation performed to obtain good statistics for $`\stackrel{~}{p}(N)`$. The next step is to infer the distribution $`p(N)`$ we actually seek by unfolding the effects of the multicanonical preweighting. This is achievable because knowledge of the preweighting function tells us exactly by what degree the relative probabilities of the states of various $`N`$ were altered in the simulation with respect to the true Boltzmann statistics. One therefore need only divide out the relative probability enhancements from $`\stackrel{~}{p}(N)`$ to yield $`p(N)`$. This is done for each value of $`N`$ in the range of interest by the simple reweighting:

$$p(N|V,\beta ,\mu )=e^{\eta (N)}\stackrel{~}{p}(N|\beta ,\mu ,\eta (N)).$$
(8)
The details of the practical implementation of this procedure are described in Sec. III C.

### B Histogram reweighting

Histogram Reweighting is the second ingredient in our simulation procedure. It rests upon the observation that histograms of observables accumulated at one set of model parameters (in our case $`\beta `$ and $`\mu `$) can be analysed to provide estimates of histograms appropriate to other values of these parameters. Consider the joint probability distribution of energy and particle number fluctuations at some particular parameter values $`\beta =\beta _0`$ and $`\mu =\mu _0`$. Formally this is given by

$`p(N,E|V,\beta _0,\mu _0)`$ (9)
$`={\displaystyle \frac{1}{Z_0}}{\displaystyle \underset{i}{\overset{N}{}}}\left\{{\displaystyle _V}𝑑r_i\right\}\delta (EE(\{r_i\}))e^{\beta _0_0},`$ (10)
where $`_0(\{𝐫\},N)E(\{𝐫\})+\mu _0N`$. It is easy to show  that an estimate for the form of $`p(N,E)`$ at some other parameters $`\beta =\beta _1,\mu =\mu _1`$ can be obtained from the measured pdf by the simple reweighting:

$$p(N,E|V,\beta _1,\mu _1)=\frac{Z_1}{Z_0}e^{(\beta _1_1\beta _0_0)}p(N,E|V,\beta _0,\mu _0),$$
(11)
where the ratio $`Z_1/Z_0`$ is an unimportant constant which is effectively absorbed into the normalisation. If desired, this joint distribution can then be marginalised to yield a uni-variable distribution, eg. the particle number pdf at $`\beta _1,\mu _1`$:

$$p(N|V,\beta _1,\mu _1)=𝑑Ep(N,E|V,\beta _1,\mu _1).$$
(12)
Hence, in principle at least, a single simulation at one state point in the phase diagram suffices to obtain information concerning all other state points. Unfortunately the reality of the situation is less auspicious. Owing to finite sampling time, it is not possible in practice to reweight a single histogram obtained from a simulation at some $`\beta _0,\mu _0`$ to arbitrary values of $`\beta _1,\mu _1`$. Instead the parameters to which one extrapolates must be close (in a sense we shall describe) to those at which the simulation was actually performed, otherwise the procedure loses accuracy.

The problem is traceable to the fact that the reweighting represented by eq. 11 may drastically modify the relative statistical weights of the various members of the set of configurations that contribute to the spectrum of measurements. Specifically, difficulties arise with that subset of sampled configurations having very low Boltzmann weight at the simulation parameters $`\beta _0,\mu _0`$, members of which consequently occur only very rarely in the sample. For expectation values of observables calculated at $`\beta _0,\mu _0`$, contributions from this rare subset of configurations do not contribute disproportionately to statistical uncertainties. However, under the reweighting to $`\beta _1,\mu _1`$, the configurations in question may be assigned a much greater statistical weight, one which does not reflect their actual representation in the overall sample. This has the effect of magnifying the overall statistical error on measured expectation values and is manifest as a reweighted histogram that appears ‘ragged’ in its extremal regions.

One way of dealing with this problem is to perform a sequence of separate simulations at strategic intervals across the range of model parameters of interest. Typically the intervals are chosen such that histogram of some observable (eg. the energy) accumulated at neighbouring state points in the sequence overlap within some region of their domain. The role of Histogram Reweighting to then to interpolate into the regions of parameter space between the simulation points. In this context it should be noted that it is possible to combine (in a self consistent fashion) the results of a number of different simulations at different model parameters and perform Histogram Reweighting on the aggregate data. For a description of this more sophisticated procedure we refer the reader to refs. .

### C Stitching together the pieces

So how does one implement the above formalism in a GCE simulation of a fluid? In fact the task falls naturally into two parts: performing the actual multicanonically preweighted simulation and the subsequent data analysis. We consider them in turn.

The business of implementing the multicanonical preweighting within a GCE simulation is basically quite straightforward. Assuming an appropriate set of multicanonical weights has been found, all one need do is read them in and store them in an array. The simulation then proceeds as outlined in Box 1., except that the Metropolis acceptance probabilities for insertion and deletion (eqns. 2 and 4) are modified to read:

$`p_{\mathrm{acc}}(NN+1)`$ (13)
$`=\mathrm{min}[1,{\displaystyle \frac{\eta (N)}{\eta (N+1)}}{\displaystyle \frac{zV}{(N+1)}}\mathrm{exp}(\beta (\mathrm{\Delta }E^I))],;`$ (14)
$`p_{\mathrm{acc}}(NN1)`$ (15)
$`=\mathrm{min}[1,{\displaystyle \frac{\eta (N)}{\eta (N1)}}{\displaystyle \frac{N}{zV}}\mathrm{exp}(\beta (\mathrm{\Delta }E^D))].`$ (16)
Consider now the simulation quantities we seek to obtain, namely the probability distributions of $`N,E`$ and any other observables of interest. Obviously it is tempting to accumulate these distributions in the form of histograms built up in the course of the simulation by simply bining successive measurements into an array. But for continuous variables such as the energy, this strategy necessitates a prior choice for the histogram bin-width. Should a choice be made that subsequently turns out to be unsatisfactory, then one is faced with little alternative but to repeat the whole simulation. A superior approach, retaining complete information, involves decoupling the data analysis from the simulation by recording the full history of raw data measurements. This data is then postprocessed by a separate analysis program. Although such an approach can make for large simulation output files, it has the overriding advantage of ensuring maximum flexibility in terms of data analysis.

To facilitate the post processing approach, the raw data should be accumulated in the form of a list. Suppose we perform a GCE MC simulation of the LJF at a given $`\beta =\beta _0`$ and $`\mu =\mu _0`$, employing some chosen weight function $`\eta (N)`$. As the simulation proceeds we make a succession of measurements (performed at regular spaced intervals of time) of the observables $`E`$ and $`N`$, together with any other quantities of interest. Successive measurements of these observables are gathered into a list by appending them to a file, viz:

$`\{E_0,N_0,O_0,..\}`$

$`\{E_1,N_1,O_1,..\}`$

$`\{E_2,N_2,O_2,..\}`$

.

.

$`\{E_j,N_j,O_j,..\}`$

.

.

$`\{E_M,N_M,O_M,..\}`$

where $`j`$ indexes the series of $`M+1`$ measurements and $`O`$ denotes an observable of interest. The data list is analysed following the actual simulation by a separate post-processing program, the task of which is three fold. Firstly it should remove from the data the unwanted effects of the preweighting in order to recover the desired Boltzmann distributed statistical properties. Secondly, it should output pdf’s of the observables of interest in the form of histograms. Thirdly, it should (if desired) reweight the data to obtain estimates of histograms appropriate to parameter values different from those at which the simulation was performed.

It turns out that the operations of histogram reweighting and bias removal can be accomplished simultaneously because their mathematical structures are very similar. To do this, one simply runs through the data list assigning each entry $`j`$ a statistical weight

$$w_j=e^{[(\beta _1\beta _0)E_j+(\beta _1\mu _1\beta _0\mu _0)N_j+\eta (N_j)]},$$
(17)
where $`\beta _1,\mu _1`$ are the parameters to which one wishes to extrapolate. The complete set of $`M+1`$ weights $`w_0,w_1,,,w_M`$ is then used to construct the reweighted histogram for some measured observable of interest $`O`$:

$$H(O|V,\beta _1,\mu _1)=\underset{j=0}{\overset{M}{}}w_j\delta (OO_j).$$
(18)
This histogram, once suitably normalised, constitutes a discrete estimate for the probability distribution $`p(O^{}|V,\beta _1,\mu _1)`$. Implicit in its construction is the specification of the bin-width, which may need to be tuned to strike a balance between resolution and data smoothness. But with the raw list safely stored in a file, this is something that can be done very quickly.

## IV Tracing coexistence curves

We now turn to the actual simulation procedure by which a fluid phase boundary can be determined. For illustrative purposes, we shall remain with our prototype example, the liquid-gas transition of the LJF. The essential approach is nevertheless rather general and can be applied to other types of phase transition such as demixing transitions in fluid mixtures.

Clearly one of the key components of our procedure is multicanonical preweighting, use of which generally entails a degree of preliminary effort in order to determine a suitable preweighting function. It transpires, however, that provided one begins tracing the phase boundary from near the critical point which marks its terminus, no additional effort need be expended determining preweighting functions. Instead these can be obtained for ‘free’ by virtue of histogram reweighting as we shall now describe.

One starts off by gauging the approximate location of the critical point from a series of short runs on a small system. This is not time consuming and can be performed interactively at the computer, provided one arranges that the data list (see section III C) is delivered directly to the screen. One starts by nominating a temperature $`T`$ and a large negative value of $`\mu `$ at which a short simulation is performed without multicanonical preweighting. The fluctuating number density will typically settle down quite rapidly, and its average value can be estimated visually from the data output. One then repeats this procedure for a succession of progressively larger $`\mu `$ values. As $`\mu `$ is increased, the average particle number will be observed to increase steadily. However on traversing the hysteresis-shifted phase boundary, a sudden jump will occur in the density. If this jump is large, eg. for the LJ fluid, from $`\rho =0.05`$ to $`\rho =0.6`$ (in units of $`\sigma ^3`$), then the temperature is well below the critical temperature and one should increase $`T`$ somewhat and begin again. If the jump is somewhat smaller e.g. from $`\rho 0.2`$ to $`0.4`$, then one has obtained an estimate for a near-critical point on the phase boundary $`\mu _\sigma (\beta )`$. Of course, if no jump in density is observed at all then the temperature probably exceeds the critical temperature and should be reduced.

One next performs a longer run for some larger system of interest at the estimated near-critical phase boundary point, let us call it $`\mu _\sigma (\beta _0)`$. Since the surface tension and the associated barrier to inter-phase crossings are low near criticality, it should be possible to accumulate an accurate form for $`p(N,E)`$ (including information on states ‘between the peaks’) without resort to multicanonical preweighting. Having done this, the next step is to reweight the data accumulated from this run to obtain an estimate of the form of $`p(N)`$ at some lower temperature point on the phase boundary . This is achieved by first choosing an extrapolation temperature $`\beta _1`$ inside the range of reliable reweighting (so that the reweighted distribution $`p_{ex}(N|\beta _1,\mu _1)`$ appears smooth). One then tunes $`\mu `$ within the histogram reweighting scheme until $`p_{ex}(N|\beta _1,\mu _1)`$ is bimodal with equal area under each peak. This tuning procedure can be easily automated within the analysis program to deliver precise values of the phase boundary chemical potential $`\mu _1=\mu _\sigma (\beta _1)`$.

The reweighted phase boundary histogram $`p_{ex}(N|\beta _1,\mu _1)`$ will, on account of the lower temperature, be more strongly peaked than that from which it derives. Thus multicanonical preweighting will probably be necessary for a new simulation at $`\beta _1,\mu _1`$. Fortunately however, a suitable preweighting function is already to hand –it is just the extrapolated function $`p_{ex}(N|\beta _1,\mu _1)`$. Thus all one need do is set $`\eta (N)=p_{ex}(N|\beta _1,\mu _1)`$ and perform a multicanonical simulation at $`\beta _1,\mu _1`$ to obtain the actual distribution $`p(N,E|\beta _1,\mu _1)`$.

One then simply iterates this procedure: histogram reweighting of $`p(N,E|\beta _1,\mu _1)`$ is used to estimate a phase boundary point $`\mu _2=\mu _\sigma (\beta _2)`$ at temperature $`\beta _2`$, together with the extrapolated distribution $`p_{ex}(N|\beta _2,\mu _2)`$. The latter serves as a preweighting function for a further simulation at $`\mu _2,\beta _2`$, and so on. In this way one steps down the coexistence curve, obtaining at the same time the locus of the phase boundary $`\mu _\sigma (\beta )`$ and the associated set of number density pdf’s. Clearly the maximum feasible step size is set by the range of reliable histogram extrapolation, which does decreases with increasing system size. In practice, however, quite large strides can be made even for large systems. For example, in reference , a system of approximately $`600`$ LJ particles was studied and $`7`$ steps were required to reach a subcritical temperature of $`T=0.85T_c`$.

The results of implementing this procedure for the LJF  are depicted in fig. 5. The forms for $`p(N)`$ shown have equal weight in each peak and hence lie on the phase boundary. The associated phase diagram $`\mu _\sigma (\beta )`$ is shown in figure 6. The values of the coexisting densities can be simply read off from the peak positions in fig. 5(a). The enormity of the probability barrier that multicanonical preweighting allows one to negotiate is revealed by plotting these distributions on a log scale, fig. 5(b).

Finally in this section we point out that the measured coexistence form of the density distribution $`p(\rho )`$ permits an estimate of the surface tension $`\gamma `$. For a cubic system of volume $`L^3`$, this is found simply from the ratio of values of $`p(\rho )`$ at the peak and at the trough between the peaks :

$`\gamma ={\displaystyle \frac{1}{2\beta L^2}}\mathrm{ln}\left({\displaystyle \frac{p_{\mathrm{max}}}{p_{\mathrm{min}}}}\right)`$Reference  describes a recent application of this relation in a study of the surface tension of the Lennard-Jones fluid.

## V Discussion and conclusions

To summarise, we have described a scheme whereby information on the locus of a fluid phase boundary is obtainable via the dual mean of multicanonical preweighting and histogram reweighting, married to a standard grand canonical simulation algorithm. The method begins tracing the phase boundary from near the critical point where the free energy barrier to inter-phase transitions is small. Histogram reweighting of the data thus obtained is used to provide an estimate of the location of the phase boundary at some lower temperature, together with a suitable form of the requisite preweighting function. A new multicanonical simulation is performed at this lower temperature phase boundary point and the process is repeated. In this way it is possible to stride down the phase boundary, obtaining pdfs of observables such as the number density from which, in turn, the coexistence properties can be inferred. Such use of multicanonical preweighting completely eliminates hysteresis at first order phase transition.

In cases where one doesn’t wish to start tracing a coexistence curve from near the critical point, it is necessary to bootstrap the procedure described above by obtaining an initial phase boundary preweighting function. A variety of techniques exist for achieving this, ranging from simple extrapolation of the weight function into the unsampled region, to more sophisticated analyses of MC transition probabilities. Most of the techniques in common use can be straightforwardly automated. It is beyond the scope of this article to describe these methods in detail, and for further details we refer the interested reader to the literature .

Although we have illustrated our approach solely in the context of a simple fluid model, it should be noted that it is equally applicable to complex fluids such as molecules or polymers. In these systems, however, the GCE insertion probability is often small, so it is necessary to supplement the standard algorithm with a more intelligent insertion scheme –one which performs a biased choice of a molecular orientation favourable for the insertion. Methods such as Configurational Bias Monte Carlo  and Recoil Growth  allow one to do this. Apart from this added complication, multicanonical preweighting and histogram reweighting are implemented exactly as for a simple fluid.

The scheme we have described is generalisable to other simulation ensembles such as the constant-NpT ensemble in which density fluctuations occur by means of volume changes at constant particle number, pressure and temperature. Use of this ensemble can be more efficient than the GCE when dealing with very dense fluids where the success rate of particle insertions and deletions is small. The natural variable in which to preweight a constant-NpT ensemble simulation is the fluctuating volume, but otherwise the formalism is very similar to that described above. For other types of phase transitions, such as the liquid-liquid transitions occurring in binary fluid mixtures, a suitable variable in which to preweight is usually the order parameter for the transition eg. the concentration of one species.

It is instructive to compare the scheme described in this article with other commonly used methods for tracing phase boundaries. One such method is the Gibbs-Duhem integration method wherein pairs of single phase simulations (one for each phase) are performed at various state points along the phase boundary. The single phase results are connected thermodynamically, not by negotiating the free energy barrier at each state point, but via a phase space path that runs back along each side of the phase boundary to some independently-known reference point on the phase boundary. Successive simulation state points along the phase boundary are found from an integration scheme which turns out  to be a low order approximation to the histogram reweighting method. The Gibbs-Duhem method is versatile provided one has prior knowledge of a phase boundary reference point and is, for the purposes of obtaining a rough estimate of the phase boundary locus, doubtless faster than the method described in this article. However, its accuracy rests heavily upon the precision with which the boundary reference point is known. Moreover, as one departs from this point, integration errors can grow and provided one remains within the (rather wide) region of metastability bordering either side of the phase boundary line, there is no feedback to indicate that one has in fact strayed away from it. By contrast, the multicanonical method is self-correcting because it reconnects the two phases at each successive simulation state point on the phase boundary.

Another widely used technique is the Gibbs Ensemble MC Method . Two separate simulation boxes (one for each pure phase) are connected thermodynamically by the dual means of particle transfers between the boxes and fluctuations in the box volumes implemented such that the overall volume of the two boxes remains constant. Like Gibbs-Duhem integration, GEMC is very useful for obtaining rough estimates of phase boundaries. The method is free of integration errors since the phases remain directly linked at each phase boundary state point. However, measurements of the pressure and chemical potential must be obtained by a sampling scheme since neither is imposed in the simulation. Moreover, direct comparisons have shown GEMC to be much less efficient than the methods described in this article because of the computational complexity associated with implementing volume fluctuations . Worries have also recently been voiced that the GEMC method suffers more severely from finite size effects than does the GCE .

Finally, it should be pointed out that despite their utility in dealing with phase equilibria involving fluids, the specific multicanonical techniques we have described do not permit one to tackle phase transitions involving solids. The problem here is that when attempting to traverse mixed-phase states involving crystalline order, the simulation invariably gets caught in ‘non-ergodic traps’ identifiable with defective crystalline configurations. Recently, however, a new technique has been developed that circumvents this problem by linking the two coexisting phases without traversing mixed-phase states. Essentially the method can be thought of as leaping directly from the configuration space of one pure phase to that of the other. Again use of multicanonical sampling is necessary, but this time its role is to encourage the simulation to visit a subset of configurations (in each pure phase) from which a leap to the other phase will be accepted. This new method has recently shown its worth in studies of solid-phase free energy differences and hard sphere freezing .

## VI Suggestions for further study

The following set of problems will help to reader to become more familiar with the techniques described in this article.

1. Using Box 1. as a guide, write a simple Grand Canonical ensemble MC program to simulate the Lennard-Jones fluid with potential cut-off at $`r_c=2.5\sigma `$. Do not correct for the potential truncation. Make sure the program prints out the particle number and the energy in list form (cf. sec. III C). Further programming details can be found in ref..
2. Run the program at the near critical phase boundary parameters  $`\beta ϵ=1.1876`$, $`\beta \mu =2.778`$, saving the output data list to a file.
3. Write a post-processing program to construct the number density histogram $`p(N)`$ from the raw data list.
4. Modify your GCE acceptance probabilities to cater for multicanonical preweighting (cf. secs. III A and III C). Use your measured $`p(N)`$ as the preweighting function for a multicanonical simulation at the same $`\beta ,\mu `$ used in (2) above. Hint: At the extrema of small (large) $`N`$, there will be histogram bins in $`p(N)`$ having zero entries. Before using $`p(N)`$ as your preweighting function, set these entries to be a constant equal to the smallest non zero entry. This avoids possible division by zero in the acceptance probabilities.
5. Extend your post-processing program to unfold the effects of the preweighting (as described in Secs. III A and III C) in order to find $`p(N)`$. Check that the form of $`p(N)`$ thus obtained agrees with that found without multicanonical preweighting. Compare the correlation time for the sampling processes with and without multicanonical preweighting.
6. Further extend your post-processing program to implement histogram reweighting (as described in Sec. III B and III C). Extrapolate the data obtained at the critical temperature $`\beta ϵ=1.1876`$ to find the location of the phase boundary and the form of $`p(N)`$ at $`\beta ϵ=1.17`$. Use this extrapolation as a preweighting function in a new multicanonical simulation at the new phase boundary state point. Compare your results with those of ref.