text:
# The Spectacular BHR 71 Outflow

## 1 Introduction

In 1977 Arge Sandqvist published a catalogue of southern “dark dust clouds of high visual opacity” (Sanqvist 1977 – 95 entries, numbered 101-195), an extention of an earlier paper with Lindroos (Sandqvist & Lindroos 1976) in which they presented H2CO absorption line studies of 42 dust clouds (#1-42). Number 136 on Sandqvist’s list (Sa 136) is a very opaque Bok globule located near the Coalsack, later catalogued as DC 297.7-2.8 by Hartley et al. (1986) and as entry 71 in the globule list of Bourke, Hyland & Robinson (1995a – BHR 71). Mark McCaughrean in his opening address at this conference highlighted a number of important events that occurred in 1977, in particular IAU Symposium 75 on Star Formation whose proceedings appeared that year. It is fitting that BHR 71, which is featured in a beautiful VLT optical image in the frontpiece (& poster) of these proceedings, can trace its origins in the literature to that same year.

## 2 Globule properties

The globule properties have been determined by Bourke et al. (1995b, 1997). Spatially and kinematically BHR 71 is associated with the Coalsack at an assumed distance of 200 pc, though it may be as close as 150$`\pm `$30 pc (Corradi et al. 1997). Large scale 12CO & 13CO maps of the globule give a size of $``$0.5 pc and mass 40$`M_{}`$, while C18O observations which trace high column density gas imply a size $`0.3\times 0.15`$ pc and a mass of 12$`M_{}`$. The high density (n$`>10^4`$ cm-3) core traced in ammonia is $`0.2\times 0.1`$ pc in size with a mass of 3$`M_{}`$. The globule velocity is $`V_{lsr}`$ $``$ $`4.5`$km s-1.

## 3 CO Outflow Properties

The properties of the large scale molecular outflow have been determined by B97. As can be seen in Fig. 1, the outflow lobes are well separated on the sky, and extend $``$0.3pc from their origin with an opening angle of $``$15. B97 find that the velocity structure is consistent with a steady flow with constant velocity (Cabrit et al. 1988). With this assumption the inclination of the outflow to the line-of-sight is determined to be 85. The CO excitation temperature in the line core is greatest at the outflow peaks, indicating that the ambient gas there has been heated by interactions with the outflow.

Correcting for inclination, optical depth, and emission hidden within the line core, B97 determine the mass in the lobes to be $``$1.0$`M_{}`$ (red lobe) and $``$0.3$`M_{}`$ (blue lobe). Considering the different methods used to determine the outflow momentum $`P`$, kinetic energy $`E_k`$, and mechanical luminosity $`L_{\mathrm{mech}}`$ (upper and lower limit methods) B97 find $`P=11`$$`M_{}`$km s-1, $`E_k=60`$$`M_{}`$km2s-2, and $`L_{\mathrm{mech}}=0.5`$$`L_{}`$. There is less mass in the blue lobe, which may be a result of this lobe breaking out of the globule, indicated in Fig. 1 by the conical reflection nebulosity just south of the protostars (see the beautiful colour VLT in the frontpiece of these proceedings for a more detailed view).

## 4 Two protostars - two outflows

Near-infrared (NIR) images from the AAT are shown in Figure 2 (Bourke 2001). Most of the non-stellar emission is due to the emission in the H2 $`v`$=1-0S(1) line, most likely due to shocks in the outflowing gas (Eislöffel 1997). The NIR emission is well aligned with the large scale CO outflow (Fig. 3).

Mid-infrared (MIR) emission in the ISO LW2 band is overlayed on Fig. 2(a). Two of the 7$`\mu `$m sources appear to be located at the apexes of NIR emission, strongly suggesting that they are associated with the emission. Source “1” (hereafter IRS 1) lies at the apex of the reflection nebulosity seen also in Fig. 1 and is co-incident with the position of the mm source BHR 71-mm, also known as IRAS 11590-6452 (B97). The 7$`\mu `$m flux from IRS 1 is an order-of-magnitude greater than from IRS 2. The NIR feature coincident with IRS 2 in Fig. 2(a) is non-stellar, by comparison of its PSF with stars in the same image.

A cm continuum source (indicated on Fig. 2(b)) is detected toward BHR 71 IRS 1, at both 3 and 6 cm (Wilner et al. 2001, in prep). The spectral index is consistent with a flat or rising spectrum due to free-free emission, a signpost of protostellar origin (Rodríguez 1994). Corporon & Reipurth (1997) discovered two Herbig-Haro associations in BHR 71 – HH 320 and HH 321, and their locations are shown on Fig. 2(b). It can be seen that HH 320 (HH 321) is coincident with the NIR emission associated with IRS 2 (IRS 1).

Bourke (2001) has shown that IRS 2 also drives a CO outflow which is more compact and much less energetic than the IRS 1 outflow. The northern part of the IRS 2 outflow is blue-shifted (and associated with HH 320) which is the opposite of the IRS 1 outflow and allows them to be separated spatially. The red lobe is confused by the IRS 1 outflow, though it is probably seen in the NIR (arrowed emission in Fig. 2(b)). Bourke (2001) suggested that IRS 1 & 2 may form a binary protostellar pair (separation $``$3400AU) though the kinematic evidence for or against is lacking.

## 5 Outflow Chemistry

The BHR 71 IRS 1 outflow is one of a handful that show significant abundance enhancements in molecules such as SiO and CH3OH (G98). Figure 3 shows the spatial distribution of SiO and CO in the outflow, compared to the NIR H2 emission (the CO data is of lower spatial sampling than Fig. 1). Figure 4 shows spectra at two locations, the red lobe (as indicated by the box in Fig. 1) and at the position of IRS 1.

The spectral line profiles in the outflow and the velocity of the outflowing gas ($`<30`$km s-1) indicate that C-shocks dominate the flow (G98). The shocks are sufficiently strong to release molecules and atoms into the gas phase via evaporation of icy grain mantles (e.g., CH3OH) and sputtering of grain cores or grain-grain collisions (e.g., Si, which rapidly forms SiO). Other molecules detected in the outflow include CS, H2CO, SO, HCN, HNC, HCO+ with SEST and H2O with SWAS (Bourke et al. 2001, in prep). SiO is removed from the gas phase in about 104 years indicating that the outflow is quite young.

In the red lobe G98 determined abundance enhancements of $``$350 in SiO and $``$40 in CH3OH. One particularly striking feature is the spatial distribution of SiO compared to CO in the outflow. Because SiO is the result of Si liberation it is usually only detected at the ends of outflows (where the shock interaction is greatest) or as a narrow jet along the outflow axis possibly due to interactions in a turbulent boundary layer (Garay 2000). The wide-spread distribution of SiO in the BHR 71 outflow is unique. This suggests that the SiO enhancement takes place in a shell-like structure produced by the dynamical interaction between the ambient cloud and an underlying wide-angle wind or wind driven shell (Garay 2000), or perhaps by a wandering jet. However, it has not been shown that an interaction between a wind and the ambient material can produce sufficient Si for this to be a viable explanation. If sufficient Si-bearing species are present in grain mantles then the wind model becomes attractive (Schilke et al. 1997).

## 6 The VLT Image

An optical composite image taken with the VLT is shown in the frontpiece of these proceedings. This image hints at the spectacular results we can expect from the VLT in the coming years. There is evidence in this image of both a wind component and a jet component to the IRS 1 outflow in the blue lobe. Extending from the reflection nebulosity which protrudes from the globule, one can trace out an elongated bubble, with its edges defined by enhanced extinction. This is characteristic of a wide-angle wind component. In addition, enhanced extinction is also seen along the axis of the bubble and extending beyond its southern tip. This may be an indication of the underlying jet which is probably driving this young outflow. Modelling of this one image may help answer some of the remaining questions about the spectacular BHR 71 outflow.

I thank my many collaborators on this project, in particular Guido Garay. A big hug to João for letting me present my work on this beautiful object.


format: text

text:
# 1 Introduction: The CORALS Survey

## 1 Introduction: The CORALS Survey

It is widely recognised that DLAs represent one of the most promising methods of probing the high redshift galaxy population. However, it is similarly acknowledged that the presence of dust in these galaxies could impose a significant bias on our current database of DLAs, which are mostly drawn from optically bright QSO samples. Observational evidence for such a concern comes, for example, from statistically steeper continuum slopes in QSOs with DLAs  and lack of metallicity evolution in DLAs from $`0.5<z_{abs}<3.5`$ .

The CORALS survey aims to address this concern by using a radio-selected sample of QSOs to search for DLAs. Full details of this sample can be found in . In brief, a complete ($`S_{2.7GHz}>0.25`$ Jy) sample of $`z_{em}2.2`$ QSOs was culled from the Parkes Catalogue. Optical counterparts have been identified for every target, resulting in a sample of 66 radio-selected QSOs with no optical magnitude limit. Follow-up spectroscopy revealed 22 $`1.8<z_{abs}<z_{em}`$ DLAs, of which 3 are excluded from our statistics due to their proximity to the QSO ($`\mathrm{\Delta }v<`$ 3000 km/s).

## 2 Results

In order to assess the extent to which previous surveys may have been biased due to dust, we calculate the number density, $`n(z)`$, and neutral gas mass density, $`\mathrm{\Omega }_{DLA}`$, of CORALS DLAs. We find that $`\mathrm{log}\mathrm{\Omega }_{\mathrm{DLA}}h_{65}=2.59_{0.24}^{+0.17}`$ ($`\mathrm{\Omega }_M=1.0,\mathrm{\Omega }_\mathrm{\Lambda }=0`$) over the redshift range $`1.8<z_{abs}<3.5`$, in good agreement with the latest results from , see Figure 1. Perhaps most importantly, we do not find any evidence for a previously undetected population of high column density absorbers. Similarly, the number density of CORALS DLAs $`n(z)=0.31_{0.08}^{+0.09}`$ at a mean absorption redshift $`z_{abs}=2.37`$ is 50% larger than previous surveys (e.g. ). However, the difference is significant only at around the 1$`\sigma `$ level. Overall, we conclude that previous surveys may have underestimated $`n(z)`$ and $`\mathrm{\Omega }_{DLA}`$ by at most a factor of two.


text:
# Numerical simulation of a binary communication channel: Comparison between a replica calculation and an exact solution

## 1 Introduction

The extraction of sensory information by the brain from a stream of multi-dimensional data may be understood as a process of optimisation of mutual information $`(MI)`$  or of redundancy . The $`MI`$ measures the statistical dependence between two random variables . In our case, they correspond to an $`N`$-dimensional input signal $`\stackrel{}{\xi }`$ provided by the sensory receptors, and a $`P`$-dimensional output $`\stackrel{}{v}`$. More precisely, the $`MI`$ indicates the amount of knowledge about $`\stackrel{}{\xi }`$ that can be extracted from $`\stackrel{}{v}`$ (See e. g. ref.). With respective probability densities $`p_\xi `$ and $`p_v`$, this is given by

$$I[p_v,p_\xi ]=<\mathrm{log}\frac{p_{v,\xi }(\stackrel{}{v},\stackrel{}{\xi })}{p_v(\stackrel{}{v})p_\xi (\stackrel{}{\xi })}>_{v,\xi },$$
(1)
here $`\mathrm{log}(x)\frac{\mathrm{ln}(x)}{\mathrm{ln}2}`$. The average $`<\mathrm{}>_{v,\xi }`$ is over the joint probability distribution $`p_{v,\xi }(\stackrel{}{v},\stackrel{}{\xi })`$. For instance, if $`\stackrel{}{v}`$ and $`\stackrel{}{\xi }`$ are independent, we have $`p_{v,\xi }=p_vp_\xi `$ and so $`I=0`$.

The problem of learning the statistical properties of a set of $`N`$-dimensional correlated Gaussian inputs is well understood for a linear channel, even in the presence of noise (which is actually necessary to regularise $`MI`$)-. A non-linear continuous channel has also been studied in the low-noise limit for rather general transfer functions . It was shown that maximisation of the $`MI`$ leads to a factorial code. Threshold-linear networks  (treated with the replica technique) have also been considered. On the other hand, the binary channel, where the outputs take discrete values (say $`v_\mu =\pm 1`$), is not well understood, although the problem has been studied using replica-symmetric ($`RS`$) statistical-mechanical techniques -. Most interestingly, an analytical solution has been found, and the existence of a large order phase transition as the number of output neurons increases has been suggested . The relevant parameter to describe this transition’s occurrence is the ratio between the number of output and of input units, $`\alpha =P/N`$. The analytical solution holds up to a value of $`\alpha `$ of order one, beyond which it is not longer correct. On the contrary, the $`RS`$ solution does not exhibit any transition and gives good approximations at both the small and large $`\alpha `$ regimes. It has been proved that below some $`\alpha 1`$ the analytical and the $`RS`$ solutions are very close. In fact, an expansion in powers of $`\alpha `$ shows that the two solutions are identical up to $`O(\alpha ^2)`$. In spite of the fact that from the third order the corresponding expansions differ, the numerical agreement up to $`\alpha 1`$ is excellent (a relative difference of less than 0.9% up to $`\alpha =0.1`$). This is due to intriguing cancellations between higher orders.

Here we present simulations with the aim of providing numerical evidence on the validity of the analytical solution. In addition, since the order of the transition is large and the $`RS`$ solution seems to interpolate well between the small $`\alpha `$ and the asymptotic regimes, we also compare numerical simulations at several values of $`\alpha `$ with the replica theory prediction .

## 2 The Binary Channel

We consider a single-layer perceptron, or $`channel`$, with $`N`$ continuous input neurons whose states $`\xi _i`$ define a vector $`\stackrel{}{\xi }=\{\xi _i\}_{i=1}^N`$ representing the $`signal`$ received from the environment. The output layer has $`P`$ binary neurons the values $`v_\mu =\pm 1`$ of which compose the vector $`\stackrel{}{v}=\{v_\mu \}_{\mu =1}^P`$, that represents the $`code`$. Between the signal and the code there is an encoder, given by a set of synaptic couplings $`𝒥\{𝒥_{i\mu }\}`$.

The inputs take values drawn from an $`N`$-dimensional Gaussian probability distribution, unbiased ($`<\xi _j>=0`$) and with correlation matrix $`𝐂_X^{ij}=<\xi _i\xi _j>`$,

$$p_\xi (\stackrel{}{\xi })=e^{\stackrel{}{\xi }𝐂_X^1\stackrel{}{\xi }/2}/\sqrt{2\pi C_X},C_X=|𝐂_X|,$$
(2)
which using a convenient shorthand can be expressed as $`\stackrel{}{\xi }N(\stackrel{}{0},𝐂_X)`$. (Here $`X`$ is a correlation parameter between input neurons, to be defined more precisely later). The transfer function is deterministic, so that $`v_\mu =\text{sign}(h_\mu )`$, where

$$h_\mu =\underset{i=1}{\overset{N}{}}J_{i\mu }\xi _i\stackrel{}{J}_\mu \stackrel{}{\xi }$$
(3)
and the $`\stackrel{}{J}_\mu `$’s denote the synaptic weight vectors linking the signal $`\stackrel{}{\xi }`$ to each output neuron $`\mu `$. They form a set of independent random vectors $`\{\stackrel{}{J}_\mu \}_{\mu =1}^P`$, each distributed according to an $`N`$-dimensional Gaussian probability, with mean $`<J_{i\mu }>=0`$ and correlation matrix $`𝚪_\mu `$, whose elements are $`\mathrm{\Gamma }_\mu ^{ij}=<J_{i\mu }J_{j\mu }>`$. This means that $`\stackrel{}{J}_\mu N(\stackrel{}{0},𝚪_\mu )`$.

We are mainly interested in computing the averaged mutual information per input unit in the thermodynamic limit, i.e. $`ilim_N\mathrm{}\frac{1}{N}<I[p_v,p_\xi ]>_𝒥`$. A useful result to have in mind is that the $`MI`$ in eq.(1) is the difference

$$I[p_v,p_\xi ]_𝒥=H[p_v]H[p_{v|\xi }],$$
(4)
where $`H[p_v]=<\mathrm{log}p_v>_v`$ is the entropy of the output, while $`H[p_{v|\xi }]=<<\mathrm{log}p_{v|\xi }>_{v|\xi }>_\xi `$ is the conditional entropy of the output given the input (averaged over the input). The code $`\stackrel{}{v}`$, given a fixed signal $`\stackrel{}{\xi }`$, has the conditional probability $`p_{v|\xi }=p_{v,\xi }/p_\xi `$. Since a deterministic channel clearly has zero conditional entropy, in that case the $`MI`$ reduces to the output entropy.

## 3 Analytical results for an example

As an example we consider a Gaussian input distribution with two-point correlation $`𝐂_X`$ given by $`C_X^{ij}=X^{|ij|}`$. The input neurons are then less (more) spatially correlated if $`X1`$ ($`X1`$). For $`X=1`$, $`𝐂_X`$ is a singular matrix, while for $`X0`$ $`𝐂_X`$ tends towards the identity matrix. The correlations between two synapses converging to the same output, $`\mathrm{\Gamma }_\mu ^{ij}`$, are chosen to be equal for all output neurons and normalised to $`1`$, $`i.e.`$ $`\mathrm{\Gamma }_\mu ^{ij}=\mathrm{\Gamma }^{ij}=\delta _{ij},\mu `$.

For small $`\alpha `$, the MI can be expanded as $`i(\alpha )a_0+a_1\alpha +a_2\alpha ^2+a_3\alpha ^3`$. One trivially obtains that $`a_0=0`$ and $`a_1=1`$. Both the $`RS`$ solution (see ) and the analytical solution (see ) give the same value of $`a_2`$,

$$a_2(X)=\frac{\gamma }{\pi ^2\mathrm{ln}(2)}$$
(5)
where $`\gamma \frac{1+X^2}{1X^2}`$. On the other hand, the two techniques disagree in their predictions for $`a_3`$. The $`RS`$ method yields

$$a_3^{RS}(X)=\frac{(612/\pi )\gamma ^22}{3\pi ^3\mathrm{ln}(2)},$$
(6)
while (following the methods of ) one can easily verify that the analytical technique gives:

$$a_3^{An}(X)=\frac{6\gamma ^22}{3\pi ^3\mathrm{ln}(2)}.$$
(7)
The $`RS`$ solution also gives rise to predictions for strongly correlated inputs. For $`X1`$ one obtains:

$$iK(1X)^{1/3}\alpha ^\nu ,\nu 2/3.$$
(8)
Finally, the same solution shows logarithmic behaviour for large $`\alpha `$:

$$i\mathrm{log}\alpha +\frac{1}{2\mathrm{ln}2}+\mathrm{log}0.72+\frac{1}{2N}\text{Tr}[\mathrm{log}\frac{𝐂_𝐗}{N}].$$
(9)
## 4 Method

The first step of the simulation entails choosing a coupling sample $`𝒥`$ at random. A given choice of $`𝒥`$ will be labelled as the sample $`s`$, and the total number of samples will be denoted by $`S`$. Next, signals $`\stackrel{}{\xi }`$ are drawn in order to obtain several codewords $`\stackrel{}{v}`$. A histogram $`p_v^s(\stackrel{}{v})`$ is then constructed from these states, representing the probability $`p_v(\stackrel{}{v})`$ and allowing us to estimate the output entropy:

$$i_s=\frac{1}{N}\underset{\{v_\mu =\pm 1\}}{}p_v^s\mathrm{log}p_v^s.$$
(10)
Clearly, this approximation to the true code entropy $`H[p_v]`$ will improve as the number of drawn signals increases. In practice, we evaluated $`p_v^s(\stackrel{}{v})`$ using about $`100`$ different input states $`\stackrel{}{\xi }`$. The final step is to calculate the average over an ensemble of $`𝒥`$, to get $`i=\frac{1}{S}_{s=1}^Si_s`$.

## 5 Results from simulations

Here we present studies of three relevant regions in parameter space $`(\alpha ,X)`$: (1) $`\alpha `$ small ; (2) intermediate values of $`\alpha `$ and strong correlations, $`X1`$; and (3) $`\alpha `$ large. The results are compared with the theoretical predictions, eqs.(5)-(9). The self-averaging property of the MI is also analysed.

In the small $`\alpha `$ case, we calculated $`i(\alpha )`$ for values of $`\alpha `$ ranging from $`0.05`$ to $`0.2`$: we fixed $`P=10`$ and took a variable number of inputs in the interval $`N=50,\mathrm{},200`$. We repeated this process for several values of the correlation parameter $`X=0.1,0.2,\mathrm{},0.9`$. Using that $`i(\alpha =0)=0`$ when $`N\mathrm{}`$, we performed a linear regression on $`i`$, $`i(\alpha ,X)/\alpha =a_1(X)+a_2(X)\alpha `$, and obtained the coefficients $`a_1,a_2`$ as a function of $`X`$. For all values of $`X`$, $`a_11`$, which is in agreement with the theory. The function $`a_2(X)`$ is plotted in fig.1, in comparison with the theoretical prediction (eq.(5)). The good agreement indicates that the thermodynamic limit solution is a good estimation even for $`N`$ not very large, as long as one is in the low-loading expansion.

In order to evaluate $`a_3`$, we set $`a_11`$ (to avoid increasing the error through a larger number of parameters to be fitted). We analysed the MI only for the value $`X=0.5`$. From eq. 5 the theoretical value $`a_2=0.244`$ can be obtained, and from eqs. 7 and 6 we have $`a_3^{An}=0.227`$ and $`a_3^{RS}=0.063`$ respectively. The simulation itself was done using several values of $`P`$, and for each of them a linear regression was performed using $`[i(\alpha )/\alpha 1]/\alpha =a_2+a_3\alpha `$.

In this case the corrections due to finite size effects are noticeable. To observe convergence to the asymptotic regime we considered several values of $`P`$, namely $`P=3,5,7,10,12`$ and $`15`$. We averaged for each over, respectively, $`S=2500,1000,200,100,50`$ and $`20`$ samples of $`𝒥`$. The number of samples $`S`$ was taken larger as the number of outputs $`P`$ became smaller, to compensate for the lack of statistics. Numerical evaluation of $`a_3`$ is extremely costly in computational time; this is the reason why we restricted ourselves to a single value of $`X`$ and did not consider values of $`P`$ larger than $`15`$.

The result is that our simulation is in agreement with the analytical prediction. As we see in fig.2b, $`a_3^P`$ converges to $`a_30.28`$ as $`P`$ increases. Given the error bars, this result is compatible with $`a_3^{An}`$ but excludes $`a_3^{RS}`$. For the sake of comparison we have included fig. 2a where the same numerical analysis is done for $`a_2`$.

The results for strongly correlated inputs are shown in fig. 3a. There is good agreement between the simulation and the exponential behaviour of eq.(8): $`\nu =0.6622/3`$.

The results for the large $`\alpha `$ limit are presented in fig.3b. The logarithmic behaviour predicted in eq.(9) is observed.

To verify that the MI is self-averaging we calculated its mean-square deviation $`\mathrm{\Delta }(i)`$ over the samples and made a fit to the form $`\frac{1}{\sqrt{N}}`$. The good agreement with this expression can be seen in fig. 3c. This shows that the methods of statistical mechanics are appropriate to studies of the information in binary channels.

## 6 Conclusions

Our main result refers to the comparison between the $`RS`$  and analytical  solutions. The difference between them can be seen in a small-$`\alpha `$ expansion. Numerical simulation confirms that the two solutions coincide up to second order. At the next order the solutions are different (see eqs. 6-7), and the simulation excludes the replica calculation while it is in agreement with the analytical one (see fig. 2).

We have also verified the conjecture that the $`RS`$ solution is a good interpolation between the small and the large $`\alpha `$ behaviors . In particular, the simulation shows that for intermediate values of $`\alpha `$ and strongly correlated inputs the MI behaves as $`i\alpha ^{2/3}`$ (see fig.3a). Moreover, in fig.3b we see that the expected logarithmic behaviour for the large $`\alpha `$ case fits the MI very well.

## 7 Acknowledgments

This work was supported by a Spanish grant PB 96-47. Antonio Turiel is financially supported by an FPI grant from the Comunidad Autónoma de Madrid, Spain. M. Maravall was supported by a Beca de Colaboración from the Spanish Ministry of Education. D.R.C. Dominguez thanks the K.U. Leuven for a research fund (grant OT/94/9).


text:
# The distance to the Vela pulsar gauged with HST parallax observations.1footnote 11footnote 1Based on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555.

## 1 Introduction

Assessing distances to Isolated Neutron Stars (INSs) is a very challenging task which has been pursued using different techniques in different regions of the electromagnetic spectrum. Distances to pulsars with no glitching activity can be obtained through radio timing techniques (Bell, 1998). However, only millisecond radio pulsars allow for positional accuracies high enough to be used for parallax measurements. Toscano et al. (1999) count 6 such cases in their list of 12 pulsars, the distances to which have been determined through parallax. For the remaining 6 objects, which are all classical pulsars with periods of few hundreds msec, a VLBI approach was used, requiring a suitable nearby calibrator. The difficulties of the VLBI technique, including the need to account for changing ionospheric conditions, are apparent from the significant revisions already published for two of the six parallax values. The parallax of PSR B0919+06 went from 0.31 $`\pm `$ 0.14 mas (Fomalont et al., 1999) to 0.83 $`\pm `$ 0.14 mas (Chatterjee et al., 2001), while for PSR B0950+08 the parallax went from 7.9 $`\pm `$ 0.8 mas (Gwinn et al., 1986) to 3.6 $`\pm `$ 0.3 mas (Brisken et al., 2000).
Even if limited to twelve objects, e.g. $`<`$1% of the pulsar family (Camilo et al., 2000), determining model-independent distances of nearby pulsars is a rewarding exercise. As summarized by Campbell et al. (1996) and Toscano et al. (1999), this allows to trace the local interstellar medium. A distance value, coupled with the pulsar dispersion measure, yields the electron density along the line of sight, to be compared with the model of the galactic $`N_e`$ distribution (Taylor & Cordes, 1993). Such a model is used to derive the distances to all of the remaining pulsars ($`>`$99% of the population). Moreover, a distance value transforms the pulsar proper motion into a firm transverse velocity, to be compared with the average pulsar 3-D velocities obtained by Lyne & Lorimer (1994) and Cordes & Chernoff (1998) on larger samples.
X-ray astronomy provides hints to the distance of the score of pulsars detected so far (Becker & Trümper, 1997) through the measurement of the absorption of their soft X-ray emission. Unfortunately, distances derived from X-ray absorption are as uncertain as those derived by the dispersion measure. In general, the radio distances are greater than the X-ray ones.
With the detection of pulsars in the optical (see e.g. Caraveo 2000), distance measurements became possible using classic optical astrometry techniques, based on parallax measurements. Measuring tiny parallactic displacements is never easy, and the task can become really challenging when the targets are intrinsically faint, like isolated neutron stars. It requires high angular resolution coupled with high sensitivity, rendering HST the instrument of choice to measure proper motions and parallaxes of faint objects. The astrometric capabilities of HST were used to obtain model free measurements of the distance to Geminga (Caraveo et al., 1996) and RXJ 1856$``$3754 (Walter, 2001), two radio quiet INSs for which radio astronomy could not provide any input.
In this paper we address the Vela pulsar (PSR B0833$``$45) which provides a compelling case of nearby INS with a highly disputed value of the distance and a relatively bright optical counterpart ($`V23.6`$). The value of 500 pc, tentatively obtained after the original pulsar discovery (e.g. Milne 1968; Prentice & ter Haar 1969), has been assumed as a reference for both the pulsar and its surrounding SNR. As such, it is still quoted in radio pulsar catalogues 222see e.g. http://www.atnf.csiro.au/Research/pulsar/psr/archive/
or http://pulsar.ucolick.org/cog/pulsars/catalog/, in spite of the doubts raised by several independent investigations carried out at different wavelengths, targeting both the pulsar and the SNR. On the basis of their analysis of the ROSAT data, Page et al. (1996) placed the Vela pulsar at 285 $`\pm `$ 120 pc, while Pavlov et al. (2001a), using Chandra data fitted by a two component model, obtain a distance of $`210\pm 20`$ pc. On the other hand, studying the Vela remnant Cha et al. (1999) and Bocchino et al. (1999) find a distance of $`250\pm 30`$ pc and $``$ 280 pc, respectively. A different view is proposed by Gvaramadze (2001) who, on the basis of the pulsar scintillation velocity and on its very uncertain interpretation, sees no reasons to revise the canonical 500 pc distance.
We started our observations aimed at the measurement of the Vela pulsar parallax in 1997, during HST observing Cycle 6, but we had to wait until Cycle 8 to complete our program. Meanwhile, our data have also been used to reassess the proper motion of the Vela pulsar (De Luca et al., 2000a) and to quantify the overall reliability of our astrometric approach (De Luca et al., 2000b). In the following we shall report the analysis of our HST data leading to the measurement of the Vela pulsar annual parallax and proper motion. The impact of our result on the current understanding of the Vela pulsar is also discussed. While HST was collecting the appropriate set of images, VLBI in the southern hemisphere came into existence and Vela was one of the targets. The preliminary radio results (Legge, 2000) will be compared to our optical ones.

## 2 The observations

The measurement of the annual parallax of a star through optical astrometry techniques requires a set of at least three observations of the field, preferably taken at the epochs of the maximum parallactic elongation. For the Vela pulsar ($`\alpha _{J2000}=08^h35^m20^{\prime \prime }.6`$, $`\delta _{J2000}=45^{}10^{}35^{\prime \prime }.1`$), the epochs of maximum parallactic elongation coincide with days 118 and 303 of the year for right ascension, and days 20 and 204 for declination, with relative parallax factors (see e.g. Murray, 1983) at the maximum elongation of $`0.97`$ and $`0.90`$, respectively (see Figure 1). Although the relative parallactic factor at maximum elongation is somewhat larger in right ascension, observations at one of the corresponding epochs turned out to be difficult to schedule, due to their very tiny visibility window. For this reason, the observations of Vela were scheduled close to the epochs of the maximum parallactic displacements in declination.
Our program was originally approved for HST Cycle 6 but, unfortunately, only two observations of the planned triplet were executed. The whole program had to go through a new approval cycle and was rescheduled and successfully completed in Cycle 8. We thus obtained a total of five observations of the field with the WFPC2 between June 1997 and July 2000. The complete journal of the observations is summarized in Table 1. At each epoch, two exposures of the field were acquired with the WFPC2 “V” filter $`F555W`$ ($`\lambda =5252`$Å,$`\mathrm{\Delta }\lambda =1225`$Å) and with similar integration times. In order to maximize the angular resolution, in all cases the pulsar optical counterpart was centered on the Planetary Camera chip (PC) of the WFPC2 (pixel size of 0.045 arcsec).

## 3 Data reduction and analysis

The data reduction was been performed using the IRAF/STSDAS package. After the standard pipeline processing of the frames (debiasing, dark subtraction, flatfielding), which was performed using the most recent reference files and tables, each couple of coaligned images was co-added for a first filtering from cosmic ray hits. Residual hits were later rejected using specific cosmic ray subtraction algorithms in IRAF. Figure 2 shows the resulting image for July 2000.
The cleaned images have then been used for the definition of a relative reference frame to be used as a starting point for our astrometric procedure. Since all the observations have been taken with different telescope roll angles and with small relative offsets, the definition of a relative reference frame must rely on a very accurate image superposition. Following the approach applied successfully in previous astrometric works (e.g. Caraveo et al. 1996; Caraveo and Mignani 1999; De Luca et al. 2000a; Mignani et al. 2000a, Mignani et al. 2000b), this is done by computing a linear coordinate transformation (i.e. accounting for 2 independent translation factors, 2 scale factors and a rotation angle) between a set of reference objects.
The selection of the reference grid objects is critical. They must be point-like, present in all observations (but not too close to the field edges), bright enough to allow for an accurate positioning (but not saturated). A set of 26 such objects (labeled in Figure 2) was selected in the reference frame of the PC.
Since the shape of the PC Point Spread Function (PSF) is known to be position-dependent (Krist, 1995), we did not use simulated PSFs for source fitting. On the other hand, the number of good reference stars was not sufficient to compute a template PSF directly from the images. Thus, the reference object coordinates were computed by fitting a 2-D gaussian function to their intensity profiles, using optimized centering areas. This yielded positional uncertainties of the order of $`0.01÷0.05`$ pixel, depending on the objects’ brightness and position on the chip. Special care was devoted to characterizing the errors involved in the centroid determination. Following De Luca et al. (2000b), we addressed both statistical errors (i.e. due to each object’s S/N) and systematic ones (due, e.g., to the telescope jitter, defocussing of the camera, charge transfer in the CCD, background fluctuations, etc.). Moreover, we checked that our results were not biased by the algorithm used for the fitting. The coordinates were corrected for the “34th row error” (Anderson & King, 1999) and for the significant, well-known, instrument geometrical distortion using the most recent mapping of the PC field of view (Casertano & Wiggs, 2000). The centroids of the Vela pulsar optical counterpart were obtained in the same way, yielding errors ranging between 0.02 and 0.04 pixel.
Having secured a reference grid, we registered all the frames on the June 1999 one, taken as a reference and previously aligned along right ascension and declination according to the telescope roll angle. The rms of the residuals on the reference object coordinates were $`<0.05`$ pixel in right ascension and $`<0.04`$ pixel in declination. The overall accuracy of the frame registration, accounting for errors in the centroid determination as well as in the geometric distortions mapping, and the accuracy of the fit, has been discussed in detail by De Luca et al. (2000b).
To ensure that our procedure was not affected by any displacement of our 26 reference objects (due either to proper motions or parallaxes), we have repeated the frame registration 26 times, i.e. each time excluding one of the objects from the fit. In addition, to exclude any possible bias due to the arbitrary choice of the reference frame, we have repeated the whole procedure cycling it over the five epochs. In all cases, we obtained statistically undistinguishable results. We are thus confident that our procedure is correct and free of systematics. Last, we have applied the coordinate transformations to the positions of the Vela pulsar; the resulting relative positions are shown in Figure 3.
If fitted with a simple proper motion, all the points in Figure 3 are seen to deviate from the straight line. Their residuals wrt the proper motion fit, however, are not randomly distributed. Rather, they follow the trend expected for an object affected also by parallactic displacement.

## 4 Analysis of the pulsar displacements

The geocentric right ascension and declination $`(\alpha _G(t),\delta _G(t))`$ of the pulsar at a given epoch $`t`$ can be expressed in the form

$$\{\begin{array}{cc}\alpha _G(t)=\alpha _B(t_0)+\mu _\alpha cos(\delta )(tt_0)+\pi P_\alpha (t)\hfill & \\ \delta _G(t)=\delta _B(t_0)+\mu _\delta (tt_0)+\pi P_\delta (t)\hfill & \end{array}$$
(1)
where $`(\alpha _B(t_0),\delta _B(t_0))`$ are the barycentric coordinates at a reference time $`t_0`$, $`(\mu _\alpha ,\mu _\delta )`$ are the right ascension and declination components of the proper motion, $`\pi `$ is the annual parallax and $`(P_\alpha (t),P_\delta (t))`$ are the parallactic factors shown in Figure 1. A system of equations like ( 1) can be written for each of the five epochs corresponding to our observations.
Setting $`t_0`$=1999, June 30 as the reference epoch, one can obtain five pairs of equations relating the observed coordinates of the pulsar to the 5 unknowns $`\pi `$, $`\mu _\alpha `$, $`\mu _\delta `$, $`\alpha _B(t_0)`$ and $`\delta _B(t_0)`$. A least squares fit yields for the parallax and the proper motion the following values: $`\pi =3.4`$ mas, $`\mu _\alpha cos(\delta )=45.0`$ mas yr-1, $`\mu _\delta =25.8`$ mas yr-1, corresponding to a reduced $`\chi ^2`$ of 0.11 (5 degrees of freedom).
To evaluate the uncertainties and the confidence levels on the fitted parameters, we ran a MonteCarlo simulation for a theoretical source featuring proper motion and parallax values equal to our best fit ones. Each synthetic data set was obtained by perturbating a coordinate set representing the expected source geocentric positions at the epochs of our observations. The experimental error (per coordinate) for each data point in each simulated data set was computed from a normal distribution with a standard deviation equal to the overall uncertainty (per coordinate) affecting the pulsar positioning estimated in sect. 3 (i.e. $``$2 mas per coordinate in 1999, June 30; $`2.5÷2.9`$ mas per coordinate in the remaining epochs).
After $`10^5`$ simulations, we estimate the 1 $`\sigma `$ error bar, to be attached to the parallax value, at 0.7 mas, while the uncertainties on the proper motion are 1.1 mas in $`\alpha `$ and 1 mas in $`\delta `$. Thus, our best evaluation of the Vela pulsar displacements is as follows:

$`\pi =3.4\pm 0.7`$ mas

$`\mu _\alpha cos(\delta )=45.0\pm 1.1`$ mas yr-1

$`\mu _\delta =25.8\pm 1.0`$ mas yr-1

corresponding to a position angle of $`299.8^{}\pm 1.2^{}`$.
The Vela path in the sky, as predicted from the best fitting proper motion and parallax, is plotted in Figure 4, together with the measured pulsar positions. The agreement between the expected positions and the measured ones is remarkably good.

### 4.1 The proper motion

Our best fit to the pulsar proper motion improves the result of De Luca et al. (2000b) by confirming the overall value but by further reducing its associated errors. Since we measure the Vela proper motion with respect to a set of reference stars in the field, our value is sensitive to the peculiar motion of the Sun in the Local Standard of Rest. This induces an apparent annual displacement of the pulsar in the anti-Apex direction, which should be corrected for in order to evaluate the pulsar motion with respect to the structures located in its immediate surroundings (such as the synchrotron nebula detected in X-rays). A word of caution is required, since the correction depends critically on a reliable definition of the LSR. Following the analysis of Dehnen & Binney (1998), based on Hypparcos data, we assumed a peculiar motion of the Sun of $`13.4\pm 0.8`$ km/s in the direction $`l=28\pm 3^{}`$, $`b=32\pm 2^{}`$. The components of the proper motion of the Vela pulsar with respect to the LSR become $`\mu _\alpha cos(\delta )=37.2\pm 1.2`$ mas yr-1, $`\mu _\delta =28.2\pm 1.3`$ mas yr-1, for a resulting position angle of $`307.2\pm 1.6^{}`$. We note that a very similar result ($`\mu _\alpha cos(\delta )=38.3\pm 1.2`$ mas yr-1, $`\mu _\delta =27.4\pm 1.3`$ mas yr-1, P.A.=$`305.6\pm 1.6^{}`$) is obtained using the correction for the Sun peculiar motion adopted by Bailes et al.(1989).

### 4.2 The annual parallax

The value of the pulsar annual parallax which best fits our data is $`3.4\pm 0.7`$ mas. The probability of obtaining such a value in absence of any parallactic displacement is of the order $`10^5`$. From our parallax value we can derive the distance to the Vela pulsar, which turns out to be:

$$294_{+76}^{50}pc$$

We can now compare our results on the Vela proper motion and parallax with the preliminary radio ones obtained by Legge (2000), e.g., $`\mu _\alpha cos(\delta )=49.8\pm 0.2`$ mas yr-1, $`\mu _\delta =30.5\pm 0.1`$ mas yr-1 and $`\pi =3.13\pm 0.33`$ mas. While his parallax value is similar to ours, his proper motion values are not compatible with those obtained using the HST both for the study of ”pure”, parallax-free, proper motion (De Luca et al., 2000b) and for our parallax fit. However, no details are available on the calibrator used, nor on the ionospheric correction approach, nor on the correction for galactic rotation (on top of the Sun peculiar motion), so that we cannot assess the accuracy of the radio data analysis. Since the optical approach we have used is free from all such systematics, we conclude that our distance determination based, as it is, on direct measurement, with well known error determination procedure, is to be taken as the most reliable estimate to the true distance of this celestial object.

## 5 Discussions

After Geminga (Caraveo et al., 1996) and RXJ1856-3754 (Walter, 2001), this is the third measurement of the optical parallax of an isolated neutron star, to be added to the list of a dozen radio parallaxes summarized by Toscano et al. (1999). At variance with the two previous cases, that had no firm distance estimates, the new value for the distance of the Vela pulsar is significantly smaller than the traditionally accepted one, confirming the earlier claims by Page et al. (1996), Bocchino et al. (1999), Cha et al. (1999) and Pavlov et al. (2001a).
While adding an important piece of information to the distribution of electrons in the local environment, a distance of $``$300 pc has several implications for the Vela pulsar physics, as well as for its kinematics. In what follows we shall address each of these points.

### 5.1 Local Interstellar Medium

Toscano et al. (1999) used the 12 radio pulsars with a measured parallax to map the local interstellar medium and deduce electron densities along their different lines of sight. The Vela pulsar was added to the sample on the basis of the distance of 250 pc, suggested by Cha et al. (1999) for the SNR. Coupling the Vela DM (67 cm-3) with such a distance, Toscano et al. (1999) computed an $`N_e`$ of 0.270 cm-3, to be compared with the average density of 0.02 cm-3 of the Taylor & Cordes (1993) model for the local region. Indeed, the $`N_e`$ towards Vela is the maximum in their sample and it is 5 to 10 times higher than the values found for four other pulsars in the 3rd galactic quadrant. These pulsars already show $`N_e`$ values systematically higher than a comparable number of objects in the first galactic quadrant. Our distance value, while slightly lowering the $`N_e`$ value to 0.23 cm-3, confirms the extremely high electron density towards Vela, possibly pointing towards the existence of ionised clouds in the Vela region direction, at the edge of the so called $`\beta `$ CMa tunnel. This problem has been recently addressed by Cha et al. (2000) who, doing the ”astronephography” of the region (3-D mapping of interstellar clouds), found three distinct absorption systems with different velocities, consistent with an enhanced density of ionized material towards PSR B0833-45.

### 5.2 Multiwavelength emission

For twenty years, up to 1993, the Crab and Vela pulsars were the only neutron stars detected throughout the entire electromagnetic spectrum, from radio to optical to high-energy $`\gamma `$-rays. Remarkably similar in high energy $`\gamma `$-rays, the behavior of the two pulsars is very different at all other wavelengths (see e.g. Thompson. 2001). Furthermore, our current understanding of the pulsar emission mechanisms at wavelengths other than radio has been biased by the phenomenology of Crab and Vela, since theories have been shaped to account for the multiwavelength behavior of these two objects. The significant downsizing of the Vela pulsar luminosity (by a factor of $`3\pm 1`$) presented here, implies some revision on the current view of pulsars multiwavelength behavior. However, lowering the distance to the Vela pulsar does not have the same impact in the various spectral domains.
Vela has a very special place in the high energy $`\gamma `$-ray sky where it is by far the brightest source, outshining the Crab by a factor of $``$4. Its radiation is totally pulsed (Kambach et al., 1994) and its spectral shape varies throughout the light curve. Of course, to convert the measured flux into a luminosity one needs to know the pulsar distance, together with its beaming factor. While the $`\gamma `$-ray community never disputed the traditional 500 pc value for the Vela distance, the beaming solid angle remains unknown, ranging from the size of a neutron star polar cap to $`4\pi `$ for isotropic emission (see Thompson et al. 1999 and references therein). For simplicity, a 1 sr conical beam is generally assumed, yielding a beaming factor of $`1/4\pi `$. With such a beaming value, the high energy luminosity of the Vela pulsar is now $`7\times 10^{33}`$ erg sec-1, comparable to that of the much older PSR B1055$``$52 and 10 times lower that that of PSR B1706$``$44, a pulsar remarkably similar to Vela in its $`P`$, $`\dot{P}`$ and overall energetics. Moreover, the efficiency with which Vela converts its rotational energy loss into $`\gamma `$-rays becomes 0.001, similar to that of the much younger Crab and 20 times smaller than that of PSR B1706$``$44. Not surprisingly, the Vela pulsar is now well below any best fitting line correlating the $`\gamma `$-ray luminosity to pulsar parameters such as the number of accelerated particles (Thompson et al., 1999) or the value of the open field line voltage (Thompson, 2001). This will require a critical re-examination of the position of Vela in the $`\gamma `$-ray emitting pulsar family. Vela may become an underluminous pulsar, an apparent paradox considering its brightness in the $`\gamma `$-ray sky.
In the X-ray domain, the situation is different. The X-ray emission phenomenology already prompted a number of authors (Page et al. 1996; Pavlov et al. 2001a; Helfand et al. 2001) to opt for a distance smaller than the “canonical” one. Indeed, for a distance of 500 pc, the soft X-ray flux, most probably originated by the hot surface of the neutron star, pointed towards an emitting area ($`R^{\mathrm{}}=3÷4`$ km) far too small to be compatible with the whole neutron star surface (Ögelman et al., 1993). On the other hand, a very small emitting hot spot is obviously incompatible with the shallow pulsation seen at these energies. Reducing the pulsar distance to half its “canonical” value, Pavlov et al. (2001a) find that for a pure black-body emission, a radius of the emitting area of just $`1.9÷2.4`$ km could explain the flux observed by Chandra. Introducing a modified black body with a pure H atmosphere eases the emitting area problem, allowing for the whole surface of a 13 km radius neutron star at 210 pc to emit at a temperate significantly smaller than that obtained in the pure black body case. Our new independent distance determination now freezes one of the parameters of the model, allowing for a better determination of both the neutron star radius and its temperature.
The parallax distance fixes the overall Vela X-ray luminosity (0.2-8 keV) to be about 5.5 1032 erg/sec, divided between a thermal component, dominating at low energies, and a non thermal one, emerging only at higher energies. If compared to the available rotational energy loss, such an X-ray emission accounts for 8 10-5 of the pulsar reservoir, significantly less than the average ratio of $``$ 10-3 found for the majority of the young X-ray emitting pulsars (Becker & Trümper, 1997). Moreover, such an already low conversion efficiency shrinks to 10-6 if one considers only the non-thermal component, confirming the severe underluminosity of PSR 0833-45 in the X-ray domain. The interpretation of such an underluminosity is, however, complicated by the composite nature of the Vela emission. It is certainly dominated by thermal emission for E $``$ 1.8 keV, with a weaker magnetospheric component emerging only at E $``$ 1.8 keV Pavlov et al. (2001a). Because of such a composite nature, the Vela pulsar X-ray emission cannot be readily compared to other purely magnetospheric cases (e.g the Crab and PSR B0540$``$69).
The X-ray data can also be used to further map the local ISM. The comparison between distance and Hydrogen column density, $`N_H`$, responsible of the X-ray absorption, provides an assessment of the reliability of this parameter as a distance indicator for nearby sources.
In the optical domain, the reduction in the pulsar luminosity has important implications. Here the emission is certainly magnetospheric - as witnessed by its double-peaked light curve (Wallace et al., 1977; Manchester et al., 1978; Gouiffes, 1998) and by its flat $`3600÷8000`$ Å spectrum (Nasuti et al., 1997; Mignani & Caraveo, 2001). However, the revised value of the B-to-U optical luminosity $`L_{BU}`$ $``$ 5.5 1028 erg sec-1 is now quite low with respect to the prediction of the classical Pacini & Salvati (1983) model, based on synchrotron emission at the light cylinder. Such an underluminosity could neither be ascribed to the spectral shape of the optical emission, which is remarkably similar to the Crab one (Mignani and Caraveo, 2001), nor to the shape of the light curve (see e.g. Gouiffes, 1998 ), which covers a broader phase interval than the Crab one, resulting, if anything, in an increase of the optical output.
The behaviour of the Vela pulsar in the optical domain should be considered in the light of its “transition” position between the group of the young energetic pulsars, characterized by pure magnetospheric emission and the middle-aged ones characterized by composite spectra and lower emission efficiency (see e.g. Caraveo, 1998; Mignani and Caraveo, 2001).
In the X-ray domain Vela behaves as a middle-aged object (see e.g. Pavlov et al, 2001a) while in the optical and gamma-ray energy ranges its phenomenology is reminiscent of that of the younger pulsars but with a low emission efficiency. The multiwavelength luminosities of the Vela pulsar, stemming from the newly determined distance, will help to discriminate the emission mechanisms at work in the different energy domains and to assess their evolution throughout the life of the pulsars.

### 5.3 Proper motion

The parallax fit yields also, as a by-product, the best proper motion value obtained so far at optical wavelengths, corresponding to a transverse velocity of $``$ 65 km sec-1. The space direction of the Vela proper motion has been correlated with the axis of symmetry of the Chandra X-ray structure by a number of authors (Mignani et al. 2000b, Helfand et al. 2001 and Pavlov et al. 2001b). Since the X-ray nebula axis of symmetry should trace the pulsar rotational axis, an alignment between such axis of symmetry and the pulsar proper motion would have important bearings of the kick mechanism responsible for the neutron star ejection. Lai et al. (2001), noting that proper motion-axis of symmetry alignment seems to be present also for the Crab pulsar (Caraveo & Mignani, 1999), considered different spin-kick alignment mechanisms. Also in view of the pulsars’ low velocities, they concluded that both natal kicks, be they hydrodynamic or neutrino-magnetic, as well as postnatal ones, such as the Harrison and Tademaru (1975) electromagnetic kick, could account for the alignments.
Here, we reassess the case for the Vela proper motion X-ray jet alignment by comparing our ”corrected” proper motion position angle (p.a.) of $`307.2^{}\pm 1.6^{}`$ with the value of $`307^{}\pm 2^{}`$ computed by Pavlov et al. (2001b) for the axis of symmetry of the structure seen by Chandra. Assuming an uniform probability distribution for the two vectors in space, we estimate the chance occurrence probability to be 1.2 $`10^3`$, thus strengthening the case for an alignment between the Vela rotation axis and its proper motion. Figure 5 shows an update of Figure 2 of Mignani et al. (2000b), with our corrected proper motion vector superimposed to the high resolution Chandra image of the Vela nebula. Taking advantage of such an alignment we can use the angle of 53.2, computed by Helfand et al (2001) between the axis of the X-ray torus and the line of sight, to convert our measured transverse velocity into the pulsar 3-D velocity. Thus, the space velocity of the Vela pulsar turns out to be 81 km/sec. While hardly compatible with the Lyne & Lorimer (1994) mean value of 400 km sec-1, our value is also in the low side of the low velocity component of the bimodal fit by Cordes & Chernoff (1998), classifying Vela as a rather slow moving pulsar.

## 6 Conclusions

We have presented the first direct measurement of the distance to the Vela pulsar at optical wavelengths. Our value is significantly smaller than the canonical one, implying a downgrading of the luminosity of the Vela pulsar, positioning Vela among the underluminous pulsars. High energy $`\gamma `$-ray astronomy is the most affected by our result since it pertains to the brightest, hence most conspicuous, source in the $`\gamma `$-ray sky. The optical domain is also seriously affected by our overall luminosity reduction, while the X-ray side is only marginally altered. The distance determination also confirms the anomalously high density of electrons towards the Vela pulsar, pointing to the presence of ionized clouds in the local interstellar medium.
As a byproduct of the parallax fit, we have obtained a very accurate proper motion value, which appears nicely aligned along the axis of symmetry of the X-ray nebula, hence the pulsar rotation axis. Moreover, our proper motion, together with our distance determination, fixes the transverse velocity of Vela to 65 km sec-1. Using the inclination of the X-ray torus, as measured by Helfand et al (2001), we can convert our transverse velocity into the 3-D one, yielding a value of 81 km sec-1. For a pulsar, this is certainly a low velocity.

## 7 Acknowledgements

This work was supported by the Italian Space Agency (ASI).
We thank the referee, David Chernoff, for his suggestions.


text:
# Acknowledgements

## Acknowledgements

It is a pleasure to thank the organisers of the Second Winter Workshop on Mathematical Physics, University of Queensland, July 2001, for providing a stimulating and multidisciplinary environment. JDB acknowledges the support of the University of Tasmania Institutional Research Scheme, grant number RDOJ0011872. PDJ thanks David Penny and Michael Hendy for hospitality and discussions on a visit at the Centre for Biomolecular Sciences, Massey University. Finally we thank Bob Delbourgo and Vladimir Rittenberg for fruitful comments.


text:
# The DLAs contribution to the metal content of the high z Universe

## 1. Introduction

Much effort has been devoted to the study of chemical abundances in the DLAs in the last few years (Lu et al. 1996, Centurión et al 2000, Prochaska $`\&`$ Wolfe 1999, Ellison et al. 2001, Molaro et al. 2001). Most of these works is new observational data with subsequent analysis, which is based in comparison to known patterns of abundance and abundance ratios, mainly the ones observed in our Galaxy. These studies have already revealed some important characteristics of the DLAs which gives us some clues about their nature: their metallicity, given generally by \[Zn/H\], lies in the range -2 to 0 (Vladilo 1998, Lauroesch et al. 1996) suggesting different epochs for DLAs formation or different star formation histories; there is no sign of evolution in the metallicity of DLAs from N(HI)-weighted Fe abundance (Prochaska $`\&`$ Wolfe 1999, Prochaska $`\&`$ Wolfe 2000), altought the metallicity content does decrease with redshift (Savaglio 2000); the \[Zn/Fe\] ratio is mostly above solar indicating the presence of dust depletion in these systems (Pettini et al. 1997); the \[$`\alpha `$/Fe\] pattern observed in DLAs roughly resembles that of metal-poor stars of our Galaxy (Pettini et al. 2000, Ellison et al. 2001), even though, there are some systems with values similar to those predictions by chemical evolution models for dwarf galaxies (Centurion et al. 2000, Molaro et al. 2001).

Among these works, however, just in a few of them chemical evolution models are used in the comparison to the observational data (Prantzos $`\&`$ Bossier 2000, Matteucci et al. 1997). The predictions of the chemical evolution models are compared to absolute abundances (Prantzos $`\&`$ Boissier 2000) and abundance ratios (Matteucci et al. 1997). The abundance ratios are better suited for chemical evolution studies because, unlike absolute values, they do not depend heavily on model assumptions, but mainly on the stellar nucleosynthesis and on the adopted IMF. The \[$`\alpha `$/Fe\] and \[N/$`\alpha `$\] ratios can be used to settle constraints for star formation timescales due to the different timescales for the formation of these elements. While the $`\alpha `$ elements are produced mainly in type II Supernovae (SNe II) in short timescales, the iron peak ones and N are produced in SNe Ia and in intermediate massive stars (IMS) in a larger timescale. There is a controversy about the observed pattern of the \[$`\alpha `$/Fe\] ratio in DLAs: while some authors detected patterns similar to that of the metal-poor stars of our Galaxy (Pettini et al. 2000, Ellison et al. 2001), others found values that resemble the chemical evolution of dwarf galaxies (Centurion et al. 2000, Molaro et al. 2001). This discrepancy may be associated with dust depletion.

We use, in this work, one-zone chemical evolution models to represent the DLAs, which predictions are compared to \[$`\alpha `$/Fe\] and \[N/$`\alpha `$\] ratios. The abundance of Zn, Fe and Si are corrected for dust depletion and robust statistical methods are used to establish the trends of the data. In this way, we try to make a complete study of the DLAs chemical evolution embracing the most important aspects that might affect these kind of analysis.

## 2. The Chemical Evolution Models

Due to uncertainties about the nature of DLAs, two types of one-zone model are used in a first study about chemical enrichment in these systems. A model with infall of pristine gas is used to represent the disk systems and another one with galactic winds for dwarf galaxies.

### 2.1. The disk model

The disk model belongs to the class of the models with one infall episode of pristine gas. We used the basic equations for chemical evolution, with a star formation rate (SFR) similar to the one used by Matteucci $`\&`$ François (1989), and a Salpeter IMF. The normalisation $`\stackrel{~}{\nu }`$ of the star formation law is made with the solar position at 8 kpc and a galaxy age of 13 Gyr. The stellar yields for intermediate mass stars (IMS) come from Renzini $`\&`$ Voli, 1981 (hereafter RV81), models with $`\alpha _c`$ = 0 and 1.5. The yields from van den Hoek $`\&`$ Groenewegen (1997) for the IMS were also tested but the models with these yields produce less nitrogen than both $`\alpha _C`$ = 0 and 1.5 models of RV81 and consequently are not applied in this DLAs study. The infall rate decreases with time and the infall timescale and superficial mass density profile vary with the galactic radius.

We run several models with different values of r: 2, 4, 8, 14, 18 kpc.

### 2.2. The dwarf galaxy model

The dwarf galaxy model belongs to the class of models with galactic winds. The majority of the models used here are similar to the classic ones, in which a wind is established in a time $`t_w`$ when the thermal energy of the gas exceeds the potential energy.

A galaxy with 1 Gyr is represented by the models here applied, but there is also a model which represents a galaxy with $`10^7`$ yr. The galaxy is inside a dark halo with a mass three times greater than the initial galaxy’s baryonic mass. The star formation is continuous with a SFR proportional to the present mass of gas and characterised by the specific star formation rate $`\nu `$ (the inverse of the star formation timescale), which varies in each model from 0.1 to 19 $`Gyr^1`$. It is adopted a Salpeter IMF between 0.1 and 100 M.

Besides the classic models with wind, an early wind model (Friaça 2000) is applied. Friaça (2000) proposed differential winds, due to SNeII, in dwarf galaxies in order to explain the \[$`\alpha `$ /Fe\] below solar in the intercluster medium.

## 3. The Data

The observational data about the DLAs were collected from the literature and include chemical abundance (represented by \[X/H\]= $`log(X/H)_{obs}`$ \- $`log(X/H)_{}`$) of a variety of elements ranging from N to Zn, observed in the redshift interval $`0.6<z<4.4`$ (for more details and references see Lanfranchi $`\&`$ Friaça 2001, submited).

## 4. Analysis of the data

The observational data are analysed with robust statistical methods and compared to predictions of the chemical evolution models. In the statistical analysis we followed the approach described by Chiappini et al. (1999). We obtain, with the summary statistics of Cleveland $`\&`$ Kleiner (1975), three lines that summarise the trend of the data: MM (midmean), LSMM (lower semi-midmean) and USMM (upper semi-midmean). The final plots represent locally weighted regression (lowess) smoothes of the summary lines.

The trend lines and the points distribution are then compared to the predictions of the models described above focusing the \[$`\alpha `$/Fe\] and \[N/$`\alpha `$\] ratios. Due to the \[Zn/Fe\] generally above solar in our sample, what may be an indication of the presence of dust in the DLAs, we corrected the abundances of Zn, Cr, Fe and Si by dust depletion. The correction is based in the depletion pattern observed in our Galaxy together with the observed \[Zn/Fe\] in each system analysed (see Lanfranchi $`\&`$ Friaça 2001). The \[Si/Fe\] corrected for dust is then used as a tracer of the \[$`\alpha `$/Fe\] ratio in the comparison of the data to the predictions of the models, while \[N/S\] is used to represent the \[N/$`\alpha `$\] ratio. In this way, we obtain dust free (or almost free) tools to be compared with the models.

## 5. Application of the one-zone models to DLAs

### 5.1. Disk galaxy model

The disk models predictions for \[$`\alpha `$/Fe\] and \[N/$`\alpha `$\] ratio, using the RV81 yields with $`\alpha _C`$ = 1.5, for different radius, r = 2, 4, 8, 14, 18 kpc compared to the observed values in DLAs are shown on the figure 1.

The observed \[N/S\] ratios with absolute values are centralised mainly between the lines for r = 4 and 14 kpc, with a large number of systems around the r = 8 kpc line, as indicated by the statistical trend lines. The lines trend are also similar to the models predicted behaviour: the ratio increases with \[$`\alpha `$/H\] rapidly in low metallicities and then becomes flatter.

The comparison between the observed \[N/$`\alpha `$\] in DLAs and the disk model predictions indicates that the observed values of this ratio represent a long timescale for the star formation in these systems for regions around the solar neighbourhood ($`4<r<14kpc`$).

There is a discrepancy, however, when the corrected \[Si/Fe\] ratio is used in the comparison to the models. The region of the models with $`r`$ 14 - 18 kpc fits reasonably well the majority of the systems. The same behaviour is found for the statistical lines. They exhibit a trend similar to the models predictions and are located in the region around $`\mathrm{\hspace{0.33em}14}18kpc`$. Consequently, the comparison indicates that the observed values correspond to external regions ($`r>14kpc`$) of disk systems, implying a long star formation timescale. A long and continuous star formation allows the iron peak group elements to be produced, giving rise to low values of \[Si/Fe\].

### 5.2. The dwarf galaxy model

On the figure 2, the comparison between the predictions of the dwarf galaxy models with different specific star formation rates and the \[$`\alpha `$/Fe\] and \[N/$`\alpha `$\] ratios are shown.

The predictions for \[N/$`\alpha `$\] are compared to the observed ratios both in DLAs and dwarf galaxies. The $`\alpha `$ elements abundance is determined only by oxygen (circles) in dwarf galaxies while in DLAs sulphur (squares) and silicon (triangules) are also used to estimate it. The statistical lines are the same as on previous figure.

The observed \[N/$`\alpha `$\] in DLAs show a great dispersion, larger than that found for Galactic stars in the same range of metallicity and for dwarf galaxies. The different values of N/$`\alpha `$ for a given $`\alpha `$/H can be explained with the difference between the formation timescale of N and $`\alpha `$ elements. This fact can be seen in models with diverse values of $`\nu `$ that go from almost quiescent star formation ($`\nu =\mathrm{\hspace{0.33em}0.1}Gyr^1`$) to an extreme rapid one ($`\nu =\mathrm{\hspace{0.33em}19}Gyr^1`$ for systems with $`10^9M_{}`$). Most dwarf galaxies seem to require a specific star formation ($`\nu `$) between 0.3 and 3 $`Gyr^1`$, while the DLAs seem to need lower $`\nu `$: $`0.1\nu 1`$ $`Gyr^1`$ for classic wind models and $`\nu `$ 1 $`Gyr^1`$ for models with early wind. The models with galactic winds reproduce all dwarf galaxy points but are unable to match some DLAs that fall totally out of the region delimited by these models.

In the case of the dwarf galaxy models too, there is a difference when the \[Si/Fe\] ratios is used in the comparison to the models. While for the \[N/S\] ratio, most of the observed values are fitted by models with $`0.3\nu 1Gyr^1`$, most systems and the statistical lines lie below the \[$`\alpha `$/Fe\] predicted by any model. Only the upper semi-midmean (USMM - dotted line) goes near the model with $`\nu `$ = 0.1 $`Gyr^1`$. This fact implies a specific star formation $`\nu 0.1Gyr^1`$, which represents a long star formation timescale.

## 6. Conclusions

The comparison between the one-zone chemical evolution models for dwarf and disk galaxies and the observational data indicates that these models are able to reproduce only partially the \[N/$`\alpha `$\] and \[$`\alpha `$/Fe\] observed in DLAs. This fact suggests that these systems may represent various types of galaxies and not a single one. The observed abundance ratios also require long star formation timescale in both models. The discrepancy between the comparisons using \[N/$`\alpha `$\] and \[$`\alpha `$/Fe\] might be related to an underestimated correction for Si.

## References

Centurión M., Bonifacio P., Molaro P. & Vladilo G. 2000, ApJ, 536, 540

Chiappini C., Matteucci F., Beers T.C. & Nomoto K. 1999, ApJ, 515, 226

Cleveland W.S. & Kleiner B. 1975, Technometrics, 17, 447

Ellison S.L., Lewis G.F., Pettini M., Sargent W.L.W., Chaffee F.H. & Irwin M.J. 1999, astro-ph 9903063

Friaça A.C.S. 2000, in Dwarf Galaxies and Cosmology, Proceedings of the XXXIII Recontres de Moriond, Y.T. TRINH, C. BALKOWSKI, CAYETTEV, J.T.T. VÂN, Paris, Editions Frontières, p.312

Lanfranchi G.A. & Friaça A.C.S. 2001, submited to MNRAS

Lauroesch J.T., Truran J.W., Welty D.E. & York D.G. 1996, PASP, 108, 641

Lu L., Sargent W.L.W., Barlow T.A., Churchill C.W. & Vogt S. 1996, ApJS, 107, 475

Matteucci F., Molaro P. & Vladilo G. 1997, AAp, 321, 45

Matteucci F., & François P., 1989, MNRAS, 239, 885

Molaro P., Bonifacio P., Centurión M., Vladilo G., D’Odorico S., Levshakov S.A. 2001, astro-ph

Pettini M., Ellison S.L., Steidel C.C., Shapley A. E. & Bowen D.V. 2000, ApJ, 532, 65

Pettini M., King D.L., Smith L.J. & Hunstead R.W., 1997, ApJ, 478, 536

Prantzos N. & Boissier S. 2000, MNRAS, 315, 82

Prochaska J.X., Naumov S.O., Carney B.W., McWilliam A. & Wolfe A.M. 2000, astro-ph 0008075

Prochaska J.X. & Wolfe A.M. 1999, ApJS, 121, 369

Prochaska J.X. & Wolfe A.M., 1997, ApJ, 487, 73

Renzini A. & Voli M. 1981, A&A, 94, 175

Savaglio S., Panagia N. & Stiavelli M. 1999, astro-ph 9912112

Van Den Hoek L.B. & Groenewegen M.A.T. 1997, AASS, 123, 305


text:
# Liquid, Glass and Crystal in Two-dimensional Hard disks

## Abstract

We study the thermodynamic and dynamic phase transitions in two-dimensional polydisperse hard disks using Monte Carlo methods. A conventional local Monte Carlo algorithm allows us to observe a dynamic liquid-glass transition at a density $`\rho ^G`$, which depends very little on the degree of polydispersity. We furthermore apply Monte Carlo methods which sample the Boltzmann equilibrium distribution at any value of the density and polydispersity, and remain ergodic even far within the glass. We find that the dynamical transition at $`\rho ^G`$ is not accompanied by a thermodynamic transition in this two-dimensional system so that the glass is thermodynamically identical to the liquid. Moreover, we scrutinize the polydispersity-driven transition from the crystal into the disordered phase (liquid or glass). Our results indicate the presence of a continuous (Kosterlitz-Thouless type) transition upon increase of the polydispersity.

For several decades now, enormous interest has been concentrated on the glass transition in classical liquids . At this transition, the viscosity of a structurally disordered liquid increases dramatically and reaches values typical of a crystalline solid. The underlying reasons for the spectacular slowdown of the system’s time evolution have been hotly debated and views very much diverge on this question . The transition is naturally related to the loss of thermal equilibrium during the quench towards low temperatures, and it is well-known that the properties of a glass strongly depend on the detailed quench procedure . The theoretical description of this phenomenon in all its complexity is a very difficult task.

Notwithstanding the intimate relationship of the glass transition with loss of thermal equilibrium (i.e. ergodicity breaking), it is nevertheless possible to describe the glass within thermal equilibrium. This means that one studies configurations of the system taken from the Boltzmann distribution (which exists in a mathematical sense at any temperature and pressure) and computes their thermodynamic and dynamic properties. For any thermodynamically stable glass, the reference to the quench history can thus be omitted, and the theoretical analysis is much simpler.

The most serious problem with the ‘equilibrium glass’ is of course that the Boltzmann distribution is experimentally inaccessible inside the glassy phase. Attempts to characterize the true equilibrium properties of a glass necessarily rely on the extrapolation of experimental data.

However, computational methods have been developed which allow to bypass this obstacle: non-local Monte Carlo (MC) algorithms can remain ergodic in the glassy phase for prominent model systems, such as hard spheres, and allow to compute their equilibrium observables directly . With this approach, we may address questions at the heart of long-standing controversies, for example about the existence of an inaccessible thermodynamic phase transition within the glass, which has been rendered responsible for the dynamic transition .

In the past, many classical statistical physics models have been considered in the study of glasses. For a long time, it was believed that the very simplest model systems (monodisperse hard spheres, hard disks, and central potential models) did also possess a glassy phase . However, more recent MC simulations have shown that monodisperse hard spheres always crystallize, and that the previously observed behavior was mainly due to finite-size effects and limited computer resources . This has also been confirmed in experiments with colloidal systems, where crystallization preempts the glass transition for monodisperse hard spheres. In order to avoid parasitic effects, some of the experiments have been performed in zero-gravity environments .

The situation changes abruptly, both in experiment  and in the simulations, in the presence of a small dispersity in size. Several investigations have shown the existence of a so-called ‘terminal polydispersity’ of a few percent (cf, e.g., ), above which crystallization can be avoided for any value of the external pressure, i.e., up to very high densities.

Polydispersity opens up opportunities to study two very interesting, yet distinct phenomena: $`(i)`$ the glass transition, which is no longer preempted by crystallization at sizeable values of the polydispersity; $`(ii)`$ the thermodynamic transition of the crystal into a disordered state (liquid or glass) upon an increase of the polydispersity.

In the present work, we study these two phenomena in two-dimensional hard disks. We investigate their dynamical and thermodynamical phase diagram as a function of density $`\rho `$ (or pressure $`P`$) and of the polydispersity for a given functional form of the size dispersion. At polydispersities for which crystallization is impossible, we find a purely dynamic liquid-glass transition at a density $`\rho ^G`$, which depends very little on the degree of polydispersity $`\epsilon `$. For moderate values of $`\epsilon `$, the glass transition density can be defined even better than in the highly polydisperse case considered earlier . Interestingly, $`\rho ^G`$ is considerably higher than the crystallization density for the monodisperse system.

Our vastly improved equilibrium simulation methods allow us to scan the complete parameter space ($`\epsilon ,\rho `$). Even far inside the glassy phase, we can compute the equation of state and the compressibility with very high precision. The polydispersity-driven transition of the crystal into the disordered phase (liquid or glass) also turns out to be particularly interesting: we find strong evidence that this disorder-driven phase transition is continuous. Our findings can hardly be reconciled with the presence of a first-order transition, and differ from what was stated in a closely related work .

In principle, for a polydisperse mixture of particles, one should choose the particle sizes $`r_i`$ independently from a given probability distribution $`p(r)`$. We would then have to perform an ensemble average over the distribution of radii. In dynamical simulations, this average corresponds to ensemble-averaging time correlation functions of a given sample. In thermodynamic calculations, the average over radii can be incorporated directly into the MC evaluation of the partition function, as some authors have done .

We have adapted a much simpler, yet practically equivalent approach, by considering a ‘fixed probability increment’ model for any probability distribution $`p(r)`$ with $`r>0`$, where we pick the $`N`$ radii $`r_i`$ e.g. according to the following equation

$$_0^{r_i}𝑑rp(r)=\frac{i}{N+1}$$
(1)
(with $`_0^{\mathrm{}}𝑑rp(r)=1`$). In this way, sample-to-sample fluctuations are completely eliminated, and the most representative sample is generated for an arbitrary distribution $`p(r)`$.

Here, we study a model with a flat distribution of radii:

$$p(r)=\{\begin{array}{cc}\text{const}\hfill & r_{min}rr_{max}\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}$$
(2)
i.e. the radii are given by $`r_i=r_1+\mathrm{\Delta }(i1)`$ (with $`\mathrm{\Delta }=(r_nr_1)/(N1)`$). As in , we set the total particle volume as $`V_{part}=N\pi /4`$ and the Boltzmann factor as $`\beta =1`$.

We define the polydispersity by the normalized width of the distribution eq.(2)

$$\epsilon =\sqrt{<\delta r^2>/<r>^2}.$$
(3)
Notice that the distribution eq.(2) satisfies $`\epsilon 1/\sqrt{3}`$, this limit being attained for $`r_{min}=0`$.

Our simulations are performed in a periodic box of size $`L\times \sqrt{3}L/2`$, as is customary , and the only difference in set-up between our thermodynamic and dynamic simulations is the choice of ensembles. The equation of state $`V(P)`$ is most naturally computed in the \[NP\] ensemble, where the size of the box is allowed to fluctuate. The dynamical calculations are done in the \[NV\] ensemble, where the time evolution of states can be best monitored. For these calculations, a conventional local MC algorithm is used. Up to a trivial rescaling of time, we expect this algorithm in the long-time limit to give the same results as a molecular dynamics method.

The equilibration of polydisperse systems is a difficult task, because the properties of the system depend sensitively on the relative arrangements of particles, which have to be averaged over. Therefore, one needs an algorithm that moves particles over large distances and averages over different possible arrangements. In the fluid phase this can in principle be achieved by a long series of local moves. In the crystal or the glass, it is however mandatory to equilibrate the system by non-local moves, as in the pivot cluster algorithm .

A polydisperse system with a continuous distribution of radii is easier to simulate than a binary distribution, because it is often possible to exchange particles with similar radii. These swaps also implement non-trivial long-range moves, and have to be combined with a local Monte Carlo algorithm. The acceptance probability for such swaps decreases only at the most extreme densities, when each particle is so close to its neighbors that it can virtually never be replaced by a particle with slightly larger radius.

In high-density binary mixtures (cf ), the pivot cluster algorithm is clearly more appropriate, because the direct particle exchanges freeze out. In the present polydisperse case, however, we have obtained identical results with both algorithms for systems of $`N=256`$ and $`N=1024`$ particles even slightly below what we believe to be close packing.

Extensive tests of ergodicity (cf ) fully confirm the above statements: In very small systems ($`N=15`$), where the minimal size differences between particles are larger, we can show that the direct swap algorithm falls out of equilibrium at slightly lower densities than the pivot cluster method for the above-mentioned reasons.

We point out a completely unrelated, yet crucial, issue: To compute the equation of state at constant pressure, we need to include volume changes into the MC algorithm. Extremely time-consuming sampling of the volume by trial and error (using the Metropolis algorithm ) can be avoided .

For orientation, we show in figure 1 the schematic phase diagram as a function of density $`\rho `$ and polydispersity $`\epsilon `$ (cf eq.(3)) which results from our thermodynamic and dynamic calculation, and which we discuss in the remainder of this paper. We have studied this diagram by the routes indicated by gray lines in the figure, using either the non-local MC algorithm (for thermodynamic averages), or the local MC method (to study the dynamics).

In Figure 2a, we show the liquid and solid branches of the thermodynamic equation of state as a function of pressure $`P`$ for a fixed $`\epsilon =0.052`$. The presence of a thermodynamic transition is evident.

In these calculations at $`\epsilon =0.052`$, and at $`N=256`$, we are able to establish coexistence of the two phases as in the monodisperse system . Even there however, the $`N\mathrm{}`$ limit has remained controversial, and it is not yet firmly established whether the liquid-solid transition in 2-d monodisperse disks is first-order or continuous (according to the KTNHY scenario, cf ). Notice that the transition shows up very clearly both in the equation of state and in the compressibility, shown in the main figure 2a.

In figure 2b, we show analogous data for $`\epsilon =0.19`$. There, we see no indication of a thermodynamic phase transition in the whole range of densities studied. Similar results were obtained at $`\epsilon =0.52`$ (remember that $`\epsilon <0.58`$ cf text following eq.(3)); these data confirm and considerably extend our earlier calculations . Within our flat probability distribution of the radii of particles, there exists a maximal polydispersity $`\epsilon =0.58`$, which corresponds to the limit $`r_{min}0`$ in eq.(2).

Our calculations thus confirm the existence of a terminal polydispersity, which was studied mostly in three dimensions, but also evoked in related two-dimensional work .

The horizontal sweeps in the phase diagram of figure 1 bracket the disordered-to-crystalline phase boundary, which must be essentially independent of $`\epsilon `$ at large $`\rho `$. We further studied the thermodynamic system as a function of polydispersity at constant pressure, along the vertical paths indicated in figure 1. There, our data strongly favor a continuous transition: We did not find hysteresis in $`V(P)`$, and detected no abrupt change of the volume distribution function in the transition region. We also computed the average bond-orientational order-parameter $`<\mathrm{\Psi }_6>`$ (cf, e.g.,  for definitions and a discussion). There also, the probability distribution of $`\mathrm{\Psi }_6`$ has a single peak, and is perfectly reproduced, without any indications of hysteresis. We conclude that the transition is continuous. In figure 3 we show $`<\mathrm{\Psi }_6>`$ as a function of polydispersity $`\epsilon `$ for different system sizes. A detailed, more rigorous study of the probable two-step melting process in this system is left for further work. Let us note that the numerical analysis of the standard Kosterlitz-Thouless transition has proven to be very subtle even in the prototypical XY-model .

Finally, we determine the dynamic properties of the system, in the \[NV\] ensemble, with the local MC algorithm, and at values of $`\epsilon `$ which ensure the absence of a liquid-solid transition. Using the protocol of ref. , we obtained the effective diffusion constants $`D_i`$ (cf figure 4), for $`\epsilon =0.19`$ and $`\epsilon =0.52`$. In both cases, the results are consistent with the scaling form :

$$D_i(\rho )(\rho _i^G\rho )^\alpha ,$$
(4)
where we extrapolate the diffusion constants for each particle separately. The diffusion constants $`D_i`$ strongly depend on $`i`$, but the extrapolated transition densities $`\rho _i^G`$ do not, especially for small $`\epsilon `$. We take these extrapolated values as a definition of the glass transition density $`\rho ^G`$. Only at large $`\epsilon `$ (left part of figure) can some very small particles ($`r_{min}/r_{max}=19`$) escape through slits left open in the system. This leads to some additional structure in the $`D_i`$ plot on the left half of figure 4, and to an increased value of the extrapolated $`\rho _i^G`$ for small $`i`$.

In conclusion, our study of the polydisperse hard disk system has given the homogeneous phases shown in the schematic diagram figure 1. Inhomogeneous phases  seem to play no role at the parameters studied in this paper. A tendency towards phase separation should first show up in irregular finite-size behavior, which we have not encountered. We find it particularly significant that the glass transition depends so little on polydispersity, and takes place at densities (and pressures) much higher than the crystallization density in the monodisperse system. Using our specialized MC methods, we are able to probe the thermodynamics of the hard-disk system far within the glassy phase, where we find no indications of an accompanying thermodynamic transition.

Acknowledgements: We thank P. Le Doussal for helpful discussions. L. S. acknowledges support from the Deutsche Forschungsgemeinschaft under Grant No. SA864/1-2.