text:
# 1 Introduction

## 1 Introduction

Magnetic fields are ubiquitous in the Universe, being observed in nature on scales from planetary size to galaxy cluster size . In galaxies and galaxy clusters, the typical strength is of order a few $`\mu `$Gauss, which is thought to be produced by dynamo action on a seed field. In galaxies the dynamo timescale is roughly a rotation period, $`10^8`$ Yr, and a simple calculation  based on the age of a typical galaxy shows that the seed field must have been about $`10^{20}`$ Gauss, or perhaps less in the currently favoured models with a cosmological term .

There is no shortage of ideas for generating this seed field. The more conventional astrophysical explanations are based on a Biermann battery operating at the era of reionisation (see e.g.  and the references therein). There are more speculative ideas based on various models of inflation , phase transitions  and primordial black holes .

All these mechanisms have the common feature of producing stochastic, homogeneous and isotropic magnetic and velocity fields which can be characterised by their power spectra and characteristic initial scales. Our interest here is to try and make model-independent statements about the evolution of the magnetic fields once they are generated. This article, which is based on Ref. , studies the evolution of a stochastic magnetic field generated at a phase transition, such as the confinement transition in QCD at $`t1`$ sec, or the electroweak symmetry-breaking transition at $`t10^{11}`$ sec. It therefore falls into the category of decaying 3D MHD turbulence, which has been studied before in the MHD community . Most directly comparable to our work, Biskamp and Müller  studied the energy decay in incompressible 3D magnetohydrodynamic turbulence in numerical simulations at relatively high Reynolds number, and in a companion letter  studied the scaling properties of the energy power spectrum.

We focus here on the transfer of magnetic energy from small to large scales, as necessitated by the conservation of magnetic helicity. This is important for a primordial magnetic field to reach a large enough scale with sufficient amplitude to be relevant for seeding the galactic dynamo .

We perform 3D simulations both with and without magnetic helicity, starting from statistically homogeneous and isotropic random initial conditions, with power spectra suggested by cosmological applications. We find a strong inverse cascade in the helical case, with the coherence scale of the field growing as $`t^{0.5}`$, and equivocal evidence for a weak cascade when only helicity fluctuations are present. In the helical case we also find a self-similar power spectrum with an approximately $`k^{2.5}`$ behaviour at high $`k`$. We find decay laws for the magnetic and kinetic energies of $`t^{0.7}`$ and $`t^{1.1}`$, respectively, in the helical case, and $`t^{1.1}`$ for both in the non-helical case.

## 2 Evolution of magnetic fields in the early Universe

A convenient benchmark is to assume that the field is created on the horizon scale with a power spectrum $`k^n`$, taking all the energy in the Universe, and that it is subsequently completely frozen into the plasma . If one takes the epoch of creation to be either the QCD transition or the electroweak transition, one finds that the RMS fluctuations on the scale of a protogalaxy are roughly

$$B(t_0,L=0.5\mathrm{Mpc})<\{\begin{array}{cc}10^{143n}\text{G}\hfill & \mathrm{QCD}\hfill \\ 10^{204.5n}\text{G}\hfill & \mathrm{EW}\hfill \end{array}$$

Thus we see that for causal fields ($`n2`$) there needs to be some amplification on large scales for there to be any interesting seed field for the galactic dynamo.

In fact, the field is not strictly frozen into the plasma, and in order to calculate observable effects, we must determine how the field scale length $`\xi `$ and magnetic energy $`E_M`$ evolves. There are various scaling arguments which have been put forward: for example, for ideal MHD (infinite conductivity) and with no helicity Olesen  (and later Son , Field and Carroll  and Shiromizu ) argued

$$\xi (t)t^{2/(n+5)},$$
(1)
where $`t`$ is conformal time. The effect of having a conserved helicity modifies this scaling law to

$$\xi (t)t^{2/3},E_Mt^{2/3}.$$
(2)
Early numerical experiments with a shell model of the full MHD equations  suggested $`\xi t^{0.25}`$, and focussed attention on the possibility of an inverse cascade, in which power is transferred locally in $`k`$-space from small to large scales.

## 3 MHD equations

The MHD equations in an expanding Universe are most conveniently expressed in terms of conformally rescaled fields $`𝐁`$, $`𝐮`$ and dissipation parameters $`\nu `$, $`\eta `$, and in the gauge $`A^0=\eta \mathbf{}\mathbf{}𝐮`$ .

The matter and radiation in the early Universe is modelled as an isothermal compressible gas with a magnetic field, which is governed by the momentum equation, the continuity equation, and the induction equation, written here in the form

$`{\displaystyle \frac{𝐮}{t}}`$ $`=`$ $`𝐮\mathbf{}𝐮c_s^2\mathbf{}\mathrm{ln}\rho +{\displaystyle \frac{𝐉\times 𝐁}{\rho }}+{\displaystyle \frac{\mu }{\rho }}\left(^2𝐮+{\displaystyle \frac{1}{3}}\mathbf{}\mathbf{}𝐮\right),`$ (3)
$`{\displaystyle \frac{\mathrm{ln}\rho }{t}}`$ $`=`$ $`𝐮\mathbf{}\mathrm{ln}\rho \mathbf{}𝐮,`$ (4)
$`{\displaystyle \frac{𝐀}{t}}`$ $`=`$ $`𝐮\times 𝐁+\eta ^2𝐀,`$ (5)
where $`𝐁=\times 𝐀`$ is the magnetic field in terms of the magnetic vector potential $`𝐀`$, $`𝐮`$ is the velocity, $`𝐉`$ is the current density, $`\rho `$ is the density, $`\mu `$ is the dynamical viscosity, and $`\eta `$ is the magnetic diffusivity.

In the ideal limit $`\mu =\eta =0`$, there is a conserved quantity in addition to the energy, which is the magnetic helicity $`H_\mathrm{M}`$, given by

$$H_\mathrm{M}=𝐀𝐁d^3x.$$
(6)
Helicity is known to be important in dynamo theory , where turbulence is driven. We shall also be able to confirm its importance in decaying turbulence.

An important dynamical quantity is the magnetic Reynolds number $`Re_\mathrm{M}=Lv/\eta `$, where $`L`$ and $`v`$ are the typical length scale and velocity of the system under consideration, because it measures the importance of the non-linear term in the equation for the magnetic field. In the early Universe, $`Re_\mathrm{M}`$ can be very large. This is often taken to mean that the magnetic field is frozen into the plasma, and the scale length of the field increases only with the expansion of the Universe. This is in general untrue, because turbulence can transfer energy to different length scales . Not only can there be the usual direct cascade of energy from large to small scales, but also an inverse cascade, increasing the overall comoving correlation length .

## 4 Plasma properties in the early Universe

We are considering the time before recombination, when the plasma in the early Universe has a relativistic equation of state $`p=\rho /3`$, and therefore a sound speed $`c_s=1/\sqrt{3}`$. The number species contributing to the pressure gradually decreases with the temperature $`T`$, until $`Tm_e`$, the electron mass, when only photons and neutrinos remain relativistic. At that point the number of charge carriers in the plasma reduces by a very large factor: for $`Tm_e`$, the electron number density $`n_e`$ is approximately equal to the photon number density $`n_\gamma `$, while for $`Tm_e`$ $`n_e/n_\gamma =n_B/n_\gamma 10^{10}`$, where $`n_B`$ is the baryon number density.

The transport properties of the plasma are determined from the mean free path $`l_{\mathrm{mfp}}`$ of the relevant particles involved in the transport of the quantity of interest, which typically is $`l_{\mathrm{mfp}}1/\alpha |\mathrm{log}\alpha |T`$ , where $`\alpha `$ is the fine structure constant. From this we can infer that the conductivity is

$$\sigma \{\begin{array}{cc}T/\alpha |\mathrm{log}\alpha |\hfill & Tm_e\hfill \\ (n_e/n_\gamma )(m_e^2/T)/\alpha \hfill & Tm_e\hfill \\ \chi (n_B/n_\gamma )(T/m_e)^2Te^2\hfill & T<T_{\mathrm{dec}}\hfill \end{array},$$
(7)
where $`T_{\mathrm{dec}}`$ is the temperature at photon decoupling, roughly 1 eV.

The viscosity parameter $`\nu `$ is given by

$$\nu \eta /\rho \{\begin{array}{cc}1/\alpha ^2|\mathrm{log}\alpha |T\hfill & Tm_e\hfill \\ (n_\gamma /n_e)(m_e^2/T^2)/\alpha T\hfill & Tm_e\hfill \end{array}$$
(8)
We can estimate the importance of non-linear terms relative to diffusion from the Reynolds numbers of the plasma, the hydrodynamical $`\mathrm{Re}=vL/\nu `$ and magnetic $`\mathrm{Re}_M=vL\sigma `$. Upper bounds are obtained by assumed a fluid flowing at the speed of light on a scale somewhere between the mean free path and the particle horizon $`ct`$:

$$(\alpha ,1/\alpha ^2\mathrm{log}^2\alpha )\begin{array}{c}<\hfill \\ \hfill \end{array}(\mathrm{Re},\mathrm{Re}_M)\begin{array}{c}<\hfill \\ \hfill \end{array}10^{19}(T/\mathrm{GeV})^1$$
(9)
## 5 3D MHD simulations of decaying turbulence

We solve these equations numerically with a code  using a variable third order Runge-Kutta timestep and sixth order explicit centred derivatives in space. All our runs are performed on a $`120^3`$ grid, and we use periodic boundary conditions, which means that the average plasma density $`\rho _0=\rho _0`$ is conserved during runs. Here $`\rho _0`$ is the value of the initially uniform density, and the brackets denote volume average.

We use natural units $`c=1`$, and set the unit of length by setting measure $`k_1=1`$, where $`k_1`$ is the smallest wave number in the simulation box. Hence the box has a size of $`L_{\mathrm{BOX}}=2\pi `$. The scale factor is fixed by setting $`\rho _0=1`$, and $`𝐁`$ is measured in units of $`\sqrt{\mu _0\rho _0}c`$, where $`\mu _0`$ is the magnetic permeability. We define the mean kinematic viscosity $`\nu \mu /\rho _0`$. The sound speed $`c_s=1/\sqrt{3}`$, as appropriate for a relativistic fluid.

The equations are not quite those for a relativistic gas in the early universe . However, we have checked that our results change little when using the true relativistic equations in the low velocity limit.

The initial conditions appropriate for the early Universe are to take both the velocity and magnetic fields to be homogeneous and isotropic Gaussian random fields drawn from a power-law distribution with a high wavenumber cut-off, determined by the physical mechanism generating the power spectrum. The power spectra, defined by $`P_\mathrm{M}(k)𝐁_𝐤^{}𝐁_𝐤`$, and $`P_\mathrm{V}(k)𝐮_𝐤^{}𝐮_𝐤`$ is taken to have the form

$$P_\mathrm{M}(k)=A_Mk^n\mathrm{exp}((k/k_c)^4),P_\mathrm{V}(k)=A_Vk^m\mathrm{exp}((k/k_c)^4).$$
(10)
note that in the plots it is the shell-integrated energy spectra, $`E_{\mathrm{M},\mathrm{V}}=4\pi k^2\times \frac{1}{2}P_{\mathrm{M},\mathrm{V}}(k)`$, which are shown.

Note also that causality demands that $`n2`$ and $`m0`$ . In the simulations presented we took the power laws to be the lowest consistent with causality. We also chose $`k_c=30`$, unless specified otherwise. The initial magnetic energy was taken equal to the kinetic energy, and had the value $`5\times 10^3`$ in all runs.

We are also able to control the helicity in the initial conditions, ranging from identically zero to maximal, where maximal helicity is defined as saturating the inequality

$$|H_\mathrm{M}(k)|2k^1E_\mathrm{M}(k).$$
(11)
It should be noted that if one chooses the initial field with no constraint, there are helicity fluctuations present, around a mean of zero. It is possible to completely eliminate the fluctuations but in practice they are quickly regenerated by the evolution of the field.

## 6 Results

In all runs the mean kinematic viscosity $`\nu `$ and the resistivity $`\eta `$ were chosen to be equal with values between $`\nu =\eta =5\times 10^45\times 10^5`$.

### 6.1 Magnetic energy spectrum

In Fig. 1 we show magnetic energy spectra $`E_\mathrm{M}(k)`$ for runs with unconstrained and maximal initial magnetic helicity.

Turning first to maximal initial helicity, one sees that after a short initial direct cascade, there is energy transfer towards lower $`k`$, and the peak of the power spectrum moves to the left, signifying a growth in the coherence scale of the field. This has been called an inverse cascade, although it should be noted that the term “cascade” carries the implication of a local interaction in $`k`$-space, which is not necessarily true . In fact, in simulations where the initial conditions contain power in only a single wavenumber, the $`k^4`$ spectrum is very quickly established at low $`k`$ , indicating that the power transfer to large scales is non-local in $`k`$-space .

For wavenumbers above the peak, one sees a decrease in power with quite a different form than is present in the initial conditions, which can be fitted approximately by $`E_\mathrm{M}(k)k^{2.5}`$.

When only helicity fluctuations are present in the initial conditions, the peak in the power spectrum still moves to low $`k`$, but there is very little transfer of power. It is perhaps surprising to see any transfer at all, as there are arguments which connect the evolution of the coherence length with helicity conservation . In an attempt to test this idea we also did a run with identical parameters but zero initial helicity, but saw essentially identical behaviour, with a small increase in power at small scales.

In our runs with $`E_\mathrm{K}=E_\mathrm{M}`$ initially, the kinetic energy spectrum shows no evidence of an inverse cascade at any scale. However, when the initial velocity distribution is zero the kinetic spectrum grows on all scales initially and in the low wave number region the energy continue to grow even after the high wave number modes start to decay.

### 6.2 Coherence length evolution

One length scale is the magnetic Taylor microscale $`L_\mathrm{T}`$, defined as $`B_{\mathrm{rms}}/J_{\mathrm{rms}}`$. In Fig. 2 we show $`L_\mathrm{T}`$ against time $`t`$ for a run with maximal initial helicity.

The asymptotic behaviour of the length scale is seen to grow approximately as $`L_\mathrm{T}t^{0.5}`$.

In runs with non-helical initial conditions the growth of the magnetic Taylor microscale is slower: we find approximately $`L_\mathrm{T}t^{0.4}`$.

### 6.3 Energy decay laws

Figure 3 shows the magnetic energy $`E_\mathrm{M}(t)`$ and the kinetic energy $`E_\mathrm{K}(t)`$ with maximal initial helicity.

It is seen that the asymptotic decay rate for $`E_\mathrm{M}(t)`$ is approximately $`t^{0.7}`$. The kinetic energy also decays with a power law behaviour at late times: $`E_\mathrm{K}(t)t^{1.1}`$.

In runs without initial helicity the decay rates of $`E_\mathrm{M}(t)`$ and $`E_\mathrm{K}(t)`$ are approximately the same, close to $`t^{1.1}`$.

### 6.4 Magnetic Reynolds number

The Reynolds numbers in our simulations are evaluated using the magnetic Taylor microscale. In most our simulations we typically obtain Reynolds numbers of the order of $`100200`$. In Fig. 4 we show the magnetic Reynolds number of the run with maximal helicity whose power spectrum is shown in Fig. 1.

Note that for the second half of the run the magnetic Reynolds number is approximately constant.

### 6.5 Self-similarity in magnetic power spectrum

We make the following ansatz for the energy spectrum

$$E_\mathrm{M}(k,t)=\xi (t)^qg_\mathrm{M}(k\xi ),$$
(12)
where $`\xi `$ is the characteristic length scale of the magnetic field, taken to be the magnetic Taylor microscale defined above, and $`q`$ is a parameter whose value is some real number.

Figure 5 shows $`\xi (t)^qE_\mathrm{M}(k,t)`$ versus the scaled variable $`k\xi (t)`$, with $`q=0.7`$, plotted for different values of time $`t`$. The data is seen to collapse onto a single curve given $`g_\mathrm{M}(k\xi )`$: thus the magnetic field evolves in a self-similar manner.

## 7 Discussion and conclusions

The most directly comparable simulations of decaying 3D MHD turbulence were carried out by Biskamp and Müller . They found similar results, the magnetic field evolved self-similarly, with a power-law behaviour at high $`k`$. However, their power law was $`k^{5/3}`$, much less steep than our $`k^{2.5}`$. There were a number of differences between their and our simulations: they were able to achieve larger Reynolds numbers, both by having larger grids, and by using hyper-diffusivity (a $`^4`$ magnetic diffusivity term). However, we believe that the real difference is due to their initial cut-off scale being significantly larger, at $`k_c=4`$. We have performed a run with larger initial length scale, $`k_c=5`$. In this case the magnetic energy spectrum develop into an approximate $`k^{5/3}`$ law at late times. However, this occurs only after the peak of the spectrum has left the simulation box.

There have also been renormalisation group (RG) analyses looking for an inverse cascade in driven MHD turbulence . In particular Berera and Hochberg  saw no evidence for an inverse cascade. However, it is not clear that the results are directly comparable, firstly because we are considering freely decaying turbulence, and secondly because RG analysis can only give information about late times when the system is in equilibrium with the driving force. What we have referred to as an inverse cascade is, even in the driven case , a time-dependent phenomenon characterised by a bump in the power spectrum travelling to smaller wave number.

In conclusion, we find good evidence from our numerical simulations that helical stochastic magnetic fields show an inverse cascade (in the sense explained above), and that even if only small helicity fluctuations are present initially, there is still weak inverse cascade.

We have determined growth laws for the magnetic and kinetic energies $`E_M`$ and $`E_K`$. In the helical case, $`E_Mt^{0.7}`$ and $`E_Kt^{1.1}`$, which means that that there is no equipartition of energy. This is because the extra constraint of helicity conservation forces the magnetic field to transfer power to larger scales rather than allow it to be dissipated. The importance of helicity is borne out by the fact that in the non-helical case, we find $`\xi t^{0.4}`$ and $`E_ME_Kt^{1.1}`$.

Length scales in the magnetic field increase as $`t^{0.5}`$ in the helical case, but slightly slower in the non-helical case, $`t^{0.4}`$. Note that these growth laws disagree with all theoretical predictions to date , which give $`t^{2/3}`$ in the helical case and $`t^{2/7}`$ for our power spectrum in the non-helical case.

Helical magnetic fields are found to evolve in a self-similar way, with a scaling function $`g_M(z)z^p`$ at large $`k`$, where $`p=2.5`$ for $`\mathrm{Re}10^2`$. Note that this is significantly different from both the Iroshnikov-Kraichnan and Kolmogorov spectra, $`k^{3/2}`$ and $`k^{5/3}`$, respectively.

A good theoretical understanding of these scaling laws is required before the evolution of magnetic fields in the early Universe is properly understood, as a small error in the exponent makes a large error in the prediction of the magnetic field strength when propagated over many orders of magnitude in time. For example, Vachaspati’s contribution to these proceedings  assumes a growth law of $`t^{2/3}`$ in the length scale, based on a simple argument invoking helicity conservation , to obtain seed fields of an interesting strength from the electroweak transition, which is quite different from our observed growth law of $`t^{0.5}`$.

ACKNOWLEDGEMENTS

This work was conducted on the Cray T3E and SGI Origin platforms using COSMOS Consortium facilities, funded by HEFCE, PPARC and SGI. We also acknowledge computing support from the Sussex High Performance Computing Initiative. MH thanks NORDITA for hospitality while this work was completed.


format: text

text:
# Anomalous Dispersion of Longitudinal Optical Phonons in 𝐍𝐝_1.86𝐂𝐞_0.14𝐂𝐮𝐎_{𝟒+𝜹} Determined by Inelastic X-ray Scattering

## Abstract

The phonon dispersions of $`\mathrm{Nd}_{1.86}\mathrm{Ce}_{0.14}\mathrm{CuO}_{4+\delta }`$ along the $`[\xi ,0,0]`$ direction have been determined by inelastic x-ray scattering. Compared to the undoped parent compound, the two highest longitudinal phonon branches, associated with the Cu-O bond-stretching and out-of-plane oxygen vibration, are shifted to lower energies. Moreover, an anomalous softening of the bond-stretching band is observed around $`𝐪=(0.2,0,0)`$. These signatures provide evidence for strong electron-phonon coupling in this electron-doped high-temperature superconductor.

While the coupling between electrons and phonons is known to be the driving mechanism for Cooper-pair formation in conventional superconductors, its role in the high-critical-temperature superconductors (HTcS) is the subject of intense research efforts. Recently, evidence for electron-phonon coupling has been invoked in the interpretation of inelastic neutron scattering (INS) and angle-resolved photoemission spectroscopy (ARPES) experiments. The INS studies, carried out on $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$ Pintschovius et al. (1991); Pintschovius and Reichardt (1995); McQueeney et al. (1999); Pintschovius and Braden (1999), oxygen-doped $`\mathrm{La}_2\mathrm{CuO}_{4+\delta }`$ Pintschovius and Braden (1996) and $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_{6+\delta }`$ Pintschovius and Reichardt (1995); Reichardt (1996) reveal an anomalous softening with doping of the highest longitudinal optical (LO) phonon branch, in particular along the $`q=[\xi ,0,0]`$ direction. This branch is assigned to the Cu-O bond-stretching mode Zhang et al. (1991); Pintschovius and Reichardt (1995). The observed softening has been interpreted as a signature of a strong electron-phonon coupling Pintschovius and Reichardt (1995); McQueeney et al. (1999), which has been discussed since the discovery of HTcS Bednorz and Müller (1986); Weber (1987). Furthermore, in an energy range similar to the LO bond-stretching phonon band, ARPES studies on three different families of hole-doped HTcS reveal a distinct “kink” anomaly in the quasiparticle dispersion Lanzara et al. (2001). The scenario emerging from the above INS and ARPES works suggests a strong coupling between the charge carriers and the Cu-O bond-stretching phonon modes to be ubiquitous in HTcS materials, at least for hole-doped compounds. However, its role in the pairing mechanism remains completely unclear Pintschovius and Braden (1999). At the moment, it may not even be excluded that the electron-phonon interaction is pair-breaking for the d-wave superconducting order parameter. Therefore, it appears very important to analyze the strength of the phonon anomalies in as many cuprate families as possible and to compare them with their superconducting properties. In this context, the electron-doped cuprates are of central importance due to the distinct character of the doped charges in this material: Cu $`3d_{x^2y^2}`$ (O $`2p`$) for n(p)-type cuprates Y.Tokura et al. (1989), leading to a very different electronic structure Armitage et al. (2001). Since the phonon anomalies are related to a coupling between charge fluctuations and phonons, charges with different character may induce quite different electron-phonon interactions.

In this Letter, we present an inelastic x-ray scattering (IXS) study of the phonon dispersion in the n-type cuprates $`\mathrm{Nd}_{1.86}\mathrm{Ce}_{0.14}\mathrm{CuO}_{4+\delta }`$ (NCCO). Inelastic x-ray scattering can overcome the main limitation of inelastic neutron scattering, i.e. the need for sufficiently large single crystals of high chemical and structural quality Pintschovius and Reichardt (1995). Lateral x-ray beam sizes of few tens of $`\mu `$m are routinely obtained. Moreover, at photon energies around 10-20 keV and $`Z>`$ 3, the total cross section is dominated by photoelectric absorption, and therefore the typical x-ray penetration depths for high-$`Z`$ materials is of the order of 10 - 100 $`\mu `$m. Consequently, very small samples (down to less than $`10^4`$ mm3) can be studied with signal rates comparable to typical INS experiments on cm3-sized samples. Despite these advantages, little work has been done using IXS on the HTcS compounds Burkel (1991). We choose $`\mathrm{Nd}_{2\mathrm{x}}\mathrm{Ce}_\mathrm{x}\mathrm{CuO}_{4+\delta }`$ for our IXS study, since its crystallographic structure is one of the simplest among the HTcS, and because extensive INS studies exist for its undoped parent compound $`\mathrm{Nd}_2\mathrm{CuO}_{4+\delta }`$ (NCO) Pintschovius et al. (1991); Pintschovius and Reichardt (1995). Our interest is focused on the $`[\xi ,0,0]`$ direction, where the LO branch displays its strongest anomaly for hole-doped HTcS McQueeney et al. (1999); Pintschovius and Braden (1999). The present results reveal that, near the zone center, the two highest longitudinal optical branches, assigned to the Cu-O bond-stretching and O(2) vibration modes, are shifted to lower frequencies with respect to the undoped parent compound. The interpretation of our data is supported by lattice dynamics calculations, taking into account a Thomas-Fermi screening mechanism. Furthermore, we observe an anomalous softening of the highest branch around $`𝐪=(0.2,0,0)`$. Our results demonstrate that this anomalous behavior of the high-energy LO phonon branch is a universal property of both hole- and electron-doped HTcS compounds, therefore providing further evidence that electron-phonon interactions may play an important role in high-Tc superconductivity. Furthermore, the present results are an important demonstration of IXS as a powerful tool for the study of the lattice dynamics in small, high-quality crystals of complex transition metal oxides.

The experiment was carried out at the very-high-energy-resolution IXS beam-line ID16 at the European Synchrotron Radiation Facility (ESRF). X-rays from an undulator source were monochromated using a Si (111) double-crystal monochromator, followed by a high-energy-resolution backscattering monochromator Verbeni et al. (1996), operating at 15816 eV (Si (888) reflection order). A toroidal gold-coated mirror refocused the x-ray beam onto the sample, where a beam size of $`250\times 250\mu `$m2 full-width-half-maximum (FWHM) was obtained. The scattered photons were energy-analyzed by a spherical silicon crystal analyzer 3 m in radius, operating at the same Bragg reflection as the monochromator Masciovecchio et al. (1996). The total energy resolution was 1.6 THz (6.6 meV) FWHM. The momentum transfer $`𝐐`$ was selected by rotating the 3 m spectrometer arm in the scattering plane perpendicular to the linear x-ray polarization vector of the incident beam. The momentum resolution was set to $`0.087\mathrm{\AA }^1`$ in both the horizontal and the vertical direction by an aperture of $`20\times 20`$ mm2 in front of the analyzer. Further experimental details are given elsewhere (Verbeni et al. (1996); Masciovecchio et al. (1996) and references therein). The sample is a single crystal grown by the traveling-solvent floating-zone method in 4 atm of $`\mathrm{O}_2`$ at Stanford University. It has been reduced under pure Ar atmosphere at $`920^{}`$C for 20 h, followed by a further 20 h of exposure at $`500^{}`$C to a pure $`\mathrm{O}_2`$ atmosphere. Such a procedure is necessary to produce a superconducting phase in $`\mathrm{Nd}_{1.86}\mathrm{Ce}_{0.14}\mathrm{CuO}_{4+\delta }`$, although its exact effect is not understood. Following this treatment the sample had a narrow superconducting transition with an onset temperature of $`T_c=24`$ K. The sample is of very good crystalline quality, with a rocking curve width of 0.02 (FWHM) around the $`[h,0,0]`$ direction. It was mounted on the cold finger of a closed-loop helium cryostat, and cooled to 15 K. The experiment was performed in reflection geometry, and the probed scattering volume corresponded to about $`1.5\times 10^3`$ mm3. IXS scans were performed in the -2$`<\nu <`$24 THz range, in the $`𝝉=(6,0,0)`$ and $`𝝉=(7,1,0)`$ Brillouin zones 111All reciprocal lattice vectors are expressed in r.l.u.. The data were collected along the following three lines:
I) $`𝐐=(6+\xi ,0,0)`$, in longitudinal configuration (i.e. with $`𝐪=(\xi ,0,0)𝐐`$),
II) $`𝐐=(7\pm \xi ,1,0)`$ in almost longitudinal configuration (i.e. with $`𝐪=(\xi ,0,0)`$ and $`(𝐐𝐪)/Qq`$),
III) $`𝐐=(7,1\xi ,0)`$ in almost transverse configuration ($`𝐪=(0,\xi ,0)`$ and $`(𝐐𝐪)/Q0`$).
The low temperature and high momentum transfer were chosen so as to optimize the count rate on the high-frequency optical mode while limiting the loss of contrast due to the contribution from the tails of the intense low frequency acoustic modes.

Fig. 1 shows a typical energy scan in almost longitudinal geometry at $`𝐐=(6.8,1,0)`$, corresponding to $`\xi =0.2`$. The data are shown together with the results of a fit, where the excitations where modeled by harmonic oscillators, convoluted with the instrumental resolution function. Three features can be clearly distinguished: (i) an elastic peak at 0 THz, due to chemical disorder and, possibly, strain due to the different thermal expansion coefficients between the sample and the sample support during cooling, (ii) the longitudinal acoustic phonon, centered near 2.7 THz and (iii) a weaker feature around 13 THz. In the intermediate energy region (from about 4 to 10 THz), no distinct phonon peaks are resolved due to the dominating contribution from the tails of the elastic and acoustic phonon signals. The inset of Fig. 1 emphasizes the high-energy region of the spectrum. One can clearly distinguish two phonons, centered around 12 and 13.5 THz, respectively. The weak shoulder at around 16 THz can be attributed to an admixture of the transverse optical (TO) mode, as indicated by the comparison (after normalization to the same intensity at 20 THz) with the equivalent scan in transverse geometry $`𝐐=(7,0.8,0)`$. Consequently, the two stronger peaks are unambiguously assigned to longitudinal optical (LO) modes. From the absence of other higher frequency modes up to 23 THz, we conclude that the two modes at 12 and 13.5 THz can be identified with the O(2) vibration and Cu-O bond-stretching modes, respectively (see Refs. Zhang et al., 1991; Pintschovius and Reichardt, 1995).

In order to determine the dispersion of the two highest LO branches, IXS spectra were recorded for $`0<\xi 1`$. In Fig. 2 we show three spectra taken along the $`(7\xi ,1,0)`$ direction. Close to the zone center, at $`𝐪=(0.1,0,0)`$, the highest frequency mode is observed slightly above 15 THz. At $`𝐪=(0.2,0,0)`$ the highest mode is found at the much lower energy of 13.5 THz. Finally, for $`𝐪=(0.4,0,0)`$, we again find a high-frequency mode around 15.5 THz.

The peak positions extracted from these and many other scans are summarized in Fig. 3. The highest branch exhibits a sharp dip around $`𝐪=(0.2,0,0)`$ and recovers for larger $`q`$-values. This behavior is most likely due to an anti-crossing with the second highest branch which is mainly associated with vibrations of the O(2) position in the $`\xi `$-direction. Therefore, within a standard anti-crossing framework, one should interpret the highest longitudinal intensities for $`𝐪=(0.3,0,0)`$ and above as being mainly due to O(2) vibration. Within that scenario, this second-highest branch increases its frequency in the middle of the zone as in the undoped compound, and, except for the fact that the Lyddane-Sachs-Teller (LO-TO) gap closes, seems to be insensitive to doping. The LO bond-stretching mode just above $`𝐪=(0.2,0,0)`$ is then found at quite low energies, $``$ 12 THz, but can not be unambiguously followed to larger q-values. Nevertheless, our data document that the LO bond-stretching branch in NCCO is strongly renormalized compared to undoped NCO, in particular it bends down anomalously from the zone center to $`𝐪=(0.2,0,0)`$.

In order to further validate the correctness of our assignments, we performed a lattice dynamical calculation Mirone (2001) based on a shell model. We used a common potential model for cuprates, in which the interatomic potentials have been derived from a comparison of INS results for different HTcS compounds by Chaplot et al. Chaplot et al. (1995), using a screened Coulomb potential, in order to simulate the effect of the free carriers introduced by doping. Following Ref. Chaplot et al. (1995) for metallic $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$ and $`\mathrm{YBa}_2\mathrm{Cu}_3\mathrm{O}_7`$, we replaced the long-range Coulomb potential $`V_c(q)`$ by $`V_c(q)/ϵ(q)`$, and for the dielectric function we take the semi-classical Thomas-Fermi limit of $`ϵ(q)=1+\kappa _s^2/q^2`$, where $`\kappa _s^2`$ indicates the screening vector. The results of the calculation (without screening: dot-dashed lines; with screening: solid lines) are shown as well in Fig. 3. The lattice dynamics calculations without screening have been included, since they reproduce very well the experimental dispersion of the undoped parent compound Pintschovius et al. (1991). The shift at the zone center of the high-energy phonon branches of NCCO with respect to NCO is due to the closing of a large LO-TO splitting. Indeed, the corresponding $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_3`$ branches in $`\mathrm{Nd}_2\mathrm{CuO}_{4+\delta }`$ are separated by almost 3 THz at the zone center, as observed by INS Pintschovius et al. (1991). We point out that in our case a strong softening due to Thomas-Fermi screening does not imply a higher Thomas-Fermi parameter $`\kappa _s`$: actually, we find a screening vector $`\kappa _s`$ of about $`0.39\mathrm{\AA }^1`$, which is comparable to that for $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$ Chaplot et al. (1995). Though our modified calculation reproduces the closing of the LO-TO splitting, we still observe an anomalous additional softening of the highest bond-stretching LO branch near $`𝐪=(0.2,0,0)`$, which is not reproduced by our calculations (see Fig. 3). This branch softens in frequency from $`𝐪=(0.1,0,0)`$ to $`𝐪=(0.2,0,0)`$ by about $`\mathrm{\Delta }\nu `$ 1.5 THz, which is a shift comparable to the anomalous shift observed in $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$ at slightly larger $`𝐪`$ Pintschovius et al. (1991); McQueeney et al. (1999); Pintschovius and Braden (1999). Therefore, we believe that this anomalous softening is of the same nature as the one observed in p-type $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$.

A comparison between the experimental and calculated intensities for the two highest phonon branches is shown in Fig. 4. The good agreement of the observed integrated intensities with the calculated ones for the upper branches validates the correctness of our phonon branch assignment, at least for $`\xi 0.2`$. For $`\xi >0.2`$ we would have expected an intensity exchange between the two highest branches, which seems to be not observed.

The main difference between $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$ and $`\mathrm{Nd}_{1.86}\mathrm{Ce}_{0.14}\mathrm{CuO}_{4+\delta }`$ is, besides the screening effect, that in NCCO the Cu-O bond-stretching branch is closer in energy to the out-of-plane oxygen vibration one, having almost the same energy at $`\xi =0.2`$. These two branches belong to the same symmetry and therefore cannot cross, so that for $`\xi >0.2`$ softening implies interaction with the out-of-plane oxygen vibration mode with the same symmetry. In the region between $`𝐪=(0.25,0,0)`$ and $`(0.3,0,0)`$, the two modes are poorly defined in energy, which is consistent with what is observed in $`\mathrm{La}_{1.85}\mathrm{Sr}_{0.15}\mathrm{CuO}_{4+\delta }`$ by Pintschovius and Braden Pintschovius and Braden (1999) and McQueeney et al. McQueeney et al. (1999) for $`\xi 0.250.3`$. The corresponding real space periodicity of 3 to 4 unit cells may therefore be linked to the proximity to some charge instability. We remark that a reduced vector $`\xi 0.250.3`$ approximately corresponds to the nesting vector along $`[\xi 00]`$ direction, as can be inferred from the ARPES data of Ref. Armitage et al., 2001.

In conclusion, the present results reveal that the anomalous softening previously observed in hole-doped compounds Pintschovius et al. (1991); McQueeney et al. (1999); Pintschovius and Braden (1999, 1996); Pintschovius and Reichardt (1995); Reichardt (1996), is also present in the electron-doped cuprates. This is evidenced by the comparison of the present results on doped NCCO with the previously reported ones on pure NCO Pintschovius et al. (1991). This implies that: (i) the anomaly also exists in n-type cuprates, giving strength to the hypothesis Pintschovius et al. (1991); McQueeney et al. (1999); Pintschovius and Braden (1999, 1996); Pintschovius and Reichardt (1995); Reichardt (1996) of an electron-phonon coupling origin of this feature; (ii) this is a generic feature of the high-temperature superconductors, as expected, if phonons are relevant to high temperature superconductivity. Moreover, this Letter demonstrates that high-energy resolution inelastic x-ray scattering has developed into an invaluable tool for the study of the lattice dynamics of complex transition metal oxides, allowing measurements on small high-quality single crystals which are inaccessible to the traditional method of inelastic neutron scattering.

###### Acknowledgements.

We acknowledge L. Paolasini and G. Monaco for useful discussions and H. Casalta for precious help during preliminary tests. The authors are grateful to D. Gambetti, C. Henriquet and R. Verbeni for technical help, to J.-L. Hodeau for help in the crystal orientation and J. -P. Vassalli for crystal cutting. P.K.M. and M.G. are supported by the U.S. Department of Energy under Contracts No. DE-FG03-99ER45773 and No. DE-AC03-76SF00515, by NSF CAREER Award No. DMR-9985067, and by the A.P. Sloan Foundation.


text:
# Comment on “Fun and frustration with quarkonium in a 1+1 dimension,” by R. S. Bhalerao and B. Ram [Am. J. Phys. 69 (7), 817–818 (2001)]

## A

Let me show that Eq. (13) has an infinite number of solutions if $`m=0`$. Indeed, in this case Eq. (13) reduces to

$$f(\nu )H_\nu ^2(0)2\nu H_{\nu 1}^2(0)=0.$$
(A1)
It is clear that $`\nu `$ cannot be negative, for then $`f(\nu )`$ is strictly positive. If $`\nu =0`$ or $`\nu =2k`$, $`k`$ a positive integer, then $`2\nu H_{\nu 1}^2(0)`$ vanishes (in the latter case because $`H_{2k1}(x)`$ is an odd function of $`x`$) and $`H_\nu ^2(0)>0`$. Therefore, $`f(\nu )>0`$ if $`\nu `$ is a non-negative even integer. On the other hand, $`f(\nu )<0`$ if $`\nu `$ is a positive odd integer, for then $`H_\nu (0)=0`$ and $`2\nu H_{\nu 1}^2(0)>0`$. Because $`f(\nu )`$ is a continuous function of $`\nu `$, one is led to the conclusion that $`f(\nu )`$ vanishes at least once in each interval $`(n,n+1)`$, $`n=0,1,2,\mathrm{}`$ A more refined analysis, in which one uses the identity (cf. Eq. (6))

$$H_\nu (0)=\frac{2^\nu \mathrm{\Gamma }\left(\frac{1}{2}\right)}{\mathrm{\Gamma }\left(\frac{1\nu }{2}\right)},$$
(A2)
shows that $`f(\nu )`$ in fact vanishes only once in these intervals. It follows that the energy levels in the massless case are given by

$$E=\pm \sqrt{2\nu _ng},$$
(A3)
with $`\nu _n`$ satisfying $`n<\nu _n<n+1`$, $`n=0,1,2,\mathrm{}`$


text:
# High Density Neutron Star Equation of State from 4U 1636-53 Observations

## I introduction

Oscillations in x-ray brightness during bursts have been observed from ten low mass x-ray binaries with the Rossi X-ray Timing Explorer Strohmayer . These oscillations appear to be due to rotational modulation of one or a pair of antipodal hotspots on the surface of the accreting neutron star Strohmayer et al. (1997). With this model the amplitude of the oscillations is determined by the mass-to-radius ratio of the neutron star, the compactness $`M/R`$. A more compact neutron star results in greater deflection of the x-ray photons emanating from a hot spot, making the hot spot visible to an observer for a larger portion of the rotation period and so reducing the amplitude of the brightness oscillations. The bounds on the compactness of the neutron star are particularly limiting when two antipodal hot spots are present Miller and Lamb (1998).

The source 4U 1636-53 appears to be the best candidate for providing a stringent limit on the compactness of neutron stars from x-ray burst oscillations Nath et al. . A 290 Hz subharmonic of the stronger 580 Hz brightness oscillation has been observed in five bursts from this source Miller (1999), but the subharmonic has not been confirmed in subsequent bursts Strohmayer . The existence of the subharmonic would suggest the spin frequency of the neutron star is 290 Hz, and the presence of two antipodal hot spots is responsible for the 580 Hz modulation. A pair of high frequency quasiperiodic oscillations separated by 251 Hz has also been reported for this source and interpreted in a beat frequency model to be a consequence of a spin frequency of 290 Hz with two hot spots on the neutron star Miller et al. (1998). The x-ray oscillations during the rising phase of bursts from 4U 1636-53 have been recently modeled as due to a pair of circular antipodal hot spots, each having an angular size that grows linearly in time Nath et al. . The two hot spot model constrains the neutron star compactness to $`M/R<0.163`$ at 90% confidence (using $`G=c=1`$), requiring the radius of a 1.4 $`M_{}`$ neutron star to be greater than 12.7 km. The 99% confidence level compactness limit, $`M/R<0.183`$, implies $`R>11.3`$ km when $`M=1.4M_{}`$.

During a burst event in low mass x-ray binaries the oscillation frequency is usually observed to evolve asymptotically to a limit that probably corresponds to the spin rate of the neutron star. Measurement of the asymptotic frequency at different orbital phases of the binary could reveal a Doppler modulation that would allow for a determination of the masses of the neutron star and the companion. A mass measurement together with the compactness limit would tightly constrain the radius of the star. A Doppler modulation has not been observed for 4U 1636-53 in an analysis of 26 burst oscillations Giles et al. . The absence of the modulation likely limits the neutron star mass in this system to $`M<1.6M_{}`$ Giles et al. .

The purpose of this paper is to use the 90% confidence level compactness bound on the 4U 1636-53 neutron star to estimate the high density equation of state of neutron star matter. Ideally one would have both a measured mass and radius for a number of neutron stars over a wide range of masses up to the maximum mass value. The observed mass-radius relationship could then be used to uniquely calculate the equation of state by inverting the Oppenhemier-Volkoff equation, the relativistic stellar structure solution to Einstein’s equations Lindblom (1992). Even if the mass and radius are known for only a single neutron star, the high density equation of state can still be accurately estimated up to the central density yielding the mass of this star. At low density the equation of state can be determined by fitting experimental nucleon-nucleon scattering data and properties of light nuclei to two-body and three-body interaction potentials Wiringa et al. (1988); Akmal et al. (1998). Confidence in these low density equations of state is high near normal nuclear density (energy density = 153 MeV/fm3, baryon density = 0.16 fm-3) because nuclear physics experiments allow for verification in the low density regime. Estimates of the maximum possible mass of neutron stars based on these equations of state have regarded them as experimentally verified up to twice normal nuclear density Kalogera and Baym (1996); Olson (2001). Assuming one has values of both the mass and radius for a single neutron star, the equation of state in the density range between the experimentally verified low density regime and the high central density of this star can be accurately estimated to have a polytrope form provided the matter does not undergo a phase transition in this density range Lindblom (1992). The polytrope high density equation of state is chosen to match the energy density and pressure of the low density equation of state at the upper limit of applicability of the low density equation, and the adiabatic index of the polytrope is chosen to yield a stellar radius corresponding to the observed radius of the star when the mass of the model star equals the observed mass. The resulting polytrope is an estimate of the equation of state of neutron star matter in the high density regime.

In the case of the neutron star in 4U 1636-53 we do not have measurements of both its mass and radius, but only an upper bound on the compactness and an upper limit on the mass. Nevertheless the high density equation of state is still tightly constrained by the compactness limit provided it is assumed the neutron star mass is at least $`1.4M_{}`$, a plausible assumption for an old accreting system like 4U 1636-53. In the next section the high density equation of state is estimated from the 90% confidence level compactness bound on the neutron star in 4U 1636-53 assuming the neutron star mass is at least $`1.4M_{}`$. The main result of this paper is that the A18 + $`\delta v`$ \+ UIX\* equation of state of Ref. Akmal et al. (1998) (hereafter the APR equation of state), a recently developed equation of state representative of the class of equations of state derived from fits of experimental nuclear physics data to interaction potentials, is too soft at high density to yield a star with the observed low compactness of the neutron star in 4U 1636-53. Further if the APR equation of state is used up to twice normal nuclear energy density, and any polytrope equation of state having a sound speed less than the speed of light is used above twice normal nuclear density, the model star compactness still exceeds the compactness of the 4U 1636-53 neutron star. Using the APR equation of state up to just a value in the range 1.55-1.64 times normal nuclear energy density and a stiff polytrope with an adiabatic index $`\mathrm{\Gamma }`$ of 5.5 $`\mathrm{\Gamma }`$ 7 at high density can produce a star with the low compactness of 4U 1636-53 without the speed of sound waves in the neutron star interior exceeding the speed of light.

## II estimation of the high density equation of state

The equation of state of high density neutron star matter is estimated in this section from the compactness bound on the neutron star in 4U 1636-53. The equation of state at low density is assumed to be known for energy densities less than a transition density $`\stackrel{~}{\rho }`$. Here the APR equation of state is used for the low density portions of the stellar interior, from normal nuclear density $`\rho _{nn}`$ up to $`\stackrel{~}{\rho }`$. The BPS equation of state Baym et al. (1971) is used for densities below $`\rho _{nn}`$ to model the matter outward through the neutron star crust to the surface. Recent work on determining the maximum possible mass of neutron stars has assumed the APR equation is accurate up to twice normal nuclear energy density Kalogera and Baym (1996); Olson (2001). For densities greater than $`\stackrel{~}{\rho }`$ an equation of state of polytrope form is used in this paper:

$$\frac{p}{\stackrel{~}{p}}=\left(\frac{\rho }{\stackrel{~}{\rho }}\right)^\mathrm{\Gamma }.$$
(1)
The polytrope form can be expected to be a good approximation to the actual high density equation of state of the 4U 1636-53 neutron star provided the matter has not undergone a phase change to a softer condensed phase, such as a kaon condensate or a deconfined quark phase, at some density above $`\stackrel{~}{\rho }`$ Lindblom (1992). The compactness limit on the 4U 1636-53 neutron star requires the stellar radius to be large, so it is unlikely this star contains a condensed phase.

The algorithim for determining the high density equation of state begins by choosing the transition energy density $`\stackrel{~}{\rho }`$ where one switches between the high and low density equations of state. The pressure $`\stackrel{~}{p}`$ at the transition density is computed from the low density equation of state and the transition energy density value: $`\stackrel{~}{p}=p\left(\stackrel{~}{\rho }\right)`$. The adiabatic index $`\mathrm{\Gamma }`$ in the polytrope high density equation of state Eq. (1) is treated here as a free parameter that is varied from small values (soft equations of state) to progressively larger values (stiffer equations of state) to find the set of model stars satisfying the 4U 1636-53 neutron star compactness bound.

The interior structure of each model star is assumed to be described accurately by the Oppenheimer-Volkoff solution, the solution to Einstein’s equations for the interior of a nonrotating star in general relativity:

$$\frac{dp}{dr}=\frac{\left(\rho +p\right)\left[m(r)+4\pi r^3p\right]}{r\left[r2m(r)\right]},$$
(2)
where

$$m(r)=_0^r4\pi \rho \left(\xi \right)\xi ^2𝑑\xi .$$
(3)
Rotational corrections to the Oppenheimer-Volkoff solution are small except for very rapidly rotating stars Friedman et al. (1986) and will be ignored here. Families of model stars, one family for each $`(\stackrel{~}{\rho },\mathrm{\Gamma })`$ pair, are determined by integrating the Oppenheimer-Volkoff solution Eq. (2) for a set of central energy density values. The polytrope equation of state Eq. (1) is used for densities exceeding $`\stackrel{~}{\rho }`$. The APR equation of state is used for densities below $`\stackrel{~}{\rho }`$ and above normal nuclear density. The BPS equation of state is used below $`\rho _{nn}`$.

Figure 1 shows how the neutron star mass depends on the radius for the case where $`\stackrel{~}{\rho }=2\rho _{nn}`$. Several families of model stars are represented differing by the value of the adiabatic index of the high density polytrope. Also shown is the neutron star mass-radius curve resulting when the APR equation of state is used at all densities. Lines representing the 90% and 99% confidence level 4U 1636-53 neutron star compactness bounds delineate the allowed regions of the mass-radius parameter space. The allowed regions are below these compactness limit lines and above the 1.4$`M_{}`$ line. The compactness of all these model stars exceeds the 90% confidence level 4U 1636-53 neutron star compactness bound, $`M/R<0.163`$, even for very stiff polytropes. Models with $`\mathrm{\Gamma }>7`$ are not physically realistic because the equation of state predicts the speed of sound waves,

$$v_s=\sqrt{\frac{\mathrm{\Gamma }p}{\rho }},$$
(4)
exceeds the speed of light near the stellar center. Hence model stars with the APR equation of state used up to twice normal nuclear density and a polytrope at high density are ruled out at the 90% confidence level. High density polytropes that are at least moderately stiff ($`\mathrm{\Gamma }3`$) are consistent with the 99% confidence level bound. The model using the APR equation of state at all densities is also consistent with the 99% confidence level bound, $`M/R<0.183`$, but is ruled out at 90% confidence.

Model stars satisfying the 90% confidence level compactness bound can result if the transition energy density is chosen to be $`1.64\rho _{nn}`$ or less. For $`\stackrel{~}{\rho }>1.64\rho _{nn}`$ the speed of sound at the stellar center exceeds the speed of light for all model stars having an adiabatic index value large enough to yield a compactness below the 90% confidence level bound, and therefore are ruled out. For $`\stackrel{~}{\rho }1.64\rho _{nn}`$ the 90% confidence level bound can be met without the sound speed exceeding the speed of light. Figure 2 shows neutron star mass vs. radius curves for the case where $`\stackrel{~}{\rho }=1.64\rho _{nn}`$. The radius is near 12.7 km for a 1.4 $`M_{}`$ star for the model families having $`\mathrm{\Gamma }`$ near 7, and the sound speed is subluminal throughout the star provided $`\mathrm{\Gamma }<7.023`$. For $`\mathrm{\Gamma }=7.023`$ the sound speed just reaches the speed of light at the center of the star. The central energy density and central pressure of the 1.4 $`M_{}`$, 12.7 km star are 317 MeV/fm3 and 45.1 MeV/fm3, respectively, when $`\stackrel{~}{\rho }=1.64\rho _{nn}`$ and $`\mathrm{\Gamma }=7.023`$.

Table 1 lists particular values of the pressure as a function of the energy density for the $`\stackrel{~}{\rho }=1.64\rho _{nn}`$, $`\mathrm{\Gamma }=7.023`$ equation of state. Model stars that have $`M=1.4M_{}`$, $`R`$ = 12.7 km, and subluminal sound speeds are also obtained for $`\stackrel{~}{\rho }<1.64\rho _{nn}`$. For $`\stackrel{~}{\rho }=1.60\rho _{nn}`$, $`\mathrm{\Gamma }=6.25`$ in order to make the 1.4 $`M_{}`$ star have a radius of 12.7 km. The 1.4 $`M_{}`$, 12.7 km star has $`\mathrm{\Gamma }=5.55`$ when $`\stackrel{~}{\rho }=1.55\rho _{nn}`$. Figure 3 compares the equations of state for four cases: $`\stackrel{~}{\rho }=1.64\rho _{nn}`$ with $`\mathrm{\Gamma }=7.023`$, $`\stackrel{~}{\rho }=1.60\rho _{nn}`$ with $`\mathrm{\Gamma }=6.25`$, $`\stackrel{~}{\rho }=1.55\rho _{nn}`$ with $`\mathrm{\Gamma }=5.55`$, and the APR equation of state used at all densities. The three equations of state that use a high density polytrope differ very little from each other. Hence the equation of state of Table 1 can be regarded as an estimate of the minimally stiff equation of state of neutron star matter satisfying the 4U 1636-53 90% confidence level neutron star compactness bound that also maintains a subluminal sound speed throughout the stellar interior. A stiffer equation of state is implied if the 4U 1636-53 neutron star is actually more massive than 1.4 $`M_{}`$. Figure 3 also shows the equations of state utilizing a high density polytrope are significantly stiffer than the APR equation of state at high density.

## III discussion

A new equation of state (Table 1) valid for describing neutron star matter up to the central density of a 1.4 $`M_{}`$ star has been developed in this paper based on observations of the low mass x-ray binary 4U 1636-53. The validity of this equation of state rests on three main assumptions: the actual compactness of the 4U 1636-53 neutron star does not exceed the 90% confidence level bound, the mass of the 4U 1636-53 neutron star is near 1.4 $`M_{}`$, and the APR equation of state is valid up to near 1.64 times normal nuclear energy density. The APR equation of state used at all densities is sufficient to describe the structure of the 4U 1636-53 neutron star if its compactness is actually near the 99% confidence level bound, $`M/R=0.183`$, rather than near the 90% confidence level bound value of $`M/R=0.163`$. The APR equation of state used at all densities is also adequate to meet the 90% confidence level bound if the 4U 1636-53 neutron star is 1.27$`M_{}`$ or less. (The APR equation of state is seen from Fig. 1 to yield a radius of 11.5 km for a mass of 1.27$`M_{}`$, making the compactness of this star $`M/R=0.163`$.) If the mass of the 4U 1636-53 neutron star is larger than 1.4$`M_{}`$, the equation of state of neutron star matter must be even stiffer than the Table 1 equation of state. The APR equation of state is based on a large body of experimental evidence at low density, so it is reasonable to trust its validity up to $`\stackrel{~}{\rho }1.6\rho _{nn}`$. Figure 3 shows the equation of state derived here (Table 1) is relatively insensitive to the exact transition energy density value as long as $`\stackrel{~}{\rho }`$ is near 1.6$`\rho _{nn}`$.

If the APR equation of state is used up to $`1.55\rho _{nn}\stackrel{~}{\rho }1.64\rho _{nn}`$, a stiff polytrope, having an adiabatic index of 5.5 $`\mathrm{\Gamma }`$ 7, is needed to produce a model star with a compactness $`M/R=0.163`$ without the speed of sound waves in the neutron star interior exceeding the speed of light. The adiabatic index would be small, near 5/3 for nonrelativistic nucleons and 4/3 in the relativistic limit, if the pressure is due mainly to the Fermi kinetic energy. The adiabatic index is near 2 if the equation of state is dominated by static two-body interactions, near 3 if static three-body interactions are most important, and in excess of 3 if repulsive momentum-dependent interactions are dominant Akmal et al. (1998). The large adiabatic index of the equation of state at high density implied by the 4U 1636-53 neutron star 90% confidence level compactness bound suggests the APR equation of state underestimates at high density the strength of the repulsive momentum-dependent interactions.

###### Acknowledgements.

This work was supported by the U.S. National Aeronautics and Space Administration through Grant No. NAG5-9541.


text:
# MERLIN radio detection of an interaction zone within a binary Orion proplyd system.

## 1 Introduction

It is now well established that the compact, highly ionized gaseous knots (LV 1-6) near the Trapezium stars in M42, discovered by Laques & Vidal (1979) contain young low mass stars (YSOs) still partially cocooned in their primaeval material. See Figure 1 for their locations. Many similar objects were found in the same vicinity at radio wavelengths (Churchwell et al. 1987; Garay 1987; Garay, Moran, & Reid 1987; Felli et al. 1993a; Felli et al. 1993b). The LV knots are particularly interesting both because of their close proximity to the Sun, permitting detailed observation, and because of their extreme local environment. This environment is swept by the energetic particle wind of the nearby O6 star, $`\theta ^1`$Ori C and the LV knots are irradiated by its Lyman (and near UV) photons. Their outside surfaces are consequently highly ionized, though shielded form the direct wind from $`\theta ^1`$Ori C by a standoff bowshock (Figure 3) in the photoevaporated flow.

*HST* imagery revealed the structure of these systems, now known as ‘proplyds’ (O’Dell, Wen, & Hu 1993), with startling clarity (O’Dell & Wong 1996; Bally et al. 1998). The central stars, found originally on a 2 $`\mu `$m IRCAM ‘engineering’ image with UKIRT by Meaburn (1988) and subsequently by McCaughrean & Stauffer (1994) are seen to be embedded in dense cocoons/disks of dusty molecular gas. Photoevaporated flows from the ionized surfaces form cometary tails pointing away from $`\theta ^1`$Ori C as they meet the UV flux of this star. Supersonic flow velocities of $``$ 50 km s-1 were found (Meaburn 1988; Meaburn et al. 1993; Massey & Meaburn 1995) using the Manchester Echelle spectrometer (MES - Meaburn et al. 1984) on a variety of telescopes. The origin of these ionized flows has been modelled, initially as a two-wind interaction (Henney et al. 1996, 1997) and more recently as the interaction of the flow from the proplyd disk with the ionizing Lyman photon field of $`\theta ^1`$Ori C and the nebula (Johnstone, Hollenbach, & Bally 1998; Henney & Arthur 1998; Henney & O’Dell 1999; Störzer & Hollenbach 1999). The two-wind interaction is now thought to manifest itself much further from the proplyd disk in the form of the standoff bowshocks (one of which is shown in Figures 1 and 3) observed by Hayward, Houck and Miles (1994), McCullough et al. (1995) and Bally et al. (1998). The flow from the proplyd itself is most likely caused by photoevaporation of the disk by either ionizing (EUV) or non-ionizing (FUV) photons from $`\theta ^1`$Ori C, depending on the proplyd’s proximity to it (Johnstone, Hollenbach, & Bally 1998). The process of photoevaporation, photoionization, and acceleration of the disk gas down the surrounding density gradient is found to produce flow velocities of only 2-3 times the sound speed of $``$10 km s-1 (Dyson 1968), which does not fully account for the velocities observed in the longslit spectra. These longslit spectra were complemented by the Fabry-Perot work of O’Dell et al. (1997), which compared velocity data with *HST* images and suggested the presence of a micro-jet associated with one of the proplyds.

The presence of predominantly monopolar jets from the YSOs in many of these proplyds was first suggested by the $``$ 100 km s-1, spatially compact ($``$ 1 arcsec across) spikes or knots on the MES position-velocity (pv) arrays of Meaburn et al. (1993) and Massey & Meaburn (1995). This has now been confirmed through direct *HST* imagery by Bally, O’Dell, & McCaughrean (2000), who find more than 20 jets associated with various proplyds. This apparent ubiquity of jets and the fact that it is difficult to get high velocity flows from density gradients suggest that most collimated high velocity features are jets. *HST* imagery, as a consequence of its high angular resolution, along with high resolution near infra-red imagery (Petr 1998), has also revealed the binary nature of a small fraction of the proplyds.

The binary structure of the prominent proplyd originally designated LV 1 when unresolved in ground-based observations is of particular interest. It was first resolved into two sources in the 2 cm VLA A-array radio maps of Felli et al. (1993a) (Sources 5 and 20 in their Tables 2 and 3). The higher resolution *HST* image shows that this system is formed by two proplyds separated by 0.4″, which are catalogued as 168–326N and as 168–326S, by O’Dell & Wen (1994). In the subsequent full catalogue of stars and compact objects of O’Dell & Wong (1996), it appears as a single proplyd, 168–326, corresponding to Felli et al.’s (1993a) source 20. However, the near infra-red 0.13″ resolution binary star survey of Petr (1998) identifies 168–326 as a binary system consisting of two components that are referred to as 168–326E and 168–326W. Bally et al. (1998), who used the same *HST* observations as used here, label the two components as 168–326 and 168–326a (their figure 2) and mention the presence of a bright arc lying between the sources. Here, the proplyds are referred to as 168–326SE and 168–326NW, which is the clearest form of designation. In the present paper the discovery is reported of the interaction zone between the flows from 168–326SE and 326NW in both the radio (at 5 GHz with MERLIN) and optical (H$`\alpha `$ and \[O iii\] 5007 Å with the *HST*) domains. The collision of two photoevaporated flows is explored to explain the nature of this interaction. Further kinematical and morphological evidence for a jet in LV 1 is also presented.

## 2 Observations

### 2.1 Imagery

#### 2.1.1 Radio

MERLIN observations of the Orion nebula were made between 1998 December 21 and 1999 February 15 at 5 GHz (6 cm) with a bandwidth of 15 MHz. The observations, taken over 8 days, alternated every 10 minutes between a phase calibrator, 0539-057 and the target itself. 3C286 was used to get a flux for a point source at the time of the observations. The phase calibrator was used in order to remove atmospheric density variations from the target source data and also to register the astrometric alignment of the images to an accuracy approximately 0.010″. Initial data editing and amplitude calibration were done using the MERLIN-specific d-programs. Subsequent phase calibration was then carried out in the National Radio Astronomy Observatories’ (NRAO) AIPS package. The resulting data were CLEANed using the AIPS task IMAGR and finally restored with a circular beam of full width at half-maximum, FWHM = 0.080″. More details of the reduction of the radio data can be found in a parallel paper on LV 2 (Henney et al. 2001).

#### 2.1.2 Optical

The optical observations presented here were obtained with the *HST*’s Planetary Camera (PC) under J. Bally’s observing program, GO 5469. The images used were obtained at the target named NGC1976-LV3 through the F656N (H$`\alpha `$) and F502N (\[O iii\] 5007 Å) filters. For each filter there were three images at this position, each lasting 60 s, for the H$`\alpha `$ filter and 100 s for the \[O iii\] 5007 Å filter. In both cases, the three images were combined using the NOAO IRAF package’s CRREJ routine, which removes cosmic rays during the combining process. The images were then rotated by an amount such that north was at the top of the image and east was at the left, using the image header keyword ORIENTAT, which gives the position angle of the y axis of the PC detector projected on the sky, in degrees east of north.

Images in the F658N and F547M filters were also used to make an exact flux calibration using the method of O’Dell and Doi (1999), which corrects for underlying continuum and contamination by adjacent lines.

#### 2.1.3 Alignment of Radio and Optical Images

The process of phase referencing MERLIN images results in an absolute astrometric calibration that is accurate to approximately 0.010″ and tied to the radio reference frame (ICRF). However the absolute co-ordinate system that is built in to *HST* images is known to be inferior to this by an order of magnitude. Therefore, alignment of the optical and radio images was done by shifting the co-ordinate system of the optical images and taking the radio co-ordinate system to be correct. Ideally, the optical images could be re-calibrated empirically using stars in the field with *HIPPARCOS* positions, of which there are three. Unfortunately, these three stars, $`\theta ^1`$Ori A (HIP 26220), $`\theta ^1`$Ori C (HIP 26221), and $`\theta ^1`$Ori D (HIP 26224) are badly saturated in the *HST* PC images, meaning that their positions on the PC CCD chip cannot be reliably determined.

In the vicinity of $`\theta ^1`$Ori C, MERLIN clearly detected all six LV knots at 5 GHz. These six proplyds are also obvious on the *HST* images (figure 2 of Bally et al. 1998). A comparison of the centroided optical and radio positions as given by their respective co-ordinate systems leads to a consistent difference between the two. As a result, this difference of a few hundred milliarcseconds was applied to the co-ordinate system of the *HST* images, in order to correct and align them with the MERLIN image. It is not unreasonable to expect that there is a real difference in the position of the optical and radio emission from these objects, which puts this method into question. However, this difference is likely to be a lot less than the gross error of a few hundred milliarcseconds in the *HST*’s astrometric calibration. Also, this real separation in optical and radio position would be expected to be along the line joining the proplyd to the ionizing source ($`\theta ^1`$Ori C). Therefore, the fact that the five proplyds used for the alignment surround $`\theta ^1`$Ori C and are oriented at various position angles means that the real difference is somewhat preserved. In effect, the difference in radio and optical positions has been minimised as opposed to removed. It should be noted that LV 6 was not used because its radio position was hard to determine, due to its apparently complicated structure. $`\theta ^1`$Ori A is a strong radio source, but could not be used in the alignment as its radio position was found to be significantly displaced from its optical position. The final astrometric calibration was found to be consistent with that of Henney et al. (2001).

### 2.2 Echelle Spectroscopy

The Trapezium region of the Orion nebula, shown in Figure 1, has been explored spectroscopically with longslit spectra obtained with the MES (Meaburn 1988; Massey & Meaburn 1993; Meaburn et al. 1993; Massey & Meaburn 1995).

The pv array of \[O iii\] 5007 Å profiles shown in Figure 2 was taken at the slit position referred to as S1 in Massey & Meaburn (1995) with the MES on the 2.5m Isaac Newton Telescope.

The slit position is shown in Figure 1 and the full observational setup for the data can be found in Massey & Meaburn (1995). Stellar continuum from $`\theta ^1`$Ori C has now been subtracted from this spectrum, leaving some residual noise, but also allowing the emission from LV 1 to be more easily appreciated.

The emission from LV 1 is identified in Figure 2. It is most intense out to velocities of $`\pm 50\mathrm{km}\mathrm{s}^1`$ with respect to the nebula and extends more faintly to $``$ 100 km s-1 with respect to the systemic heliocentric radial velocity, $`V_{\mathrm{sys}}`$ = 25 km s-1. Such a high speed ‘spike’ in the pv array is indicative of a jet and although one was not originally identified in the *HST* imagery of Bally, O’Dell, & McCaughrean (2000), García-Arredondo, Henney & Arthur (2001) present evidence of a micro-jet from 168-326SE at a position angle of +85. This jet is shown in Figure 3 to be predominantly monopolar, although there is a hint of an approaching component in the PV array in Figure 2.

## 3 Radio/Optical Comparisons

The identification chart for the Trapezium region in Figure 1 shows $`\theta ^1`$Ori A-F, the six LV objects, and the bowshock that marks the interface between the low speed proplyd disk flow of LV 1 and the high speed stellar wind of $`\theta ^1`$Ori C. This interaction shields the LV 1 binary system from a direct interaction with $`\theta ^1`$Ori C’s wind and happens on a much larger scale than the interaction between the two proplyds which constitute LV 1 (168–326SE and 168–326NW).

Three distinct sources can be seen in the MERLIN 5 GHz radio contours shown in Figures 4 and 5, which compare the 5 GHz radio emission to the H$`\alpha `$ and \[O iii\] 5007 Å emission respectively.

The result of the alignment techniques described above is that the south-east radio source corresponds to 168–326SE and the north-west radio source corresponds to 168–326NW. This leads to the strong possibility that the radio source located between the binary components is the result of an interaction between them.

This interaction zone is not as prominent in the H$`\alpha `$ image. Furthermore, the peak radio to peak H$`\alpha `$ flux ratios for the two proplyds are the same, whereas this ratio for the interaction source is just over 2.56 times that for the proplyds. This could indicate that the interaction radio emission has a non-thermal origin (see Section 5).

However, the interaction zone can be seen in the light of the \[O iii\] 5007 Å line, as shown in Figure 4. This could mean that the region has a temperature high enough (i.e. $`>10^4`$ K) to radiate predominantly in \[Oiii\], as opposed to H$`\alpha `$. A different density in this region compared with the proplyds, combined with this higher temperature could therefore generate sufficiently enhanced thermal radio emission (see $`\mathrm{\S }`$ 5) to give the enhanced radio/H$`\alpha `$ ratio that is observed without invoking non-thermal radio processes.

The position of the brightest part of the interaction source in the \[O iii\] 5007 Å image is of particular interest. In this image, the tail of 168–326SE can be seen pointing almost exactly away from $`\theta ^1`$Ori C. This is a common feature of many of the proplyds (Bally et al. 1998) and leads to the conclusion that the diffuse photon flux of the nebula is causing a flow off the ‘back’ of the proplyds. The brightest point of the interaction source is on the opposite side of 168–326NW to $`\theta ^1`$Ori C and not directly on the axis joining the two binary components. A simplistic explanation for this is that the interaction is happening between material flowing off the ‘back’ of 168–326NW into the flow from the ‘front’ of 168–326SE.

## 4 Modelling the Interaction

The morphology of the binary proplyd system, LV 1, shown in Figure 5 can be summarised as follows. There is a bright interaction peak of emission at 5 GHz between the two proplyds, which is nearer the smaller proplyd, 168–326NW. This peak is the brightest part of an arc of emission in H$`\alpha `$ and \[O iii\] 5007 Å that is concave with respect to the smaller proplyd. The peak does not lie on the line connecting the two proplyd ‘heads’, but instead lies behind 168–326NW as seen from $`\theta ^1`$Ori C. This morphology seems to represent the first evidence for a hydrodynamic interaction between two proplyds. A schematic model of such an interaction is shown in Figure 6, in which the transonic photoevaporated flows of a few 10’s of km s-1 (as seen in Figure 2) from each proplyd’s disk, driven by the UV radiation from $`\theta ^1`$Ori C, are colliding to form a dense interproplyd shell bounded by two weak shocks.

The shell formed by the interaction of the disk flows from two proplyds with the stellar wind of $`\theta ^1`$Ori C is also shown. Canto, Raga, & Wilkin (1996, hereafter CRW) have studied the interaction between two hypersonic, non-accelerating, spherical winds forming an axisymmetric bowshock structure under the approximation that the interaction region forms a thin, momentum-conserving shell. The situation with the proplyds is more complicated, mainly because the individual flows are non-isotropic, which will result in a totally asymmetric bowshock. However, the following principal results of the model appear to have observational counterparts in the LV 1 system:

* There will exist a point, termed the *stagnation point*, where the shell surface is perpendicular to the velocity vectors of both flows. This forms the ‘nose’ of the bowshock and its position is determined by the condition of ram pressure balance; it will be closer to the less powerful proplyd in terms of wind momentum flux. The normal to the shell surface at the stagnation point can be considered the bowshock axis, even though it is not an axis of symmetry.
* The shell will curve back towards the less powerful flow. At each azimuth around the bowshock axis, the shell will tend asymptotically to a constant opening angle from the axis, but the anisotropy of the flows will lead to this opening angle varying with azimuth.

The stagnation point is the area of highest density in the interaction shell and will therefore show up as a peak in the intensity of any bowshock that lies between the proplyds and is therefore clearly identifiable in Figure 5. The anisotropy of the flows means that the stagnation point need not necessarily lie on a line joining the proplyd heads. The observed radio emission at the stagnation point takes the approximate form of a point source, with some evidence of extension at $`PA45^{}`$ in the lowest (3$`\sigma `$) contour level (Figures 4 & 5). If the radio emission does continue to extend in the same way as the optical, it could be below the sensitivity of these observations.

The details of deciphering the geometry of the interacting binary system and then applying the model of CRW is left to a subsequent paper (Henney 2001). However, a preliminary result of considering the 3-dimensional geometry of LV 1 leads to two possible geometrical configurations for the system, as shown in Figure 7.

Both configurations imply that the proplyds lie approximately between the observer and $`\theta ^1`$Ori C, i.e. the observer is ‘behind’ the proplyds with respect to $`\theta ^1`$Ori C. In configuration (a) the proplyds are roughly side-by-side as seen from $`\theta ^1`$Ori C, with the stagnation point forming between them. Whereas in configuration (b), the smaller proplyd is closer to $`\theta ^1`$Ori C, resulting in the tail flow from this smaller proplyd interacting with the head flow from the larger proplyd, with the stagnation point forming close to the base of the smaller proplyd’s tail, as suggested in § 3. In both cases the stagnation point will be seen to lie between the proplyds, but in configuration (a) it will be superimposed on the tail of 168–326NW. This superposition can be seen in the *HST* image in Figure 5 and therefore favours the side-by-side configuration. However, this configuration places the observer within the ‘bowshock cone’, implying that an arc would not be observed. The observed arc in the *HST* images favours configuration (b). If the observer does lie outside the bowshock cone, then the shell would be seen tangentially at some point and the fact that it is clearly resolved with the *HST* would put a lower limit on its thickness.

In order to conclusively choose between these two configurations, the morphological appearance of the shell has been modelled for both cases. For this, a model for the shape of the shell was needed, and has been provided by the CRW model. The shape of a CRW shell is characterised by its asymptotic opening angle, $`\theta _{\mathrm{}}`$, which is uniquely determined by the momentum flux ratio, $`\beta `$, of the two interacting flows. In Figure 8, this shape for $`\beta =5`$, $`\theta _{\mathrm{}}=61.9^{}`$ has been used to calculate the appearance of a shell at various viewing angles, $`\alpha `$, the angle between the line of sight and the bowshock axis for a cylindrically symmetric bowshock.

It can be seen that the morphology of the modelled shell when viewed from within the bowshock cone with $`\alpha =61`$ is close to that of the observed shell, but this modelled shell has no apparent curvature. In fact, the observed curvature is similar to a modelled shell with $`\alpha =63^{}`$, with the observer outside the bowshock cone. This discrepancy could be explained by the fact that the bowshock need not be cylindrically symmetric. Each proplyd has a distinctive cometary shape extending away from $`\theta ^1`$Ori C along an axis referred to as the proplyd axis (see Figure 7). If there is a contribution to a proplyd’s flow from the entire length of the proplyd axis, the bowshock between two proplyds is likely to be flattened in the direction perpendicular to the proplyd axis. The appearance of an asymmetric bowshock can be modelled by arbitrarily decreasing $`\beta `$ and therefore increasing $`\theta _{\mathrm{}}`$ in the direction parallel to the proplyd axis. This will result in a bowshock whose cross-section resembles an ellipse with a major axis parallel to the proplyd axis. Figure 9 is a set of modelled images for such an asymmetric bowshock, with the $`\theta _{\mathrm{}}`$ increased to 75 in the direction parallel to the proplyd axis.

With $`\alpha =0`$ the intrinsic flatness of the bowshock can be seen, giving rise to elliptical isophotes with a major axis in the horizontal direction. As $`\alpha `$ increases this horizontal elongation decreases and at $`\alpha 45^{}`$, the intrinsic elongation has been cancelled out, giving rise to approximately circular isophotes. As $`\alpha `$ continues to increase the apparent elongation is in the vertical direction and some apparent curvature is visible along the bowshock, even with the observer still inside the bowshock cone. With $`\alpha =70^{}`$, the morphology of the modelled image is a good match for the observed morphology, in terms of both smoothness and curvature (see Figures 4 & 5). In these modelled images (Figures 8 & 9), the thickness of the shell has been assumed to be unresolved, as can be seen in those images where the shell is viewed tangentially. Henney (2001) predicts that the shell’s thickness is likely to be unresolved with the *HST*. However, the presence of a magnetic field in the shell will act as an additional pressure, which will thicken the interproplyd shell. If this magnetic pressure dominates the thermal pressure in holding the shell up against the ram pressure of the two incident flows, the density in the region will be lower than that predicted by Henney (2001). A lower density will lower the cooling rate and cause the cooling length to increase. Therefore a significant magnetic field increases both the overall size of the shell and the cooling length within it. If this effect is large enough, the shell could be viewed from outside the bowshock cone, and therefore tangentially, and still be resolved. Dust (other than foreground extinction from the neutral veil) will be unimportant because the column density of the shocked shell is too low.

## 5 Nature of the Radio Emission

An interaction zone between the two components of the binary proplyd, LV 1, has been detected in both the radio and optical domains. This appears as an extended arc, containing a compact (5 GHz) knot. The morphology of the interacting proplyd binary is similar to that of other interacting binaries such as WR 147 (see e.g. Williams et al. 1997; Niemela et al. 1998). In such cases the interaction radio source is non-thermal, synchrotron emission caused by shock-accelerated, relativistic electrons. High energy electrons can be produced in this scenario, because at least one of the winds has a speed of order $`1000\mathrm{km}\mathrm{s}^1`$. In the case of the proplyds, the speed of the evaporated flows is thought to be a few times the sound speed of $`10\mathrm{km}\mathrm{s}^1`$ and the bulk of the emission from LV 1 appears to be from gas flowing at velocities of less than $`50\mathrm{km}\mathrm{s}^1`$ (Figure 2). Two interacting flows each with this velocity, would lead to two shocked layers at a post-shock temperature of $`3\times 10^4`$ K (Dyson & Williams 1997). Although this is high enough to explain the enhanced \[Oiii\] emission from the interaction region, it is not sufficient to account for the existence of relativistic electrons. Therefore, synchrotron emission is not expected from these low speed collisions.

The 5 GHz radio/H$`\alpha `$ intensity ratio is approximately 2.6 times higher for the interaction region than for the proplyds. For identical physical conditions in the proplyds and the interaction region, this could be indicative of non-thermal radio emission from the interaction zone. However, it has already been established that the post-shock material in this zone is likely to have a higher temperature than the proplyds. Post-shock material is also unlikely to have the same density as the unshocked material flowing from the proplyd. The canonical density jump for a strong shock is a factor of 4 (Dyson & Williams 1997), but for weaker shocks this factor is lower. Once the post-shock material has cooled, the final compression ratio as compared with preshocked material will be higher than this. It should also be noted that the material is flowing from the proplyd down a density, $`N`$, gradient that, for an accelerating spherical flow, would scale with radial distance from the centre of the proplyd disk, $`r`$, as follows,

$$N(r)=N(r_{\mathrm{if}})\left(\frac{r_{\mathrm{if}}}{r}\right)^2\left(\frac{V_{\mathrm{if}}}{V(r)}\right)$$
(1)
where $`r_{\mathrm{if}}`$ is the radius of the proplyd ionization front (IF), $`N(r_{\mathrm{if}})`$ is the density at this radius, $`V_{\mathrm{if}}`$ is the gas velocity at the IF ($``$ sound speed) and V(r) is the gas velocity. The interaction between 168–326SE and 168–326NW is happening at approximately 2-3 IF radii. Therefore, in the case of a spherical flow from the proplyds, a density drop of a factor of $`>4`$ would be expected in the pre-shocked proplyd flow as compared with material flowing from the IF. As a result, even after compression due to the shock, the density in the post-shock material in the interproplyd shell is likely to be lower than the density at the proplyd’s IF. This difference in both temperature and density between the IF and the interproplyd shell could provide an explanation for the observed difference in 5 GHz radio/H$`\alpha `$ intensity ratio. However, the contribution from the high temperature gas, which will rapidly cool through \[Oiii\] emission, depends on its cooling length as compared with the total length of the interproplyd shell. If this contribution from the post-shock cooling layer is insignificant, the majority of the material in the shell will be at the equilibrium nebular temperature of $`10^4`$ K. The density of this cooled material will be higher than that in the cooling layer, but still different to that at the proplyd IF.

Using the tabulated values of $`j_{\mathrm{H}\alpha }`$, the emission coefficient of the H$`\alpha `$ line provided by Storey & Hummer (1995) at various temperatures, $`T`$, the H$`\alpha `$ flux density from an object with uniform density, $`N`$, at a distance, $`D`$, for case B recombination is

$$F(\mathrm{H}\alpha )=1.850\times 10^{21}T^{0.932}N^2R^3D^2(\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1),$$
(2)
where $`R`$ is the radius of the object. Using typical values for the proplyds, equation (2) predicts flux densities of order $`10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, which is consistent with the total fluxes calculated by O’Dell (1998).

The radio flux, $`S_\nu `$, due to thermal brehmstrahlung from a source observed at a frequency, $`\nu `$, is

$$S_\nu =\frac{2k\nu ^2\mathrm{\Omega }}{c^2}T\left(1e^{\tau _\nu }\right),$$
(3)
where $`\mathrm{\Omega }`$ is the solid angle of the source and $`\tau _\nu `$ is the optical depth. This is reasonably approximated by the following relationship

$$\tau _\nu =8.235\times 10^2\left(\frac{T}{\mathrm{K}}\right)^{1.35}\left(\frac{\nu }{\mathrm{GHz}}\right)^{2.1}\left(\frac{EM}{\mathrm{pc}\mathrm{cm}^6}\right)$$
(4)
where $`EMN^2R`$ is the emission measure (Mezger & Henderson 1967). For typical values for the proplyds at 5 GHz, this is close to 1, so the proplyds can neither be treated as optically thin nor thick at this frequency. A full treatment allowing for an arbitrary optical depth is required.

Therefore, the peak H$`\alpha `$ flux ratio for two objects at temperatures, $`T_1`$ and $`T_2`$, densities, $`N_1`$ and $`N_2`$ and sizes $`R_1`$ and $`R_2`$ is given by

$$\frac{F(\mathrm{H}\alpha )_1}{F(\mathrm{H}\alpha )_2}=\left(\frac{T_1}{T_2}\right)^{0.932}\left(\frac{N_1}{N_2}\right)^2\left(\frac{R_1}{R_2}\right)$$
(5)
and the 5 GHz radio flux ratio for the two objects with optical depths, $`\tau _{5\mathrm{G}\mathrm{H}\mathrm{z}(1)}`$ and $`\tau _{5\mathrm{G}\mathrm{H}\mathrm{z}(2)}`$ is given by

$$\frac{S_{5\mathrm{G}\mathrm{H}\mathrm{z}(1)}}{S_{5\mathrm{G}\mathrm{H}\mathrm{z}(2)}}=\frac{T_1\left(1e^{\tau _{5\mathrm{G}\mathrm{H}\mathrm{z}(1)}}\right)}{T_2\left(1e^{\tau _{5\mathrm{G}\mathrm{H}\mathrm{z}(2)}}\right)}.$$
(6)
These ratios are known for 168–326SE compared to the interaction zone from the *HST* and MERLIN 5 GHz data, to be 3.4 in H$`\alpha `$ and 1.3 in the radio. Equations (5) and (6) can be solved simultaneously for $`T_2`$ and $`N_2`$, given values for $`T_1`$ and $`N_1`$, i.e. given a temperature and density for 168–326SE a temperature and density consistent with the flux ratios can be found for the interaction zone. A scale-size is needed in order to estimate the volume of material being observed. For both the proplyd and the interaction zone, this was chosen to be $`10^4`$ pc ($`20`$ AU), which is approximately the radius of the IF for a proplyd. With $`T_1=10^4`$ K and $`N_1=10^6\mathrm{cm}^3`$ for the proplyd, the equations are solved with $`T_2=25470`$ K and $`N_2=0.839\times 10^6\mathrm{cm}^3`$ for the interaction zone. These values are consistent with both the temperature predicted by the shock speed and the expectation that the density should be lower in the interaction zone. However, the values are sensitive to the chosen scale-size, initial temperature ($`T_1`$) and density ($`N_1`$). The H$`\alpha `$ fluxes could also be significantly affected by extinction. Nevertheless, this treatment shows that by changing the temperature and density in the post-shock region, the required flux ratios can be arrived at without the need for a non-thermal source of radio emission, as long as there is sufficient material at the high post-shock temperature and lower density. If there is insufficient post-shock material, as suggested by Henney (2001), the difference in density between the proplyd and the cooled material could still explain the difference in 5 GHz radio/H$`\alpha `$ intensity ratio between the two regions. The intensity ratio is also affected by the physical size of the emitting region, so the uncertainty in the shell thickness and cooling length makes it difficult to determine whether the radio emission has a thermal origin.

The spectral index of the interaction compact radio source is now needed to distinguish conclusively between a thermal or non-thermal origin. For this, observations are needed at 1.5 GHz (20 cm) with a similar angular resolution to that of the 5 GHz MERLIN observations presented here.

## 6 Conclusions

The Orion proplyd originally designated LV 1 is shown at both radio (with MERLIN) and optical wavelengths (with the HST) to be a proplyd binary system with an interaction zone separating the two components 168-326NW and 168-326SE.

Comparison of detailed hydrodynamical modelling of the interaction zone with the observations has identified the geometry of the system convincingly.

A strong 6 cm radio source coincides with the stagnation point in this interaction zone. The present observations cannot distinguish between a thermal or non-thermal origin for this radio emission. Observations at 20 cm but with a similar angular resolution (0.1 ″) are awaited to clarify this point through measurement of the radio spectral index.

The existence of a high-speed micro-jet from the binary component 168-326SE is also suggested by spectral observations and HST imagery.

Acknowledgements

MFG would like to acknowledge a studentship received from the Particle Physics and Astronomy Research Council (PPARC). MERLIN is a national facility operated by the University of Manchester on behalf of PPARC in the U.K. The work presented in this paper is based on observations made with MERLIN and the NASA/ESA Hubble Space Telescope, the latter obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.