text:
# Direct correlation functions of the Widom-Rowlinson model

## I Introduction

Fluid binary mixtures may exhibit the phenomenon of phase separation. The simplest system able to undergo a demixing phase transition is the model introduced by Widom and Rowlinson some years ago Widom and Rowlinson (1970). Consider a binary mixture of non-additive hard spheres. This is a fluid made of hard spheres of specie 1 of diameter $`R_{11}`$ and number density $`\rho _1`$ and hard spheres of specie 2 of diameter $`R_{22}`$ and number density $`\rho _2`$, with a pair interaction potential between species $`i`$ and $`j`$ that can be written as follows

$`v_{ij}(r)=\{\begin{array}{cc}\mathrm{}\hfill & r<R_{ij}\hfill \\ 0\hfill & r>R_{ij}\hfill \end{array},`$ (3)
where $`R_{12}=(R_{11}+R_{22})/2+\alpha `$. The Widom-Rowlinson (WR) model is obtained choosing the diameters of the spheres equal to 0,

$`R_{11}=R_{22}=0,`$ (4)
so that there is no interaction between like spheres and there is a hard core repulsion of diameter $`\alpha `$ between unlike spheres. The symmetry of the system induces the symmetry of the unlike correlations \[$`h_{12}(r)=h_{21}(r),c_{12}(r)=c_{21}(r)`$\]. The WR model has been studied in the past by exact Ruelle (1971) and approximate Ahn and Lebowitz (1973, 1974); Bergmann (1976); Karkheck and Stell (1979) methods and it has been shown that it exhibits a phase transition at high density. More recently, additional studies have appeared and theoretical predictions have been confirmed by Monte Carlo (MC) computer simulations H. Luo, C. Hoheisel, and J. Karkheck (1976); Dickman and Stell (1995); Shew and Yethiraj (1996); G. Johnson, H. Gould, J. Machta, and L. K. Chayes (1997)

In this paper we will study the three dimensional symmetric Widom-Rowlinson mixture for which $`\rho _1=\rho _2=\rho /2`$, where $`\rho `$ is the total number density of the fluid, and

$`h_{11}(r)`$ $`=`$ $`h_{22}(r),`$ (5)
$`c_{11}(r)`$ $`=`$ $`c_{22}(r).`$ (6)
Moreover we know from (3) that the partial pair correlation function $`g_{ij}=h_{ij}+1`$ must obey

$`g_{ij}(r)=0\text{for }r<R_{ij}.`$ (7)
Our main goal is to focus on the direct correlation functions (dcf) of the WR model as a simplified prototype of non-additive hard spheres (NAHS) systems. The reasons to focus on the dcf’s is twofold: on the one hand, they are easier functions to model and fit. On the other hand, they play a central role in approximate theories like the Percus-Yevick approximation or mean spherical approximation (MSA) Hansen and McDonald (1986). We hope that a better understanding of the dcf’s properties in the WR model, could help in developing accurate analytical theories for the NAHS systems.

We calculate through Monte Carlo simulations the like $`g_{11}^{(MC)}(r)`$ and unlike $`g_{12}^{(MC)}(r)`$ pair distribution functions for a system large enough to allow a meaningful determination of the correspondent partial direct correlation functions $`c_{11}^{(MC)}(r)`$ and $`c_{12}^{(MC)}(r)`$, using the Ornstein-Zernike equation Hansen and McDonald (1986). We compare the results for the unlike direct correlation function with the results of the Percus-Yevick (PY) analytic solution found by Ahn and Lebowitz Ahn and Lebowitz (1973, 1974). In the same spirit as the work of Grundke and Henderson for a mixture of additive hard spheres Grundke and Henderson (1972), we propose a fit for the functions $`\mathrm{\Delta }_{11}^c(r)=c_{11}^{(MC)}(r)`$ and $`\mathrm{\Delta }_{12}^c(r)=c_{12}^{(MC)}(r)c_{12}^{(PY)}(r)`$.

At the end of the paper we also show the results from two Monte Carlo simulations on a mixture of non-additive hard spheres with equal diameter spheres $`R_{11}=R_{22}=R_{12}/2`$ and on one with different diameter spheres $`R_{11}=0`$ and $`R_{22}=R_{12}`$ to study the effect of non zero like diameters on the WR dcf’s.

## II Monte Carlo simulation and PY solution

The Monte Carlo simulation was performed with a standard NVT Metropolis algorithm Allen and Tildesley (1987) using $`N=4000`$ particles. Linked lists Allen and Tildesley (1987) have been used to reduce the computational cost. We generally used $`5.2\times 10^8`$ Monte Carlo steps where one step corresponds to the attempt to move a single particle. The typical CPU time for each density was around 20 hours (runs at higher densities took longer than runs at smaller densities) on a Compaq AlphaServer 4100 5/533.

We run the simulation of WR model at $`6`$ different densities $`\overline{\rho }=\rho \alpha ^3=0.28748`$, $`0.4`$, $`0.45`$, $`0.5`$, $`0.575`$, and $`0.65`$. Notice that the most recent computer simulation calculations Shew and Yethiraj (1996); G. Johnson, H. Gould, J. Machta, and L. K. Chayes (1997) give consistent estimates of the critical density around 0.75. Our data at the highest density (0.65) are consistent with a one phase system.

The Monte Carlo simulation returned the $`g_{ij}(r)`$ over a range not less than $`9.175\alpha `$ for the densest system. In all the studied cases the pair distribution functions attained their asymptotic value well inside the maximum distance they were evaluated. Thus, it has been possible to obtain accurate fourier transforms of the correlation functions \[$`h_{ij}(k)`$\]. To obtain the $`c_{ij}(r)`$ we used Ornstein-Zernike equation as follows

$`c_{11}(k)`$ $`=`$ $`{\displaystyle \frac{h_{11}(k)\left[1+\frac{\rho }{2}h_{11}(k)\right]\frac{\rho }{2}h_{12}^2(k)}{\left[1+\frac{\rho }{2}h_{11}(k)\right]^2\left[\frac{\rho }{2}h_{12}(k)\right]^2}}`$ (8)
$`c_{12}(k)`$ $`=`$ $`{\displaystyle \frac{h_{12}(k)}{\left[1+\frac{\rho }{2}h_{11}(k)\right]^2\left[\frac{\rho }{2}h_{12}(k)\right]^2}}`$ (9)
From the $`h_{ij}(k)`$ and $`c_{ij}(k)`$ we get the difference $`\gamma _{ij}(k)=h_{ij}(k)c_{ij}(k)`$ which is the fourier transform of a continuous function in real space. So it is safe to transform back in real space \[to get $`\gamma _{ij}(r)`$\]. Finally, the dcf’s are obtained from the differences $`h_{ij}(r)\gamma _{ij}(r)`$.

While for a system of non-additive hard sphere in three dimensions a closed form solution to the PY approximation is still lacking, Ahn and Lebowitz have found an analytic solution of this approximation for the WR model (in one and three dimensions).

The PY approximation consists of the assumption that $`c_{ij}(r)`$ does not extend beyond the range of the potential

$`c_{ij}(r)=0\text{for }r>R_{ij}.`$ (10)
Combining this with the exact relation (7) and using the Ornstein-Zernike equation we are left with a set of equations for $`c_{ij}(r)`$ and $`g_{ij}(r)`$ which have been solved analytically by Ahn and Lebowitz.

Their solution is parameterized by a parameter $`z_0`$. They introduce the following two functions of $`z_0`$ (which can be written in terms of elliptic integrals of the first and third kind)

$`I_1`$ $``$ $`{\displaystyle _{z_0}^{\mathrm{}}}{\displaystyle \frac{dz}{z\sqrt{z^3+4z/z_04}}},`$ (11)
$`I_2`$ $``$ $`{\displaystyle _{z_0}^{\mathrm{}}}{\displaystyle \frac{dz}{\sqrt{z^3+4z/z_04}}},`$ (12)
and define $`z_0`$ in terms of the partial densities $`\rho _1`$ and $`\rho _2`$ as follows

$`\eta 2\pi \sqrt{\rho _1\rho _2}={\displaystyle \frac{(I_2/2)^3}{\mathrm{cos}I_1}}.`$ (13)
They then define the following functions (note that in the last equality of equation (3.76) in Ahn and Lebowitz (1974) there is a misprint)

$`\overline{c}_{12}(k)`$ $``$ $`{\displaystyle \frac{2}{\sqrt{\rho _1\rho _2}}}\sqrt{{\displaystyle \frac{1+Y}{z_0^3Y^3+4Y+4}}}`$ (14)
$`\times \mathrm{sin}\left[{\displaystyle \frac{1}{2}}\sqrt{z_0^3Y^3+4Y+4}{\displaystyle _1^{\mathrm{}}}{\displaystyle \frac{dz}{(z+Y)\sqrt{z_0^3z^3+4z4}}}\right],`$
$`\overline{h}_{12}(k)`$ $``$ $`\overline{c}_{12}(k)[1\rho _1\rho _2\overline{c}_{12}(k)],`$ (15)
where $`Y(2k/I_2)^2`$.

We also realized that some other misprint should be present in the Ahn and Lebowitz paper since we have found empirically that the PY solution (with $`k`$ in units of $`\alpha `$) should be given by

$`c_{12}(k)=\overline{c}_{12}(ks),`$ (16)
where $`s`$ is a scale parameter to be determined as follows

$`s=[\overline{h}_{12}(r=0)]^{1/3}`$ (17)
Notice that for the symmetric case $`\rho _1=\rho _2=\rho /2`$ and $`\eta =\pi \rho =0.90316\mathrm{}`$ we find $`z_0=1`$ and $`s=1`$.

In Figs. 1, 2, and 3 we show three cases corresponding to the extreme and one intermediate density. In the figures, we compare the MC simulation data with the PY solution for the partial pair distribution functions and the partial direct correlation functions. Our results for the partial pair distribution functions at $`\rho \alpha ^3=0.65`$ are in good agreement with the ones of Shew and Yethiraj Shew and Yethiraj (1996). The figures show how the like correlation functions obtained from the PY approximation are the ones that differ most from the MC simulation data. The difference is more marked in a neighborhood of $`r=0`$ and becomes more pronounced as the density increases.

## III Fit of the data

From the simulations we found that $`c_{12}^{(MC)}(r)<8\times 10^3`$ for $`r>\alpha `$ at all the densities studied. This allows us to say that $`\mathrm{\Delta }_{12}^c(r)0`$ for $`r>\alpha `$. Moreover we found that both $`\mathrm{\Delta }_{12}^c(r)`$ for $`r<\alpha `$, and $`\mathrm{\Delta }_{11}^c(r)`$ are very well fitted by Gaussians

$`\mathrm{\Delta }_{11}^c(r)`$ $``$ $`b_{11}\mathrm{exp}[a_{11}(r+d_{11})^2]\text{for all }r>0,`$ (18)
$`\mathrm{\Delta }_{12}^c(r)`$ $``$ $`b_{12}\mathrm{exp}[a_{12}r^2]\text{for }0<r<\alpha ,`$ (19)
In Figs. 4 and 5 we show the behaviors of the parameters of the best fit (18) and (19), with density. In order to check the quality of the fit, we did not use the data at $`\overline{\rho }=0.45`$ in the determination of the parameters. The points for $`a_{12}`$ and $`b_{12}`$ are well fitted by a straight line or a parabola. As shown in Fig. 4 the best parabolae are

$`a_{12}(\overline{\rho })`$ $`=`$ $`0.839+0.096\overline{\rho }1.287\overline{\rho }^2,`$ (20)
$`b_{12}(\overline{\rho })`$ $`=`$ $`0.155+0.759\overline{\rho }0.159\overline{\rho }^2.`$ (21)
Fig. 5 shows how the parameters for $`\mathrm{\Delta }_{11}^c(r)`$ are much more scattered and hard to fit. The quartic polynomial going through the five points, for each coefficient, are

$`a_{11}(\overline{\rho })`$ $`=`$ $`55.25+504.8\overline{\rho }1659.\overline{\rho }^2+2364.\overline{\rho }^31236.\overline{\rho }^4,`$ (22)
$`b_{11}(\overline{\rho })`$ $`=`$ $`171.41556.\overline{\rho }+5166.\overline{\rho }^27421.\overline{\rho }^3+3906.\overline{\rho }^4,`$ (23)
$`d_{11}(\overline{\rho })`$ $`=`$ $`128.91144.\overline{\rho }+3747.\overline{\rho }^25328.\overline{\rho }^3+2782.\overline{\rho }^4,`$ (24)
The difficulty in finding a good fit for these parameters may be twofold: first we are fitting $`\mathrm{\Delta }_{11}^c(r)`$ with a three (instead of two) parameters curve and second the partial pair distribution functions obtained from the Monte Carlo simulation are less accurate in a neighborhood of the origin (due to the reduced statistics there). This inaccuracy is amplified in the process of finding the partial direct correlation functions. Such inaccuracy will not affect significantly $`\mathrm{\Delta }_{12}^c(r)`$ which has a derivative very close or equal to zero near the origin, but it will significantly affect $`\mathrm{\Delta }_{11}^c(r)`$ which is very steep near the origin.

In order to estimate the quality of the fit we have used the simulation data at $`\overline{\rho }=0.45`$. From Fig. 4 we can see how the parabolic fit is a very good one. In Fig. 5 the point at $`\overline{\rho }=0.45`$ gives an indication of the accuracy of the quartic fit. We have also compared the pair distribution and direct correlation functions obtained from the fit with those from MC: both the like and unlike distribution functions are well reproduced while there is a visible discrepancy in the dcf from the origin up to $`r=0.5\alpha `$. However we expect that moving on the high density or low density regions (where the quartic polynomial becomes more steep) the quality of the fit will get worst. In particular the predicted negative values for $`a_{11}`$, in those regions, are completely unphysical and the fit should not be used to extrapolate beyond the range $`0.28<\overline{\rho }<0.65`$.

## IV From WR to non additive hard spheres

In order to see how the structure, and in particular the dcf’s of the Widom-Rowlinson model change as one switches on the spheres diameters we have made two additional Monte Carlo simulations. In the first one we chose $`\rho _1=\rho _2=0.65/R_{12}^3`$ and $`R_{11}=R_{22}=R_{12}/2`$. The resulting partial pair distribution functions and partial direct correlation functions are shown in Fig. 6. From a comparison with Fig. 3 we see how in this case the switching on of the like diameters causes both $`c_{12}(r)`$ for $`r<R_{12}`$ and $`g_{12}(r)`$ for $`r>R_{12}`$ to approach $`r=R_{12}`$ with a slope close to zero.

In the second simulation we chose $`\rho _1=\rho _2=0.65/R_{12}^3`$ and $`R_{11}=0`$, $`R_{22}=R_{12}`$. The resulting partial pair distribution functions and partial direct correlation functions are shown in Fig. 7. From a comparison with Fig. 3 we see how in this case the switching on of the like diameters causes both $`g_{11}(0)`$ and $`c_{11}(0)`$ to increase, and $`c_{12}(r)`$ to lose the nearly zero slope at $`r=0`$. As in the previous case $`g_{12}(r)`$ for $`r>R_{12}`$ approaches $`r=R_{12}`$ with a slope close to zero. The like $`22`$ correlation functions for $`r>R_{12}`$ vary over a range comparable to the one over which vary the like $`11`$ correlation functions of the WR model.

For both these cases there is no analytic solution of the PY approximation available and a better understanding of the behavior of the direct correlation functions may help in finding approximate expressions Gazzillo and Pastore (1989).

## V Conclusions

In this paper we have evaluated the direct correlation functions of a Widom-Rowlinson mixture at different densities through Monte Carlo simulation and we have studied the possibility of fitting the difference between MC data and the PY dcf’s. We found a very good parameterization of $`c_{12}(r)`$ for $`r<\alpha `$ \[see equations (19) and (20)-(21)\] and a poorer one for $`c_{11}(r)`$ \[see equations (18) and (22)-(24)\]. The difficulty in this last case probably arises from the necessity of using three parameters \[instead of just two needed for parameterizing $`c_{12}(r)`$\], although it cannot be completely excluded some effect of the decreasing precision of the simulation data near the origin.

In the last part of the paper we have illustrated with additional Monte Carlo data the changes induced in the WR dcf’s by a finite size of the excluded volume of like correlations. These results are meant to provide a guide in the search of a manageable, simple analytical parameterization of the structure of mixtures of non additive hard spheres which is still not available although highly desirable.

LIST OF FIGURES

* Top panel: partial direct correlation functions obtained from the Monte Carlo simulation (points) with the $`c_{12}^{(PY)}(r)`$ obtained from the PY approximation (line) at a density $`\rho \alpha ^3=0.28748`$. Bottom panel: partial pair distribution functions obtained from the Monte Carlo simulation compared with the ones obtained from the PY approximation at the same density. The open circles and the dashed line: the like correlation functions. Closed circles and the continuous line: the unlike correlation functions.
* Same as in Fig.1 at a density $`\rho \alpha ^3=0.4`$.
* Same as in Fig.1 at a density $`\rho \alpha ^3=0.65`$.
* We plot, for five different values of the density, the parameters $`a_{12}`$ (diagonal crosses) and $`b_{12}`$ (starred crosses) of the best Gaussian fit (19) to $`\mathrm{\Delta }_{12}^c(r)`$ for $`r<\alpha `$, and fit them with parabolae (lines). The parameters at $`\rho \alpha ^3=0.45`$ where not used for the parabolic fit and give an indication of the quality of the fit.
* We plot, for five different values of the density, the parameters $`a_{11},b_{11}`$ and $`d_{11}`$ (stars) of the best Gaussian fit (18) to $`\mathrm{\Delta }_{11}^c(r)`$, and draw the quartic polynomial (lines) through them. The parameters at $`\rho \alpha ^3=0.45`$ where not used to determine the quartic polynomial and give an indication of the quality of the fit.
* Monte Carlo results at a density $`\rho =\rho _1=\rho _2=0.65/R_{12}^3`$ for the partial direct correlation function (on top) and the partial pair distribution function (below) of a mixture of non additive hard spheres with $`R_{11}=R_{22}=R_{12}/2`$. The open circles denote the like correlation functions. The closed circles denote the unlike correlation functions.
* Monte Carlo results at a density $`\rho =\rho _1=\rho _2=0.65/R_{12}^3`$ for the partial direct correlation function (on top) and the partial pair distribution function (below) of a mixture of non additive hard spheres with $`R_{11}=0`$ and $`R_{22}=R_{12}`$. The open circles denote the like $`11`$ correlation functions. The open triangles denote the like $`22`$ correlation functions. The closed circles denote the unlike correlation functions.


format: text

text:
# CHANDRA OBSERVATIONS OF SNR 1987A

## INTRODUCTION

With a known age, distance, and the progenitor type, supernova (SN) 1987A provides an unprecedented opportunity for the study of the early evolution of a supernova remnant (SNR). The optical inner ring, as observed by the Hubble Space Telescope (HST) is believed to be the dense circumstellar medium (CSM) produced by the stellar winds from the massive progenitor (Luo and McCray, 1991). Inside of this inner ring is an HII region generated by the UV radiation from the progenitor (Chevalier and Dwarkadas, 1995). The development of optically bright spots along the inner ring (Pun et al., 1997; Garnavich et al., 1997) was interpreted as emission from the radiative shock where the blast wave begins to strike inward protrusions of the dense inner ring (Michael et al., 2000). In X-rays, SN 1987A was detected with the Röntgen Satellite (ROSAT) as an unresolved source and the observed X-ray flux was linearly increasing (Hasinger et al., 1996). The high angular resolution images from the Chandra observations have revealed a shell-like structure of the X-ray remnant (Burrows et al. 2000). Monitoring observations with Chandra have shown a continuous development of the X-ray bright spots, which are in general spatially coincident with the optical and radio spots; a steady increase of the soft X-ray flux; and evidence of the radial expansion of the SNR (Park et al., 2002). The thermal X-ray spectrum was fitted with an electron temperature of $`kT`$ $``$ 2.5 keV. The broadening of the atomic emission lines as detected with the dispersed spectrum indicated a shock velocity of $``$3400 km s-1, which provided direct evidence for electron-ion non-equilibrium behind the shock front (Michael et al., 2002). We have recently performed two new observations with Chandra, which bring the total of monitoring observations of SNR 1987A up to six. We here report an update on the results from all six monitoring Chandra observations of SNR 1987A.

## X-RAY IMAGES

The broad-band Chandra/ACIS images from six monitoring observations of SNR 1987A are presented in Figure 1. The overlay contours are from HST H$`\alpha `$ data. The observation dates of the HST data are presented in the parentheses. The X-ray remnant is shell-like, as previously reported (Burrows et al., 2000; Park et al., 2002). Although the asymmetric surface brightness between the east and the west is persistent for all six images, the developements of new X-ray spots in the western half are evident in the latest images. These developing X-ray spots are spatially correlated with the optical spots. These “simultaneous” developments of the X-ray and optical spots support that the X-ray bright spots are where the blast wave is decelerated as it approaches the dense protrusions of the inner ring. We may thus expect the emergence of a complete X-ray ring in the near future as the shock front eventually strikes the main body of the entire inner ring. Park et al. (2002) have reported energy-dependent X-ray morphologies of SNR 1987A; i.e., the soft X-ray images (0.3 $``$ 1.2 keV) were correlated with the optical images and the hard band image (1.2 $``$ 8.0 keV) was consistent with the radio image. Figure 2 shows the development of X-ray spots in both soft and hard bands for the past 17-month period. These images confirm the previous findings that the soft X-ray spots are generally coincident with the optical spots while the hard X-ray spots are consistent with the radio bright lobes. These correlations of X-ray features with optical and radio images are consistent with the interpretation that the X-ray emission is from slow and fast shocks propagating in the HII region toward the dense inner ring.

The combined line profile as detected with the ACIS/HETG has indicated a shock front velocity of $``$3400 km s-1 (Michael et al., 2002), which was consistent with the SNR expansion rate ($``$3500 km s-1) as measured with the radio data (Gaensler et al., 2000) and with the blast wave velocity inferred by hydrodynamic models (Borkowski et al., 1997). We can also estimate the SNR expansion rate by directly measuring the average radius of the X-ray emission as a function of time. We fit the radial count distribution with a Gaussian and use the mean radius of the best-fit Gaussian as the radial extent of the X-ray shell. Fitting these data to a linear rate, we obtain an X-ray expansion rate of 4954$`\pm `$1057 km s-1 (Figure 3).

Fig. 1. Exposure-corrected deconvolved X-ray images of SNR 1987A overlaid with the HST contours.

Fig. 2. Subband images of SNR 1987A.

Fig. 3. Long-term variation of the mean radius of the X-ray count distribution. The solid line is the best-fit linear rate (4954$`\pm `$1057 km s-1).

## SPECTRUM AND LIGHTCURVE

The undispersed X-ray spectrum from SNR 1987A, as extracted from the latest observation taken on May 2002, is presented in Figure 4. The spectra from five other observations show similar features to those in Figure 4. The X-ray spectrum is thermal and shows broad emission line features for the elemental species O, Ne, Mg, Si, and S. The observed spectrum is fitted with a plane-parallel shock with an electron temperature of $`kT`$ $``$ 2.5 keV and an ionization timescale of $`n_et`$ $``$ 5 $`\times `$ 1010 cm-3 s. The fitted metal abundances are typically subsolar, which is consistent with the LMC ISM (Russell and Dopita, 1992) and ring abundances (Lundqvist and Fransson, 1996). The fitted electron temperature is consistent within uncertainties for all six observations. The volume emission measure ($`EM`$) has constantly increased over the past $``$30 months as the SNR has brightened. In the bottom panels of Figure 5, we present the $`EM`$ (left) and the derived soft/hard X-ray fluxes (right) for the overall spectrum between January 2000 and May 2002. We have also separately performed the spectral analysis between eastern and western halves and present the best-fit $`EM`$ for each half in the upper panels of Figure 5 (because of the limited photon statistics, only the latest four observations were feasible for this spectral analysis by halves). While the overall $`EM`$ has been steadily increasing, the latest data begin to show different $`EM`$ evolutions between the east and the west. This is an interesting observation and may indicate some complex density variations of the CSM along the inner ring. Follow-up observations will be necessary to verify this differential developments of $`EM`$ between the east and the west. We derive the 0.5 $``$ 2.0 keV flux from each observation (Table 2) in order to monitor the long-term variations of the soft X-ray flux. The latest lightcurves are presented in Figure 6. The X-ray fluxes in Table 2 have been corrected for the ACIS quantum efficiency (QE) degradation. For comparison, the derived soft X-ray fluxes before (triangles) and after (squares) the QE correction are presented in Figure 6b. We note that the effects of the QE degradation appear to become significant since the 2000 December observation. In Figure 6a, we present the radio flux variation, which shows a deviation from the linear increase rate since 1997 (around Day 4000). The solid lines in Figure 6b represent the linear rates fitted to the first four ACIS data points as presented in Park et al. (2002). It is evident that the ACIS rate increase has turned up and that the ACIS data can no longer be fitted with a linear rate. A quadratic increase rate for the combined (ROSAT \+ ACIS) lightcurve cannot fit the data (the dotted curve in Figure 6b). The X-ray flux has nearly tripled for the last $``$30 months.

Fig. 4. X-ray spectrum of SNR 1987A as observed on 2002-05-15.

## POINT SOURCE

As of 2002 May, we still observe no direct evidence of a point source within SNR 1987A. With the photon statistics of the latest data, a 90% upper limit of the embedded point source counts is 8% of the total SNR counts in the 2 $``$ 8 keV band. This implies an observed upper limit of $``$5.5 $`\times `$ 1033 ergs s-1 in the 2 $``$ 10 keV band.

Fig. 5. Emission measure/flux variations of SNR 1987A between 2000-01 and 2002-05.

Fig. 6. (a) The long-term lightcurve of SNR 1987A in the 4.7 GHz (ATCA). (b) The long-term X-ray lightcurve.

of six observations now, we detect newly developing X-ray spots. The soft X-ray spots are consistent with the optical spots while the hard X-ray spots are coincident with the radio lobes. These characteristics of the X-ray spots support the interpretation that the soft spots are the X-ray emission from the slow shock encountering dense protrusions of the inner ring and that the hard spots are emission from a fast forward shock propagating through the HII region. The emergence of new X-ray spots in the western side in addition to the eastern side indicates that the X-ray emission from SNR 1987A will be a complete ring in near future. With two new ACIS observations, we confirm the previously reported radial expansion rate of the X-ray SNR. The spectrum is fitted with a plane-parallel shock model with $`kT`$ $``$ 2.5 keV and low abundances, which is consistent with the previous results indicating X-ray emission primarily from the shocked CSM. The latest observation suggests a differential evolution of the volume emission measure between the east and the west, which might be an indication of complex density structure of the CSM along the inner ring. The current data are however insufficient to make a firm conclusion and follow-up observations should be helpful to understand this feature. We obtain a 90% upper limit on the 2 $``$ 10 keV band X-ray luminosity $`L_X`$ = 5.5 $`\times `$ 1033 ergs s-1 for any embedded point source.

## ACKNOWLEDGEMENTS

The authors thank P. Challis and the SINS collaboration for providing the HST data. We also thank B. Gaensler and the Australian Telescope Compact Array for providing the radio data. This work has been supported by NASA under the contract NAG8-01128 and SAO grants GO1-2064B and GO2-3098A.


text:
# 2-Dimensional Polymers Confined in a Strip

## I Introduction

The behaviour of flexible polymers in a good solvent confined to different geometries and in the presence of walls or other obstacles have been studied for many years deGen ; Eisenriegler93 . A particularly simple geometry is the space between two parallel walls. For simplicity we shall only discuss here the case of walls without energetic effects, i.e. the walls play a purely geometric role.

An important theoretical prediction is that near such a wall the monomer density profile increases as $`z^{1/\nu }`$, where $`z`$ is the distance from the wall ($`zD`$, and $`D`$ is the width between the two parallel walls) and $`\nu `$ is the Flory exponent deGen . This should hold in any dimension of space $`d`$. On the other hand it is intuitively obvious that the force exerted by the polymer onto the wall is proportional to the monomer density near the wall. The ratio between the two can be expressed in terms of a universal amplitude ratio which has been calculated by Eisenriegler Eisenriegler2  (using conformal invariance results of Cardy et al. cardy ) in $`d=2`$, and in any $`d4`$ by Eisenriegler Eisenriegler  by means of an $`ϵ`$-expansion. Attempts to verify these detailed predictions by Monte Carlo simulations in three dimensions Webman ; Ishinabe ; Milchev ; Joannis  have so far been without very convincing results. As far as we know, no attempt was made yet to verify them in $`d=2`$, and that is where the present paper sets in.

We study single polymer chains confined to a 2-d strip. They are described by self-avoiding random walks (SAWs) of $`N`$ steps on a square lattice between two hard walls with distance $`D`$ as shown in Fig. 1. More precisely, monomers are supposed to sit on lattice sites and $`D`$ is the number of rows accessible to monomers, i.e. the walls are placed at $`y=0`$ and at $`y=D+1`$, and the monomers can be at $`y=1,\mathrm{}D`$. We only consider the case where the Flory radius of a free chain of length $`N`$, $`R_FN^\nu `$, is much larger than $`D`$. When using a chain growth algorithm, the polymer has then to grow , after a short initial phase of $`D^{1/\nu }`$ steps, in either the positive or negative $`x`$-direction without possibility to change its orientation. This allows us to use an additional wall at $`x<x_0`$ which forces all chains to grow into the positive $`x`$-direction. For $`ND^{1/\nu }`$ this will essentially reduce the partition sum by a constant factor, without affecting any of the scaling laws or any of the detailed comparisons with theoretical predictions. On the other hand, it simplifies the subsequent discussion.

The force exerted onto the wall is most straightforwardly expressed in terms of the work done when moving one of the walls, i.e. by the dependence of the free energy – and thus also of the partition sum – on $`D`$,

$$F=k_BT\frac{\mathrm{ln}Z_N}{D},$$
(1)
where we have introduced a dummy temperature $`T`$ which can take any positive value. The partition sum $`Z_N`$ is just the number of $`N`$-step SAWs in the strip starting from a given $`x`$, but summed over all values $`1y_0D`$ of the $`y`$-component of the starting point.

The partition sum of a free SAW in infinite volume scales for $`N\mathrm{}`$ as $`Z_N=\mu _{\mathrm{}}^NN^{\gamma 1}`$ with $`\mu _{\mathrm{}}`$ being the critical fugacity per monomer, and with $`\gamma =43/32`$ being a universal exponent. In contrast, the partition function on a strip scales as

$$Z_N\mu _D^N,$$
(2)
without the power correction and with $`\mu _D`$ scaling for large $`D`$ as

$$\mu _D\mu _{\mathrm{}}aD^{1/\nu }.$$
(3)
with $`a`$ being a non-universal amplitude. The force per monomer is then obtained as

$$f=F/Nk_BTD^{11/\nu }.$$
(4)
Standard fixed-length Markov Chain Monte Carlo simulations do not give estimates of the partition sum or of the free energy, so that Eqs. (1-4) cannot be used directly. This has led to algorithms specifically designed for estimation of forces dickman , but employing the pruned-enriched-Rosenbluth method (PERM) g97  one can use Eqs. (1-4) directly.

Using PERM with $`k`$-step Markovian anticipation Frauenkron98 ; Frauenkron99 ; Caracciolo , we measured the partition sum directly and estimated the dependence of the monomer fugacity on the width $`D`$. In the same simulations also the monomer density profile, the end-to-end distance along the strip, and the density profile of chain ends are measured.

Details of the simulation method are given in the next section, results and their comparison with theoretical predictions are discussed in Sec. III.

## II Algorithm: PERM with $`k`$-step Markovian anticipation

PERM g97  is a chain growth algorithm with population control. Polymer chains are built like random walks by adding one monomer at each step. We use a Rosenbluth like bias for self-avoidance. As usual, this bias is compensated by a weight factor rosenbluth , i.e. each sample configuration should be given a weight. But actually we use a stronger bias which in addition suppresses dense configurations and samples more finely relatively open chain configurations, called Markovian anticipation in Frauenkron99 .

In $`k`$-step Markovian anticipation, we let the additional bias in the next step depend on the last $`k`$ steps made before. Let us denote the $`2d`$ directions on a $`d`$-dimensional hypercubic lattice by $`s=0,\mathrm{},2d1`$. All possible $`(k+1)`$-step configurations ($`k`$ previous steps $`i=1,2,\mathrm{}k`$ plus one future step $`i=0`$) are then indexed by

$$𝐒=(s_k,\mathrm{},s_0)=(𝐬,s_0).$$
(5)
In addition, choose an integer $`m100`$ (the precise value is not important). Either during an auxiliary run or during the early stages of the present run, we obtained a histogram $`H_m`$ such that $`H_m(𝐒)`$ is the accumulated weight at chain length $`n+m`$ of those chains which had configuration $`𝐒`$ during steps $`nk,nk+1,\mathrm{}n`$. Since we sample uniformly, $`H_0(𝐒)`$ is independent of $`𝐒`$. Thus $`H_m(𝐒)`$ with $`m>0`$ indicates how “successful” is a configuration $`𝐒`$ after $`m`$ more steps. In importance sampling we want each chosen direction to have in average the same later success. Therefore we choose the next step with probability

$$p(s_0|𝐬)=\frac{H_m(𝐬,s_0)}{_{s_0^{}=0}^{2d1}H_m(𝐬,s_0^{})},$$
(6)
In our simulations we choose $`k=9`$ and $`m=100`$. We accumulate contributions to $`H`$ only for $`n+m>300`$, and we apply Eq. (6) only for chain lengths $`>k`$ (for chain lengths $`<k`$ there is not yet enough history to condition upon).

Accumulating the histogram only for $`n1`$ is suggested by the fact that only for large $`n`$ the anisotropic bias is fully developed Frauenkron99 . This anisotropy makes the histogram strongly dependent on $`D`$. We found that using for all $`D`$ only the histogram obtained for free chains, i.e. for $`D=\mathrm{}`$, gives nearly the same efficiency. This is quite different from the case where the anisotropy is not due to geometry, but is due to stretching of the polymer. In the latter case, simulations with markovian anticipation become much more efficient with increased stretching unpub . This is not the case for the present problem where the anisotropy is due to geometric constraints, for which PERM with Markovian anticipation is however still the most efficient known simulation method by far.

## III Results

Before presenting our results, let us stress that we have several possibilities for checking our algorithm. For very large $`D`$ we can compare our estimates of $`\mu `$ with the very precise estimate $`\mu _{\mathrm{}}=0.37905228`$ Guttmann88 ; Guttmann96 . For $`D11`$ we can compare with exact transfer matrix calculations of Burkhardt . And for $`D2`$ we can even solve the problem analytically.

For $`D=1`$, the polymer can only grow in a straight configuration, giving $`\mu _{(D=1)}=1`$. For $`D=2`$, each step can be either up (u), down (d), or to the right (r). After an ‘u’ or ‘d’ move, the next step has to be ‘r’, while any move is possible after ‘r’. Lumping moves ‘u’ and ‘d’ together into a vertical move (‘v’), we see that the set of all possible configurations forms a regular language hopcroft  with the associated graph shown in Fig. 2. The partition sum for chains of length $`N`$ is just twice the $`(N+1)`$-st Fibonacci number,

$$Z_N(D=2)=2F_{N+1}(N1),$$
(7)
where $`F_0=F_1=1`$ and $`F_N=F_{N1}+F_{N2}`$. Thus the critical monomer fugacity for $`D=2`$ is the inverse of the golden mean,

$$\mu _{D=2}=g^1\frac{\sqrt{5}1}{2}=0.61803\mathrm{}.$$
(8)
In this case, we can also show that markovian anticipation gives the optimal bias. Markovian anticipation corresponds in this case to $`p_r:p_v=g:1`$ if a vertical move is allowed, and $`p_r:p_v=1:0`$ else. This choice leads to weights which oscillate between two values, thus no population control (pruning/cloning) is needed footnote . All this is verified in our simulations, which serves thus as a test for our algorithms.

For $`D>2`$, we used simulations. We simulated strip widths up to $`D=320`$ and chain length between $`N=3000`$ (for $`D=2`$) and $`N=125,000`$ (for $`D=320`$). Critical fugacities are determined by plotting $`Z_N\mu _D^N`$ against $`\mathrm{log}N`$ and demanding that these curves become horizontal for large $`N`$. Results are shown in Fig. 3, where we plot $`\mu _D\mu _{\mathrm{}}`$ with $`\mu _{\mathrm{}}=0.37905228`$ Guttmann96 . They are in perfect agreement with the theoretical prediction of Eq. (3), and provide the estimate $`a=0.7365\pm 0.0007`$. In addition, $`\mu _D`$ can be compared for $`D=3`$ to $`12`$ with the transfer matrix results of Burkhardt . For all $`D`$, the values agree for at least six digits.

As we said in the introduction, we expect the force $`f`$ onto the wall to be proportional to the monomer density $`\rho (y)`$ near the wall. The precise relation is given by Eisenriegler Eisenriegler

$$\underset{y0}{lim}k\frac{\rho (y)}{y^{1/\nu }}=B\frac{f}{k_BT}.$$
(9)
Here the non-universal amplitude $`k`$ relates the end-to-end distance of a free SAW to the chain length, $`k=R_x^{1/\nu }/N=0.5297\pm 0.0002`$ LiMadrasSokal  for the square lattice. On the other hand, $`B`$ is a universal number. For ideal chains $`B=2`$, while for chains with excluded volume in $`4ϵ`$ dimensions one has $`B2(1b_1ϵ)`$ with $`b_1=0.075`$ Eisenriegler . The latter is of course of dubious value for $`d=2`$ (where it would give $`B1.7`$), but conformal invariance leads to the supposedly exact value $`B=2.01`$ in $`d=2`$ cardy ; Eisenriegler2 .

The monomer density is normalized such that $`_{y=1}^D\rho (y)=1`$. According to Eq. (9) it should scale as $`y^{1/\nu }`$ near the walls (this holds in any dimension) with $`\nu =3/4`$. Surprisingly, we found that the simplest ansatz generalizing this power law to all $`y[1,D]`$,

$$\rho (y)=\frac{1}{D+1}f_\rho (\xi )\frac{1}{D+1}A(\xi (1\xi ))^{4/3},\xi =\frac{y}{D+1},$$
(10)
with $`A=\mathrm{\Gamma }(14/3)/\mathrm{\Gamma }^2(7/3)=10.38`$ gives already an excellent fit (Fig. 4a) in the limit $`D\mathrm{}`$. There are small deviations invisible in Fig. 4a but clearly seen when plotting $`\rho (y)/f_\rho (\xi )`$ (see Fig. 4b), but they seem to vanish slowly in the limit $`D\mathrm{}`$ (see Fig. 4b).

Assuming the latter, i.e. assuming that Eq. (10) becomes exact asymptotically, the universal number $`B`$ could be estimated by using Eqs. (1), (2), (3), and (9) which together would give

$$A=\underset{D\mathrm{}}{lim}\underset{y0}{lim}\frac{D^{7/3}\rho (y)}{y^{4/3}}=\frac{4}{3}\frac{Ba}{k\mu _{\mathrm{}}}.$$
(11)
Inserting the above numbers gives $`B=2.122\pm 0.002`$ which is definitely larger than the value predicted in cardy ; Eisenriegler2 , by some 50 standard deviations.

In an alternative scenario we could assume that the deviations from Eq. (10) seen in Fig. 4b do not vanish in the limit $`D\mathrm{}`$. In that case one should also allow for a modified scaling variable

$$\xi _\delta =(y\delta )/(D+12\delta )$$
(12)
with an unknown small (non-universal) parameter $`\delta `$. Values of $`(D+12\delta )\rho (y)/f_\rho (\xi _\delta )`$ for $`D=128`$ and three different values of $`\delta `$ are plotted in Fig. 4c. Similar results are obtained for other $`D`$. Combining them, we see that best scaling (i.e. least dependence of $`D`$ for small $`y`$) is obtained for $`\delta 0.02`$. For this value of $`\delta `$ we have

$$\underset{y0,D\mathrm{}}{lim}D^{7/3}\rho (y)/y^{4/3}=(0.95\pm 0.02)\times A.$$
(13)
The large uncertainty in this estimate reflects the rather steep slope of the central curve in Fig. 4c as $`\xi _\delta 0`$ and the related uncertainty in the best estimate of $`\delta `$ footnote2 . If we accept Eq. (13), we obtain $`B=2.04\pm 0.04`$, which is in perfect agreement with cardy ; Eisenriegler2 . Thus we have two scenarios. Both imply very large corrections to scaling. While the a priori simpler scenario would be in conflict with the theoretical prediction, this prediction suggests that indeed the second scenario is correct, which is our preferred solution.

In Fig. 5, we plot the end-to-end distance per monomer $`x/N`$ versus $`N`$ for various widths $`D`$. These curves become horizontal as $`N\mathrm{}`$, i.e. $`x`$ increases indeed linearly with $`N`$, $`lim_N\mathrm{}x/N=\mathrm{\Delta }x`$. In order to find how $`\mathrm{\Delta }x`$ scales with $`D`$, we plot it in Fig. 6 on a doubly logarithmic scale. As indicated by the dashed line, it is fitted perfectly by the theoretical prediction Eisenriegler

$$x/ND^{11/\nu }=D^{1/3}.$$
(14)
We also estimated the density of wall contacts $`\rho _b`$ (number of monomers at $`y=1`$ or at $`y=D`$, divided by $`2x`$). For each fixed value of $`D`$, this becomes independent of $`N`$ as $`N\mathrm{}`$. The asymptotic values, obtained from plots similar to Fig. 5, are also plotted in Fig. 6. The full line corresponds to the very simple prediction Eisenriegler

$$\rho _bD^2$$
(15)
which is independent of $`\nu `$ and indeed holds also for Gaussian chains. Eq. (15) can be easily understood, in terms of the pressure exerted onto the wall:

$$\rho _bp=\frac{Nf}{x}D^{11/\nu +1/3}=D^2.$$
(16)
Finally, we show in Fig. 7 the distribution $`\rho _{end}(y)`$ of chain ends. We see that $`\rho _{end}(y)`$ scales for large $`D`$ (and for $`ND^{1/\nu }`$, of course). Theoretically it is predicted that Cardy-Redner

$$\rho _{end}(y)y^{25/48}$$
(17)
for $`yD`$. We see that this is indeed verified (the full lines in Fig. 7 correspond to the prediction). But in contrast to the monomer density which was described fairly well for all $`y`$ by the product of the power laws holding near the two walls, the same definitely does not hold for $`\rho _{end}(y)`$. There the function $`g_\rho (\xi )=const(\xi (1\xi ))^{25/48}`$ does not describe the behaviour away from the walls.

## IV Summary

We have shown that we could simulate 2-d polymers, modelled as self-avoiding walks, with chain length up to 125000 on strips of widths up to 320. This was possible using the PERM algorithm with Markovian anticipation. The fact that PERM gives by default very precise estimates of free energies allowed us to measure precisely the forces exerted onto the walls, by measuring how the critical fugacities depend on the width of the strips. We verified all critical scaling laws predicted for this problem, including the scaling of monomer and end point densities near the walls and the scaling of the total pressure with chain length and with strip width.

The only prediction for which we found possibly disagreement is for the universal amplitude ratio $`B`$ defined in Eq. (9). A scenario based on some minimal assumption about scaling functions and corrections to scaling gives an estimate higher than the prediction by some fifty standard deviations. But a different scenario, maybe less plausible a priori but not very unlikely either, gives perfect agreement with the prediction. This illustrates again that one should be very careful about corrections to scaling, and that even very precise simulations do not always give unique answers when their analysis is not guided by a reliable theory.

Previous simulations of 3-d polymers between two parallel planar walls had indicated that also there the value of $`B`$ might be larger than predicted, but those simulations had very large uncertainties. Using PERM we can simulate fairly easily much longer chains with rather high statistics. Results of such 3-d simulations will be given elsewhere.

Acknowledgements: We thank Profs. Erich Eisenriegler and Ted Burkhardt for valuable discussions, and Walter Nadler for carefully reading the manuscript.