text:
# Models of the Solar Wind Interaction with Local Interstellar Cloud

## 1 INTRODUCTION

The Local Interstellar Cloud (LIC) is a cloud of partly ionized plasma surrounding the Solar System. The plasma component of LIC interacts with the solar wind plasma and forms the heliospheric interface (Figure 1). The heliospheric interface is a complex structure, where the solar wind and interstellar plasma, interplanetary and interstellar magnetic fields, interstellar atoms of hydrogen, galactic and anomalous cosmic rays (GCRs and ACRs) and pickup ions play prominent roles.

Although a space mission into the Local Interstellar Cloud is becoming now more realisable, there are no yet direct observations inside the heliospheric interface. Therefore, at the present time the heliospheric interface structure and local interstellar parameters can be derived only from remote experiments and measurements. Currently, backscattered solar Ly-$`\alpha `$ radiation, pickup ions, anomalous cosmic rays, and Voyager measurements of distant solar wind are the major sources of information on the heliospheric interface structure and position of the termination shock . kHz emission detected by Voyager can put some constraints. Recently, it was shown that study of Ly-$`\alpha `$ absorptions toward nearby stars can serve as remote diagnostics of the heliospheric interface and, in particular, the hydrogen wall around the heliopause (e.g.,  \- ). In the foreseeable future, remote diagnostics will be also possible with images of heliospheric energetic neutrals (ENAs) . To reconstruct the structure of the interface and physical processes inside the interface on the basis of remote observations, a theoretical model should be employed.

Theoretical studies of the heliospheric interface were performed over more than four decades after pioneering papers by Parker  and Baranov et al. . However, a complete theoretical model of the heliospheric interface has not been constructed yet. The difficulty in doing this is connected with the multi-component nature of both the LIC and the solar wind. The LIC consists of at least five components: plasma (electrons and protons), hydrogen atoms, interstellar magnetic field, galactic cosmic rays, and interstellar dust. The heliospheric plasma consists of original solar particles (protons, electrons, alpha particles, etc.), pickup ions and the anomalous cosmic ray component. The pickup ion component is a result of ionization of those interstellar H atoms that penetrate into the heliosphere through the heliospheric interface. A part of the pickup ions is accelerated to high energies of ACRs. ACRs may also modify the plasma flow upstream of the termination shock and in the heliosheath. Spectra of ACRs can serve as remote diagnostics of the termination shock. For a recent review on ACRs see .

To construct a theoretical model of the heliospheric interface, one needs to choose a specific approach for each interstellar and solar wind component. Interstellar and solar wind protons and electrons can probably be described as fluids. At the same time interstellar H atom flow requires kinetic description. For pickup ion and cosmic ray components, the kinetic approach is also required. However, for interpretations that are not directly connected to pickup ions and ACRs, a cruder model can be used.

This paper focuses on the models of the global heliospheric interface structure. Under global models I understand those models that study the whole interaction region, including the termination shock, the heliopause and possible bow shock. In this sense, this paper should not be considered as a complete review of progress in the field. Many different approaches were used to look into different aspects of the solar wind interaction with LIC connecting with pickup ion transport and acceleration, with the termination shock structure under influence of ACRs and pickup ions. For more complete overview see recent reviews ,.

The structure of the paper is the following: The next section briefly describes our current knowledge of the local interstellar and solar wind parameters. Section 3 discusses theoretical approaches to be used for the interstellar and solar wind components. Section 4 gives an overview of heliospheric interface models. Section 5 describes basic results of the Baranov-Malama model of the heliospheric interface and its future developments. In section 6 we demonstrate possible analyses of space experiments on the basis of a theoretical model of the heliospheric interface. Section 7 underlines current problems in the modeling of the global heliosphere and discusses future perspectives.

## 2 BRIEF SUMMARY OF OBSERVATIONAL KNOWLEDGE

Choice of an adequate theoretical model of the heliospheric interface depends on boundary conditions, i.e. on undisturbed solar wind and interstellar parameters.

### 2.1 Solar wind observations

At the Earth’s orbit the flux of interstellar atoms is quite small, and the solar wind can be considered undisturbed. Measurements of pickup ions and ACRs also show that these components do not have dynamical influences on the original solar wind particles at the Earth’s orbit. Therefore, solar wind parameters at the Earth’s orbit can be taken as inner boundary conditions.

It has been shown by many authors that pickup and ACR components dynamically influence the solar wind at large heliocentric distances. Observable evidence of such influence is, for example, deceleration of the solar wind detected by Voyager . Table 1 presents estimates of dynamic importance of the heliospheric plasma components at small and large heliocentric distances. The table shows that pickup ion thermal pressure can be up to 30-50 % of the dynamic pressure of solar wind.

### 2.2 Interstellar parameters

Local interstellar temperature and velocity can be inferred from direct measurements of interstellar atoms of helium by Ulysses/GAS instrument . Atoms of interstellar helium penetrate the heliospheric interface undisturbed, because of the small strength of their coupling with interstellar and solar wind protons. Indeed, due to small cross sections of elastic collisions and charge exchange with protons, the mean free path of these atoms is larger than the heliospheric interface. Independently, the velocity and temperature in the Local Interstellar Cloud can be deduced from analysis of absorption features in the stellar spectra . However, this method provides mean values along the line of sight in the LIC. A comparison of local interstellar temperatures and velocities derived from stellar absorption with those derived from direct measurements of interstellar helium shows quite good agreement (see Table 2).

Other local parameters of the interstellar medium, such as interstellar H atom and electron number densities, and strength and direction of the interstellar magnetic field, are not well known. In the models they can be considered as free parameters. However, measurements of interstellar H atoms and their derivatives as pickup ions and ACRs provide important constraints on local interstellar densities and total pressure. The neutral H density in the inner heliosphere depends on filtration the neutral H atoms in the heliospheric interface due to charge exchange. Since interstellar He is not perturbed in the interface, local interstellar number density of H atoms can be estimated from the neutral hydrogen to the neutral helium ratio in the LIC, $`R(HI/HeI)_{LIC}`$: $`n_{LIC}(HI)=R(HI/HeI)_{LIC}n_{LIC}(HeI)`$. The neutral He number density in the heliosphere has been recently determined to be very likely around $`0.0130.018`$ cm-3 ( \- ). Interstellar ratio HI/HeI is likely in the range of 10-14. Therefore, expected interstellar H atom number densities are in the range of $`0.130.25`$ cm-3. It was shown by modeling ,  that the filtration factor, which is the ratio of neutral H density inside and outside the heliosphere, is a function of interstellar plasma number density. Therefore, the number density of interstellar protons (electrons) can be estimated from this filtration factor . Independently, the electron number density in the LIC can be estimated from abundances ratios of ions of different ionization states .

Note that there are other methods to estimate interstellar H atom density inside the heliosphere, based on their influence on the distant solar wind  or from ACR spectra . Recent estimates of the location of the heliospheric termination shock using transient decreases of cosmic rays observed by Voyager 1 and 2 also provide constraints on the local interstellar parameters . However, simultaneous analysis of different types of observational constraints has not been done yet. Theoretical models should be employed to make such analysis. Table 2 presents a summary of our knowledge of local interstellar parameters. Using these parameters, we estimate local pressures of different interstellar components (Table 3). All pressures have the same order of magnitude. This means that theoretical models should not neglect any of these interstellar components. Dynamical pressure of interstellar H atoms is larger than all other pressures. A part of H atoms, ACRs and GCRs penetrate into the heliosphere, which makes their real dynamical influence on the heliospheric plasma interface difficult to estimate.

## 3 OVERVIEW OF THEORETICAL APPROACHES

In this section we consider theoretical approaches for components involved in the dynamical processes in the heliospheric interface.

Generally, any gas can be described on a kinetic or a hydrodynamic level. In the kinetic approach, macroscopic parameters of a gas of $`s`$-particles (or, briefly, $`s`$-gas) can be expressed through integrals of velocity distribution function $`f_s(\stackrel{}{r},\stackrel{}{w},t)`$: $`n_s=f_s𝑑\stackrel{}{w}`$, $`\stackrel{}{V}_s=(\stackrel{}{w}f_s𝑑\stackrel{}{w})/n_s`$, $`P_{s,ij}=m_s(w_iV_{s,i})(w_jV_{s,j})f_s𝑑\stackrel{}{w}`$, $`\stackrel{}{q}_s=0.5m_s(\stackrel{}{w}\stackrel{}{V_s})^2(\stackrel{}{w}\stackrel{}{V}_s)f_s𝑑\stackrel{}{w}`$, where $`n_s`$ is the number density of $`s`$-gas, $`\stackrel{}{V}_s`$ is the bulk velocity of $`s`$-gas, $`P_{s,ij}`$ are components of the stress tensor $`\widehat{P}_s`$, $`\stackrel{}{q}_s`$ is the thermal flux vector, $`m_s`$ is the mass of individual $`s`$-particle. In the hydrodynamic approach, some assumptions should be made to specify the stress tensor $`\widehat{P}_s`$, and the thermal flux vector, $`\stackrel{}{q}_s`$ to make hydrodynamic system closed. For example, these values can be calculated by the Chapman-Enskog method, assuming $`Kn=l/L<<1`$, where $`l`$ and $`L`$ are the mean free path of the particles and characteristic size of the problem, respectively. The zero approximation of the Chapman-Enskog method gives local Maxwellian distribution, and the gas can be considered as an ideal gas, where the stress tensor reduces to scalar pressure $`P`$ and $`\stackrel{}{q}=0`$.

### 3.1 H atoms

Interstellar atoms of hydrogen form the most abundant component in the circumsolar local interstellar medium (see, Table 2). These atoms penetrate deep into the heliosphere and interact with interstellar and solar wind plasma protons. The cross sections of elastic H-H, H-p collisions are negligible as compared with the charge exchange cross section . Charge exchange with solar wind/interstellar protons determines the properties of the H atom gas in the interface. Atoms, newly created by charge exchange, have the local properties of protons. Since plasma properties are different in the four regions of the heliospheric interface shown in Figure 1, the H atoms can be separated into four populations, each having significantly different properties. The strength of H atom-proton coupling can be estimated through the calculation of mean free path of H atoms in plasma. Generally, the mean free path (with respect to the momentum transfer) of s-particle in t-gas can be calculated by the formula: $`l=m_sw_s^2/(\delta M_{st}/\delta t)`$. Here, $`w_s`$ is the individual velocity of s-particle, and $`\delta M_{st}/\delta t`$ is individual s-particle momentum transfer rate in t-gas.

Table 4 shows the mean free paths of H atoms with respect to charge exchange with protons. The mean free paths are calculated for typical atoms of different populations at different regions of the interface in the upwind direction. For every population of H atoms, there is at least one region in the interface where the Knudsen number $`Kn0.51.0`$. Therefore, the kinetic Boltzmann approach must be used to describe interstellar atoms in the heliospheric interface.

The velocity distribution of H atoms $`f_\mathrm{H}(\stackrel{}{r},\stackrel{}{w}_\mathrm{H},t)`$ may be calculated from the linear kinetic equation introduced in :

$`{\displaystyle \frac{f_\mathrm{H}}{t}}+\stackrel{}{w}_\mathrm{H}{\displaystyle \frac{f_\mathrm{H}}{\stackrel{}{r}}}+{\displaystyle \frac{\stackrel{}{F}}{m_\mathrm{H}}}{\displaystyle \frac{f_\mathrm{H}}{\stackrel{}{w}_\mathrm{H}}}=f_\mathrm{H}{\displaystyle |\stackrel{}{w}_\mathrm{H}\stackrel{}{w}_p|\sigma _{ex}^{\mathrm{HP}}f_p(\stackrel{}{r},\stackrel{}{w}_p)𝑑\stackrel{}{w}_p}`$ (1)
$`+f_p(\stackrel{}{r},\stackrel{}{w}_\mathrm{H}){\displaystyle |\stackrel{}{w}_\mathrm{H}^{}\stackrel{}{w}_\mathrm{H}|\sigma _{ex}^{\mathrm{HP}}f_\mathrm{H}(\stackrel{}{r},\stackrel{}{w}_\mathrm{H}^{})𝑑\stackrel{}{w}_\mathrm{H}^{}}(\nu _{ph}+\nu _{\mathrm{impact}})f_\mathrm{H}(\stackrel{}{r},\stackrel{}{w}_\mathrm{H}).`$
Here $`f_\mathrm{H}(\stackrel{}{r},\stackrel{}{w}_\mathrm{H})`$ is the distribution function of H atoms; $`f_p(\stackrel{}{r},\stackrel{}{w}_p)`$ is the local distribution function of protons; $`\stackrel{}{w}_p`$ and $`\stackrel{}{w}_\mathrm{H}`$ are the individual proton and H atom velocities, respectively; $`\sigma _{ex}^{\mathrm{HP}}`$ is the charge exchange cross section of an H atom with a proton; $`\nu _{ph}`$ is the photoionization rate; $`m_\mathrm{H}`$ is the atomic mass; $`\nu _{\mathrm{impact}}`$ is the electron impact ionization rate; and $`\stackrel{}{F}`$ is the sum of the solar gravitational force and the solar radiation pressure force. The plasma and neutral components interact mainly by charge exchange. However, photoionization, solar gravitation, and radiation pressure, which are taken into account in equation (1), are important at small heliocentric distances. Electron impact ionization may be important in the heliosheath (region 2). The interaction of the plasma and H atom components leads to the mutual exchanges of mass, momentum and energy. These exchanges should be taken into account in the plasma equations through source terms, which are integrals of $`f_H(\stackrel{}{r},\stackrel{}{w},t)`$.

### 3.2 Solar wind and interstellar electron and proton components

Basic assumptions necessary to employ a hydrodynamic approach for space plasmas were reviewed in . In particular, it was concluded in the paper that interstellar and solar wind plasmas can be treated hydrodynamically. Indeed, the mean free path of the charged particles in the local interstellar plasma is less than 1 AU, which is much smaller than the size of the heliospheric interface itself. Therefore, the local interstellar plasma is collisional plasma, and a hydrodynamic approach can be used to describe it. Solar wind plasma is collisionless, because the mean free path of the solar wind particles is much larger than the size of the heliopause. Therefore, the heliospheric termination shock (TS) is a collisionless shock. A hydrodynamic approach can be justified for collisionless plasmas when scattering of charged particles on plasma fluctuations is efficient (”collective plasma processes”). In this case, the mean free path $`l`$ with respect to collisions is replaced by $`l_{coll}`$, the mean free path of collective processes, which is assumed to be less than the characteristic length of the problem L: $`l_{coll}<<L`$. However, the integral of ”collective collisions” is too complicated to be used to calculate the transport coefficient for collisionless plasmas.

One-fluid description of heliospheric and interstellar plasmas is commonly used in the global models of the heliospheric interface. However, since measurements of the solar wind show different electron and proton temperatures, two-fluid approach is more appropriate. The temperatures may remain different up to the termination shock and beyond due to the weak energy exchange between protons and electrons.

Hydrodynamic Euler equations for proton and electron components, which take into account the influence of other components such as interstellar H atoms, pickup ions, cosmic rays, electric and magnetic fields, are written below. Mass balance or continuity equations are

$`{\displaystyle \frac{n_s}{t}}+(n_s\stackrel{}{V}_s)=q_{1,s},(s=e,p)`$ (2)
Index $`e`$ denotes electrons, index $`p`$ denotes solar wind protons. $`q_{1,e}=n_H(\nu _{ph}+\nu _{impact})`$, $`q_{1,p}=u\sigma _{ex}^{\mathrm{HP}}(u)f_p(\stackrel{}{w})f_H(\stackrel{}{w}_H)𝑑\stackrel{}{w}𝑑\stackrel{}{w}_H`$ are sources and sinks due to charge exchange, photoionization and electron impact ionization. Here, $`u=|\stackrel{}{w}_H\stackrel{}{w}|`$ is the relative atom-proton velocity, and $`\stackrel{}{w}_H`$ and $`\stackrel{}{w}`$ are individual velocities of H atoms and protons, respectively. Momentum balance equations are

$`{\displaystyle \frac{(n_sm_s\stackrel{}{V}_s)}{t}}+P_s+m_s(n_s\stackrel{}{V}_s\stackrel{}{V}_s)n_se_s(\stackrel{}{E}+{\displaystyle \frac{1}{c}}[\stackrel{}{V}_s\times \stackrel{}{B}])+{\displaystyle \underset{r}{}}\stackrel{}{R}_{sr}=m_s\stackrel{}{q}_{2,s}`$ (3)
where $`\stackrel{}{q}_{2,e}=(\nu _{ph}+\nu _{impact})\stackrel{}{w}_Hf_H(\stackrel{}{w}_H)𝑑\stackrel{}{w}_H`$, $`\stackrel{}{q}_{2,p}=u\sigma _{ex}^{\mathrm{HP}}(u)\stackrel{}{w}_pf_H(\stackrel{}{w}_H)f_p(\stackrel{}{w}_p)𝑑\stackrel{}{w}_H𝑑\stackrel{}{w}_p`$ and $`\stackrel{}{R}_{sr}`$ is the rate of momentum transfer from the $`r`$-gas component to the $`s`$-gas component. The symbol $``$ represents the dyadic product. The momentum transfer term $`\stackrel{}{R}_{sr}`$ can be expressed in a general form through the collision integral $`S_{sr}`$ of kinetic equation of the $`s`$-gas component: $`\stackrel{}{R}_{sr}=m_s\stackrel{}{c}_sS_{sr}𝑑c_s`$, where $`\stackrel{}{c}_s=\stackrel{}{w}_s\stackrel{}{V}_s`$. That $`\stackrel{}{R}_{sr}+\stackrel{}{R}_{rs}=0`$ is a consequence of this definition of $`\stackrel{}{R}_{sr}`$.

Heat balance equations have the following form:

$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{3}{2}}P_s\right)+\left({\displaystyle \frac{3}{2}}P_s\stackrel{}{V}_s\right)+P_s\stackrel{}{V}_s={\displaystyle \underset{r}{}}Q_{sr}+m_sq_{3,s}m_s\stackrel{}{q}_{2,s}\stackrel{}{V}_s`$ (4)
with $`q_{3,e}=\left(\nu _{ph}+\nu _{impact}\right)\frac{\stackrel{}{w}_H^2}{2}f_H(\stackrel{}{w}_H)𝑑\stackrel{}{w}_H`$,

$`q_{3,p}=u\sigma _{ex}^{\mathrm{HP}}(u)\frac{\stackrel{}{w}_p^2}{2}f_H(\stackrel{}{w}_H)f_p(\stackrel{}{w}_p)𝑑\stackrel{}{w}_p𝑑\stackrel{}{v}_H`$. $`Q_{sr}`$ is the heat source due to interactions between between particles of $`s`$ and $`r`$ components. Those terms can be expressed in a general form through the collision term $`S_{sr}`$ of the kinetic equation $`Q_{sr}=\frac{m_sc_s^2}{2}S_{pr}𝑑c_s`$. As a result of this interaction, we have a relation connecting $`Q_{sr}`$ and $`Q_{rs}`$: $`Q_{sr}+Q_{rs}=\stackrel{}{R}_{sr}(\stackrel{}{V}_s\stackrel{}{V}_r)(sr`$).

System (2)-(4) should be added by the state equations: $`P_\alpha =n_\alpha kT_\alpha `$ ($`\alpha =e,p)`$, where $`k`$ is Boltzman constant, and Maxwell equations:

$$\times \stackrel{}{E}=\frac{1}{c}\frac{\stackrel{}{B}}{t};\stackrel{}{E}=4\pi \rho _e;\times \stackrel{}{B}=\frac{4\pi }{c}\stackrel{}{j};\stackrel{}{B}=0$$
(5)
where $`\rho _e`$ is the charge density, $`e`$ is the charge of electron, and $`\stackrel{}{j}`$ is current density. The displacement current has been dropped. Note, that charge and current densities of all charged populations should be taken into account in (5). Neglecting cosmic ray charges and currents we have $`\rho _e=e(n_p+n_{pui}n_e)`$ and $`\stackrel{}{j}=e(n_p\stackrel{}{V}_pn_e\stackrel{}{V}_e+n_{pui}\stackrel{}{V}_{pui})`$. The number density of pickup ions, $`n_{pui}`$, and the bulk velocity of pickup ions, $`\stackrel{}{V}_{pui}`$ are integrals of pickup proton velocity distribution function: $`n_{pui}=f_{pui}(\stackrel{}{w})𝑑\stackrel{}{w}`$, $`\stackrel{}{V}_{pui}=(\stackrel{}{w}f_{pui}(\stackrel{}{w})𝑑\stackrel{}{w})/n_{pui}`$.

Note that in equations (2)-(4) we assume that pickup electrons are indistinguishable from original solar wind electrons, while pickup protons are considered as a separate population.

The expressions for various interaction terms $`\stackrel{}{R}_{sr}`$, $`Q_s`$ ($`s=e,p`$)must be specified. Electron-proton collision terms can be taken in the form given by Braginski :

$$\stackrel{}{R}_{ep}=\frac{m_en_e}{\tau _e}\stackrel{}{u}_e$$
(6)
$$Q_{pe}=Q_{ep}R_{ep}\stackrel{}{u}_e=\frac{3n_e}{\tau _e}\frac{m_e}{m_p}k(T_eT_p)$$
(7)
Here, $`n_e`$ is the electron number density; $`T_e`$ and $`T_p`$ are electron and proton densities, respectively; $`m_e`$ and $`m_p`$ are the electron and proton masses; $`\stackrel{}{u}_e`$ is the electron velocity relative to the proton rest frame. Parameter $`\tau _e`$ characterizes the coupling between electrons and protons and corresponds to the electron collision time for collisional plasma . In collisionless heliospheric plasma, additional assumptions are needed to determine $`\tau _e`$; otherwise, it can be considered as free parameter.

### 3.3 Pickup ions

To study pickup ion dynamical influence on the distant solar wind, the termination shock structure and, finally, on the global heliospheric interface structure, details of the process of charged particle assimilation into the magnetized plasma are needed. A newly created ion under the influence of the steady solar wind electric and magnetic fields executes a cycloidal trajectory with the guiding center, which is drifting at the bulk velocity of the solar wind. Assuming that the gyroradius is much smaller than the typical scale length, one can average velocity distribution function over the gyratory motion. Initial ring-beam distribution of pickup ions is unstable. Basic processes that determine evolution of pickup ion distribution are pitch-angle scattering, energy diffusion in the wave field generated by both pickup ions and the solar wind waves, convection, adiabatic cooling in the expanding solar wind, and injection of newly ionized particles. The most general form of the relevant transport equation to describe the evolution of gyrotropic velocity distribution function $`f_{pui}=f_{pui}(t,\stackrel{}{r},v,\mu )`$ of pickup ions in a background plasma moving at a velocity $`\stackrel{}{V}_{sw}`$ were written in , . $`f_{pui}`$ is a function of the modulus of velocity in the solar wind rest frame, and $`\mu `$ is the cosine of pitch angle.

Complete assimilation of pickup ions into the solar wind would result in a great increase in the temperature with increasing heliocentric distance, which is not observed. Therefore, the solar wind and pickup protons represent two distinct proton populations. Nevertheless, the radial temperature profile of protons measured by Voyager 2 shows a smaller decrease as compared with the adiabatic cooling. A fraction of heating of solar wind protons may be connected with pickup generated waves . Many aspects of pickup ion evolution were studied ( e.g., ; for review, see , ). However, today it still seems to be impossible to take into account all details of the assimilation process of pickup ions into the solar wind in the global models of the heliospheric interface structure. Instead, one may try to use the hydrodynamic approach. In this approach, equations (2)-(4) written for pickup ions represent the balance of their mass, momentum and energy. The right sides of the equations include sources of pickup ions due to ionization processes:

$$q_{1,pui}=n_H\nu _{ph}+u\sigma _{ex}^{\mathrm{HP}}(u)f_H(\stackrel{}{w}_H)f_p(\stackrel{}{w})𝑑\stackrel{}{w}𝑑\stackrel{}{v}_H$$

$$\stackrel{}{q}_{2,pui}=(\nu _{ph}+\nu _{impact})\stackrel{}{v}_Hf_H(\stackrel{}{v}_H)𝑑\stackrel{}{v}_H+u\sigma _{ex}^{\mathrm{HP}}(u)\stackrel{}{w}_Hf_H(\stackrel{}{w}_H)f_p(\stackrel{}{w}_p)𝑑\stackrel{}{w}_H𝑑\stackrel{}{w}_p+$$

$$u\sigma _{ex}^{\mathrm{HP}}(u)(\stackrel{}{w}_H\stackrel{}{w}_i)f_H(\stackrel{}{w}_H)f_{pui}(\stackrel{}{w}_i)𝑑\stackrel{}{w}_H𝑑\stackrel{}{w}_i$$

$$q_{3,pui}=\left(\nu _{ph}+\nu _{impact}\right)\frac{\stackrel{}{w}_H^2}{2}f_H(\stackrel{}{w}_H)f_p(\stackrel{}{w}_p)𝑑\stackrel{}{w}_p𝑑\stackrel{}{w}_H$$

$$+u\sigma _{ex}^{\mathrm{HP}}(u)\frac{\stackrel{}{w}_H^2\stackrel{}{w}_i^2}{2}f_H(\stackrel{}{w}_H)f_{pui}(\stackrel{}{w}_i)𝑑\stackrel{}{w}_i𝑑\stackrel{}{w}_H$$

To complete the model, one should also specify interaction terms $`\stackrel{}{R}_{pui,r}`$, $`Q_{pui,r}`$ ($`rpui`$). The specification of these terms for pickup ion-proton interactions requires analysis of the pickup process in detail at the kinetic level. Global models usually assume immediate assimilation of pickup ions into the solar wind (one-fluid model) or perfect co-moving of these populations $`V_p=V_{pui}`$ and no exchange of energy $`Q_{pui,p}=0`$ (two- or three-fluid models). Energy exchange term of pickup with ACRs, $`Q_{pui,acr}=Q_{acr,pui}`$, is specified in next subsection.

### 3.4 Cosmic rays

The cosmic rays are coupled to background flow via scattering with plasma waves. The net effect is that the cosmic rays tend to be convected along with the background plasma as they diffuse through the magnetic irregularities carried by the background plasma. Both galactic and anomalous cosmic rays can be treated as populations with negligible mass density and sufficient energy density. At a hydrodynamical level, the cosmic rays may modify the wind flow via their pressure gradient $`P_c`$ with the net energy transfer rate from fluid to the cosmic rays given by $`\stackrel{}{V}P_c`$. $`P_c(\stackrel{}{r},t)=\frac{4\pi }{3}_0^{\mathrm{}}f_c(\stackrel{}{r},p,t)wp^3𝑑p`$ is a cosmic ray pressure; $`f_c(\stackrel{}{r},p,t)`$ is the isotropic velocity distribution of cosmic rays.

The transport equation of these particles has the following form :

$$\frac{f_c}{t}=\frac{1}{p^2}\frac{}{p}\left(p^2D\frac{f_c}{p}\right)+(\widehat{k}f_c)\stackrel{}{V}f_c+\frac{1}{3}(\stackrel{}{V})\frac{f_c}{lnp}+S(\stackrel{}{r},p,t)$$
(8)
Here $`p`$ is the modulus of the momentum of the particle; $`D`$ is the diffusion coefficient in momentum space, often assumed to be zero; $`\widehat{k}`$ is the tensor of spatial diffusion; $`\stackrel{}{V}=\stackrel{}{U}+\stackrel{}{V}_{drift}`$ is the convection velocity; $`\stackrel{}{U}`$ is the plasma bulk velocity; $`\stackrel{}{V}_{drift}`$ is a drift velocity in the heliospheric or interstellar magnetic field; and $`S(\stackrel{}{r},p,t)`$ is the source term.

At the hydrodynamic level, the transport equation of the cosmic rays in the heliospheric interface is:

$$\frac{P_c}{t}=[\widehat{k}P_c\gamma _c(\stackrel{}{U}+U_{dr})P_c]+(\gamma _c1)\stackrel{}{U}P_c+Q_{acr,pui}(\stackrel{}{r},t)$$
(9)
Here we assume that $`D=0`$; $`U_{dr}`$ is momentum-averaged drift velocity; $`\gamma `$ is the polytropic index; and $`Q_{acr,pui}`$ is the energy injection rate describing energy gains of the ACRs from pickup ions. Chalov and Fahr (, ) suggested that $`Q_{acr,pui}=\alpha p_{pui}div\stackrel{}{U}`$, where $`\alpha `$ is a constant injection efficiency defined by the specific plasma properties . $`\alpha `$ is set to zero for GCRs since no injection occurs into the GCR component.

## 4 OVERVIEW OF HELIOSPHERIC INTERFACE MODELS

Together with Maxwellian equations (5) the Boltzman equation (1) for interstellar H atoms; sets of hydrodynamic equations (2)-(4) written for solar protons, electrons and pickup ions; and equation (9) written for anomalous and cosmic ray components form a closed system of equations, when interaction terms $`\stackrel{}{R}_{sr}`$, $`Q_{sr}`$ are specified. Possible specification of the interaction terms is given in equations (6), (7). We have to note here, that such a complete model has yet to be developed. However, in recent years, several groups have focused their efforts on theory and modeling in order to understand some effects separately from others. In particular, the influence of the interstellar magnetic field on the interface structure was studied in \-  for the two-dimensional case and in  and  for the three-dimensional case. Both interstellar and interplanetary magnetic fields were considered in  and . A comparison of these MHD models was given recently in . Latitudinal variations of the solar wind have been considered in . The influence of the solar cycle variations on the heliospheric interface was studied in the 2D case in  \-  and in  for the 3D case. In spite of many interesting findings in the papers cited above, these theoretical studies did not take into account the interstellar H atoms, or took them into account but under greatly simplified assumptions, as it was done in , where velocity and temperature of interstellar H atoms were assumed as constants in the entire interface.

Since most of the observational information on the heliospheric interface is connected with interstellar neutrals and their derivatives as pickup ions and ACRs, we will focus on the models, which include interstellar neutrals in a more appropriate way. These models can be separated into two types. Models of the first type (Table 5) use a simplified fluid (or multi-fluid) approach for interstellar H atoms. A kinetic approach was used in the models of the second type. Development of the fluid (or multi-fluid) models of H atoms was connected with the fact that fluid (or multi-fluid) approach is simpler for numerical realization. At the same time such an approach can lead to nonphysical results. Results of one of the most sophisticated multi-fluid models  were compared with the kinetic Baranov-Malama model in . The comparison shows qualitative and quantitative disagreements in distributions of H atoms. At the same time, it was concluded in  that the two models agreed on the distances to the termination shock, heliopause and bow shock in upwind, but not in positions of the termination shock in downwind.

### 4.1 One-fluid plasma models

One of the common features in the models , ,  \- , , , - is that proton, electron and pickup ion components were considered as one fluid. The great advantage of this approach is that its equations are considerably simpler than the three- and two-fluid approaches (see next subsection). A key assumption of this approach is immediate assimilation of pickup protons into the original solar protons. In other words, it is assumed that immediately after ionization one cannot distinguish between original solar protons and pickup protons. Another important assumption is that electron and proton components have equal temperatures, $`T_e=T_p`$. For quasineutral plasma ($`n_p+n_{pui}=n_e+o(n_e)`$) this means that the pressure of the electrons is equal to half of total pressure ($`P=n_ekT_e+(n_p+n_{pui})kT_p2n_ekT_e=2P_e`$). Let us denote total density $`\rho =m_en_e+m_p(n_p+n_{pui})`$, bulk velocity $`\stackrel{}{V}=(_sm_sn_s\stackrel{}{V}_s)/\rho `$ ($`s=e,p,pui`$).

Governing equations for the one-fluid approach can be obtained by summarizing equations (2)-(4) for indexes $`s=e,p,pui`$. Introducing solar wind protons and pickup ions as co-moving and taking into account that $`m_e<<m_p`$ yield one-fluid equations in their general form:

$$\frac{\rho }{t}+(\rho \stackrel{}{V})=q_1,q_1=m_pn_H(\nu _{ph}+\nu _{impact})$$
(10)
$$\frac{(\rho \stackrel{}{V})}{t}+P+(\rho \stackrel{}{V}\stackrel{}{V})\rho _e\stackrel{}{E}\frac{1}{c}[\stackrel{}{j}\times \stackrel{}{B}]=\stackrel{}{q}_2n_em_e\stackrel{}{u}_e\stackrel{}{u}_e$$
(11)
Here $`\stackrel{}{j}=_\alpha n_\alpha e_\alpha \stackrel{}{V}_\alpha `$, $`\rho _e=_\alpha n_\alpha e_\alpha `$; $`\stackrel{}{q}_2=_sm_s\stackrel{}{q}_{2,s}`$; $`\stackrel{}{u}_s=\stackrel{}{V}_s\stackrel{}{V}`$.

$`{\displaystyle \frac{}{t}}\left({\displaystyle \frac{3}{2}}P\right)+\left({\displaystyle \frac{5}{2}}P\stackrel{}{V}\right)\stackrel{}{V}P=\left({\displaystyle \frac{5}{2}}P_e\stackrel{}{u}_e\right)+q_3\stackrel{}{q}_2\stackrel{}{V}+`$ (12)
$`(\stackrel{}{j}\rho _e\stackrel{}{V})\left(\stackrel{}{E}+{\displaystyle \frac{1}{c}}\stackrel{}{V}_e\times \stackrel{}{B}+{\displaystyle \frac{m_e}{en_e}}(\stackrel{}{q}_{2,p}+\stackrel{}{q}_{2,pui})\right)`$
where $`q_3=_sm_sq_{3,s}`$. The relative electron velocity is connected with vector $`\stackrel{}{j}`$ as follows: $`\stackrel{}{j}=\rho _e\stackrel{}{V}en_e\stackrel{}{u}_e`$. Note that to derive (12) we use a generalized form of Ohm’s law and neglect the terms proportional $`m_e/m_p`$ in it:

$$n_ee\left(\stackrel{}{E}+\frac{1}{c}\stackrel{}{V}_e\times \stackrel{}{B}\right)=P_e+\stackrel{}{R}_{ep}+m_e(\stackrel{}{q}_{2,p}+\stackrel{}{q}_{2,pui}\stackrel{}{q}_{2,e})$$

If we assume momentum transfer term $`\stackrel{}{R_{ep}}`$ as in (6) and neglect the term $`\rho _e\stackrel{}{V}`$ in quasineutral plasma, Ohm’s law may be rewritten:

$`\stackrel{}{E}={\displaystyle \frac{1}{c}}[\stackrel{}{V}\times \stackrel{}{B}]+{\displaystyle \frac{m_e}{\tau _en_ee^2}}\stackrel{}{j}+{\displaystyle \frac{1}{en_e}}\left({\displaystyle \frac{1}{c}}\stackrel{}{j}\times \stackrel{}{B}P_e\right)+{\displaystyle \frac{m_e}{en_e}}(\stackrel{}{q}_{2,p}+\stackrel{}{q}_{2,pui}\stackrel{}{q}_{2,e})`$ (13)
To derive the classical system of hydrodynamic equations applied for heliospheric interface in one-fluid models, one needs to ignore terms containing magnetic and electric fields in equations (11) and (12). If we also disregard the second-order term $`n_em_e\stackrel{}{u}_e\stackrel{}{u}_e`$ in (11) and $`(5p_e\stackrel{}{u}_e/2)`$ in (12), equations (10)-(12) together with kinetic equation (1) for H atoms form a closed system of equations.

To get from (10)-(13) a system of ideally conducting MHD equations with a magnetic field frozen in plasma, one obviously needs to make other assumptions in addition to that of high conductivity, when $`R_m>>1`$. $`R_m=4\pi \sigma VL/c^2`$ is magnetic Reynolds number, with the electrical conductivity $`\sigma =n_ee^2\tau _e/m_e`$, and $`V`$ and $`L`$ characteristic velocity and length, respectively. Vanishing electron pressure gradients, the Hall term, $`\stackrel{}{j}\times \stackrel{}{B}`$ and last term of (13) connected with the charge-exchange effect, corresponding terms in the generalized Ohm’s law (13) also must be ignored. Under these conditions, the Ohm’s law has its classical form $`\stackrel{}{E}=\frac{1}{c}[\stackrel{}{V}\times \stackrel{}{B}]`$, and last term of heat flux equation (12) is equal to $`j^2/\sigma `$ and can be neglected. Ideal MHD equations with source terms $`q_1,\stackrel{}{q}_2,q_3`$, on the right-hand sides were considered in , , , .

### 4.2 Three- and two-fluid plasma models

For solar wind, the one-fluid model assumes essentially that wave-particle interactions are sufficient for pickup ions to assimilate quickly into the solar wind, becoming indistinguishable from solar wind protons. However, as discussed above, Voyager observations have shown that this is probably not the case. Pickup ions are unlikely to be assimilated completely. Instead, two co-moving thermal populations can be expected. A model that distinguishes the pickup ions from the solar wind ions was suggested by Isenberg in . Electrons were considered as a third fluid. The key assumption in the model is that pickup ions and solar wind protons are co-moving ($`V_p=V_{pui}`$). It was also assumed that there is no exchange of thermal energy between solar wind protons and pickup ions. Isenberg’s approach consists of two continuity equations (2) for solar protons and pickup ions; one momentum equation (11) and three energy equations (4) for solar wind protons, electrons and pickup ions. In (11) Isenberg neglect term $`\rho _e\stackrel{}{E}`$, the last term and assumes that $`\stackrel{}{j}=c\times \stackrel{}{B}/(4\pi )`$. In energy equations he disregards the energy exchange terms $`Q_{sr}`$. Note that Isenberg used the simplified form of source terms suggested in  and applied these equations to the spherically symmetric solar wind upstream the termination shock.

Another two-fluid approach to model solar wind protons and pickup protons was developed recently in . This model also assumes that convection speed of pickup ions is identical to that of solar wind protons. The pressure of pickup ions is calculated by assuming a rectangular shape of the pickup ion isotropic distribution function. In this case, the pressure can be expressed through pickup ion density $`\rho _i`$ and solar wind bulk velocity $`V_{sw}`$ as

$$P_{pui}=\rho _{pui}V_{sw}^2/5.$$
(14)
Therefore, governing plasma equations are i) one-fluid equations for the mixture of solar protons, pickup ions and electrons; ii) continuity equation for pickup ions; iii) two transport equations for ACRs and GCRs. The influence of cosmic ray components was taken as terms $`(P_{ACR}+P_{GCR})`$ and $`\stackrel{}{V}(P_{ACR}+P_{GCR})\alpha P_{pui}div(\stackrel{}{V})`$ in the right-hand side of the momentum and energy equations, respectively.

## 5  BARANOV-MALAMA MODEL OF THE HELIOSPHERIC INTERFACE

The first self-consistent model of the two-component (plasma and H atoms) LIC interaction with the solar wind was developed by Baranov and Malama . The interstellar wind is assumed to have uniform parallel flow in the model. The solar wind is assumed to be spherically symmetric at the Earth’s orbit. Under these assumptions, the heliospheric interface has axisymmetric structure.

Plasma and neutral components interact mainly by charge exchange. However, photoionization, solar gravity and solar radiation pressure, which are especially important in the vicinity of the Sun, are also taken into account.

Kinetic and hydrodynamic approaches were used for the neutral and plasma components, respectively. The kinetic equation (1 for neutrals is solved together with the Euler equations for one-fluid plasma (10) - (12). The influence of the interstellar neutrals is taken into account in the right-hand side of the Euler equations that contain source terms $`q_1`$, $`\stackrel{}{q}_2`$, $`q_3`$, which are integrals of the H atom distribution function $`f_H(\stackrel{}{V}_H)`$ and can be calculated directly by the Monte Carlo method . The set of kinetic and Euler equations is solved by iterative procedure, as suggested in . Supersonic boundary conditions were used for the unperturbed interstellar plasma and for the solar wind plasma at the Earth’s orbit. The velocity distribution of interstellar atoms is assumed to be Maxwellian in the unperturbed LIC. The model results are discussed below in this section.

### 5.1 Plasma

Interstellar atoms strongly influence the heliospheric interface structure. In the presence of interstellar neutrals, the heliospheric interface is much closer to the Sun than in a pure gas dynamical case (Figure 2). The termination shock becomes more spherical. The Mach disk and the complicated shock structure in the tail disappear.

The supersonic plasma flows upstream of the bow and termination shocks are disturbed. The supersonic solar wind is disturbed by charge exchange with the interstellar neutrals. The new ions created by charge exchange are picked up by the solar wind magnetic field. The Baranov-Malama model assumes immediate assimilation of pickup ions into the solar wind plasma. The solar wind protons and pickup ions are treated as one-fluid, called the solar wind. The number density, velocity, temperature, and Mach number of the solar wind are shown in Figure 3A. The effect of charge exchange on the solar wind is significant. By the time the solar wind flow reaches the termination shock, it is decelerated (15-30 %), strongly heated (5-8 times) and mass loaded (20-50 %) by the pickup ion component.

The interstellar plasma flow is disturbed upstream of the bow shock by charge exchange with the secondary atoms originating in the solar wind and compressed interstellar plasma. Charge exchange results in the heating (40-70 %) and deceleration (15-30 %) of the interstellar plasma before it reaches the bow shock. The Mach number decreases and for a certain set of interstellar parameters ($`n_{H,LIC}>>n_{p,LIC}`$) the bow shock may disappear. Solid curves on Figure 3B correspond to the small ionization degree of LIC ($`n_p/(n_p+n_H)=1/6`$). The bow shock almost disappears.

The interstellar neutrals also modify the plasma structure in the heliosheath. In a pure gas dynamic case (without neutrals) the density and temperature of the postshock plasma are nearly constant. However, the charge exchange process leads to a large increase of the plasma number density and a decrease of its temperature (Figure 3C). The electron impact ionization process may influence the heliosheath plasma flow by increasing the gradient of the plasma density from the termination shock to the heliopause . The effects of interstellar atom influence on the heliosheath plasma flow may be important, in particular, for the interpretations of kHz radio emission detected by Voyager (, ) and possible future heliospheric imaging in energetic neutral atom (ENA) fluxes .

### 5.2 Atoms

Charge exchange significantly disturbs the interstellar atom flow. Atoms newly created by charge exchange have velocities of their ion partners in charge exchange collisions. Therefore, the velocity distribution of these new atoms depends on the local plasma properties. It is convenient to distinguish four different populations of atoms depending on where in the heliospheric interface they originated. Population 1 is the atoms created in the supersonic solar wind. Population 2 is the atoms originating in the heliosheath. Population 3 is the atoms created in the disturbed interstellar wind. We will call original (or primary) interstellar atoms population 4. The number densities and mean velocities of these populations are shown in Figure 4 as the function of the heliocentric distance. The velocity distribution function of interstellar atoms $`f_H(\stackrel{}{w}_H,\stackrel{}{r})`$ can be represented as a sum of the distribution functions of these populations: $`f_H=f_{H,1}+f_{H,2}+f_{H,3}+f_{H,4}`$. The Monte Carlo method allows us to calculate these four distribution functions. The velocity distributions of the interstellar atoms in the 12 selected points in the heliospheric interface were presented in . For example, the velocity distributions at the termination shock in the upwind direction are shown in Figure 5. Note that velocity distributions of H atoms in the heliosphere were also presented in . However, different populations of H atoms cannot be considered separately in mesh particle simulations of H atoms .

Original (or primary) interstellar atoms are significantly filtered (i.e. their number density is reduced) before reaching the termination shock (Figure 4A). Since slow atoms have a smaller mean free path as compared with fast atoms, they undergo more charge exchange. This kinetic effect, called “selection”, results in a deviation of the interstellar distribution function from Maxwellian ( Figure 5A). The selection also results in $``$10 % increase of the primary atom mean velocity to the termination shock (Figure 4C).

The secondary interstellar atoms are created in the disturbed interstellar medium by charge exchange of primary interstellar neutrals and protons decelerated by the bow shock. The secondary interstellar atoms collectively make up the “H wall”, a density increase at the heliopause. The “H wall” has been predicted in  and detected toward $`\alpha `$ Cen . At the termination shock, the number density of the secondary neutrals is comparable to the number density of the primary interstellar atoms (Figure 4A, dashed curve). The relative abundances of the secondary and primary atoms entering the heliosphere vary with degree of interstellar ionization. It has been shown in  that the relative abundance of the secondary interstellar atoms inside the termination shock increases with increasing interstellar proton number density. The bulk velocity of the population 3 is about -18 -19 km/s. The sign “-” means that the population approaches the Sun. One can see that the velocity distribution of this population is not Maxwellian (Figure 5B). The reason for the abrupt behavior of the velocity distribution for $`V_z>0`$ is that the particles with significant positive $`V_z`$ velocities can reach the termination shock only from the downwind direction. The velocity distributions of different populations of H atoms were calculated in  for different directions from upwind. The fine structures of the velocity distribution of the primary and secondary interstellar populations vary with direction. These variations of the velocity distributions reflect the geometrical pattern of the heliospheric interface. The velocity distributions of the interstellar atoms can be a good diagnostics of the global structure of the heliospheric interface.

The third population of the heliospheric neutrals is the neutrals created in the heliosheath from hot and compressed solar wind protons. The number density of this population is an order of magnitude smaller than the number densities of the primary and secondary interstellar atoms. This population has a minor importance for interpretations of Ly $`\alpha `$ and pickup ion measurements inside the heliosphere. However, some of these atoms may probably be detected by Ly $`\alpha `$ hydrogen cell experiments due to their large Doppler shifts. Due to their high energies, the particles influence the plasma distributions in the LIC. Inside the termination shock the atoms propagate freely. Thus, these atoms can be the source of information on the plasma properties in the place of their birth, i.e. the heliosheath .

The last population of heliospheric atoms is the atoms created in the supersonic solar wind. The number density of this atom population has a maximum at $``$5 AU. At this distance, the number density of population 1 is about two orders of magnitude smaller than the number density of the interstellar atoms. Outside the termination shock the density decreases faster than $`1/r^2`$ where $`r`$ is the heliocentric distance (curve 1, Figure 4B). The mean velocity of population 1 is about 450 km/sec, which corresponds to the bulk velocity of the supersonic solar wind. The velocity distribution of this population is not Maxwellian either (Figure 5D). The extended “tail” in the distribution function is caused by the solar wind plasma deceleration upstream of the termination shock. The “supersonic” atom population results in the plasma heating and deceleration upstream of the bow shock. This leads to the decrease of the Mach number ahead of the bow shock.

### 5.3  Recent developments in the Baranov-Malama model

The Baranov-Malama model, the basic results of which were discussed above, takes into account essentially two interstellar components: H atoms and charged particles. To apply this model to space experiments, one needs to evaluate how other possible components of the interstellar medium influence the results of this two-component model. Recently, several effects were taken into account in the frame of this axisymmetric model.

The influence of the galactic cosmic rays on the heliospheric interface structure was studied recently in , . The study was done in the frame of two-component (plasma and GCRs) and three-component (plasma, H atoms and GCRs) models. For the two-component case it was found that cosmic rays could considerably modify the shape and structure of the solar wind termination shock and the bow shock and change the positions of the heliopause and the bow shock. At the same time, for the three-component model it was shown  that the GCR influence on the plasma flows is negligible as compared with the influence of H atoms. The exception is the bow shock, a structure that can be strongly modified by the cosmic rays. It was also found (; Alexashov, private communication) that an anomalous component does not have a significant effect on the position of the termination shock. However, ACRs may significantly reduce compression at the termination shock .

Effects of the interstellar magnetic field on the plasma flow and on distribution of H atoms in the interface were studied in  in the case of magnetic field parallel to the relative Sun/LIC velocity vector. In this case, the model remains axisymmetric. It was shown that effects of the the interstellar magnetic field on the positions of the termination and bow shocks and the heliopause are significantly smaller as compared to model with no atoms . The calculations were performed with various Alfven Mach numbers in the undisturbed LIC. It was found that the bow shock straightens out with decreasing Alfven Mach number (increasing magnetic field strength in LIC). It approaches the Sun near the symmetry axis, but recedes from it on the flanks. By contrast, the nose of the heliopause recedes from the Sun due to tension of magnetic field lines, while the heliopause in its wings approaches the Sun under magnetic pressure. As a result, the region of compressed interstellar medium around the heliopause (or ”pileup region”) decreases by almost 30 $`\%`$, as the magnetic field increases from zero to 3.5 $`\times 10^6`$ Gauss. It was also shown in  that H atom filtration and heliospheric distributions of primary and secondary interstellar atoms are virtually unchanged over the entire assumed range of the interstellar magnetic field (0 - $`3.510^6`$ Gauss). The magnetic field has the strongest effect on density distribution of population 2 of H atoms, which increases by a factor of almost 1.5 as the interstellar magnetic field increases from zero to $`3.510^6`$ Gauss.

Very recently a new non-stationary model of the solar wind interaction with two-component (H atoms and plasma) LIC was proposed in . In this model the primary and secondary interstellar atoms (populations 3 and 4) were treated as quasi-stationary kinetic gases. Population 1 of atoms originating in the supersonic solar wind was considered as zero-pressure fluid. The calculations show that the qualitative features of the non-stationary SW/LIC interaction established in  remain, but the effect of the solar activity cycle is quantitatively stronger because the interface is closer to the Sun than in the model with no atoms. The motion of the termination shock during the solar cycle on the axis of symmetry is about 30 AU. Due to the solar cycle variations of the neutralized solar wind (i.e. atoms of population 1) the region between the heliopause and the bow shock widens and the mean plasma density in the region becomes smaller than for the stationary problem.

## 6 INTERPRETATIONS OF SPACECRAFT EXPERIMENTS ON THE BASIS OF THE BARANOV-MALAMA MODEL

The Sun/LIC relative velocity and the LIC temperature are now well constrained (, -). Using the SWICS pickup He results and an interstellar HI/HeI ratio of 13 $`\pm `$ 1 (the average value of the ratio toward the nearby white dwarfs), Gloeckler et al. concluded that $`n_{LIC}(HI)=0.2\pm 0.03`$ cm -3. This estimation of $`n_{LIC}(HI)`$ is independent of the heliospheric interface model but model-dependent for determination of the number density of H atoms from pickup fluxes. Estimates of interstellar electron number density require a theoretical model of the heliospheric interface. The Baranov-Malama model was used in  to study the sensitivity of the various types of indirect diagnostics of local interstellar plasma density. The diagnostics are the degree of filtration, the temperature and the velocity of the interstellar H atoms in the outer heliosphere (at the termination shock), the distances to the termination shock, the heliopause, and the bow shock, and the plasma frequencies in the LIC, at the bow shock and in the maximum compression region around the heliopause, which constitutes the “barrier” for radio waves formed in the interstellar medium. We also searched  for a number density of interstellar protons compatible with SWICS/Ulysses pickup ion observations, backscattered solar Ly $`\alpha `$ observed by SOHO, Voyager and HST, and kHz radiations observed by Voyager. Table 1 presents the ranges of $`n_{p,LIC}`$ obtained on the basis of the Baranov-Malama model and comparable to these observations.

From analysis of the ranges, it was concluded in  that it is difficult in the frame of the model to reconcile the results obtained from all types of data as they stand now. There is a need for some modifications of the interpretations or of the confidence intervals. Two mutually exclusive solutions have been suggested: (1) It is possible to reconcile the pickup ions and Ly $`\alpha `$ measurements with the radio emission time delays if a small additional interstellar (magnetic or low-energy cosmic ray) pressure is added to the main plasma pressure. In this case, $`n_{\mathrm{p},\mathrm{LIC}}=0.07`$ cm-3 and $`n_{\mathrm{H},\mathrm{LIC}}=0.23`$ cm-3 is the favored pair of interstellar densities. However, in this case, the low frequency cutoff at 1.8 kHz does not correspond to the interstellar plasma density, and one has to search for another explanation. (2) The low-frequency cutoff at 1.8 kHz constrains the interstellar plasma density, i.e., $`n_{\mathrm{p},\mathrm{LIC}}=0.04`$ cm-3. In this case, the bulk velocity deduced from the Ly $`\alpha `$ spectral measurement is underestimated by about 30-50% (the deceleration is about 3 km s-1 instead of 5-6 km s-1). Model limitations (e.g. a stationary hot model to derive the bulk velocity) or the influence of a strong solar Ly $`\alpha `$ radiation pressure may play a role. In this case, a significant additional interstellar (magnetic or cosmic ray) pressure as compared with case (1) would be needed.

This need for an additional pressure is in agreement with the conclusions made in , which were derived from the analysis of the H wall absorption toward alpha Centauri . In their model the authors modified the equation of the state of the gas to simulate the effect of the interstellar magnetic field (IMF) and concluded that H wall absorption favors the “subsonic case”. However, the best model of these authors corresponds to a neutral H density of 0.025 cm-3 in the inner heliosphere, at least 4 times smaller than the density derived from the pickup ions. Also, the precision required to model the differences between the theoretical absorptions, namely small differences of the order of a few kilometers per second at the bottom of the lines, is of the order of the differences between the kinetic and multi-fluid model results for the same parameters in the supersonic case (see Appendix B in  ; ). Thus, an additional study of the absorption toward nearby stars for more realistic densities and models is desired.

## 7 PROBLEMS FOR FUTURE WORK

The Local Interstellar Medium interacts with the solar wind and influences the outer heliosphere in a complicated way. Several particle populations and magnetic fields are involved in this interaction. From the interstellar side, the interacting populations are the plasma (electron and proton) component, H atom component, interstellar magnetic field, and galactic cosmic rays. Heliospheric plasma consists of original solar wind protons, electrons, pickup protons, and the anomalous component of cosmic rays. A large effort has been done to study the theoretical physics of the interaction region. However, a complete, self-consistent model of the heliospheric interface has not yet been constructed, because of the difficulty connecting both the multi-fluid nature of the heliosphere and the requirements of the different theoretical approaches for different components of the interaction. Many aspects were studied and reported here in previous sections. However, some aspects require additional theoretical explorations. Most theoretical models employ the one-fluid approach for solar wind and interstellar plasmas. It has been shown that, to derive one-fluid approach equations, several assumptions are needed. A key assumption that looks reasonable is co-moving character of all components. Another assumption for a one-fluid plasma model is the immediate assimilation of the pickup ion component into the solar wind. As demonstrated by space experiments, this is not the case and it would be more natural to consider solar protons and pickup protons separately as co-moving populations. The electron component should also be treated as a distinct population. However, since the assumption of the co-moving character of these three heliospheric plasma populations looks reasonable, the one-fluid approach gives us a reasonably accurate picture of the flow pattern (positions of the shocks and heliopause) and plasma velocity distributions. Theoretical models of pickup ion acceleration and diffusion can be employed to determine the distribution of thermal energy between solar wind and pickup proton components. A similar study should be done for electrons.

Another important aspect of the solar/wind interaction is a study of the tail region of the solar wind and interstellar medium interaction. Although some studies were done (, ) it is still not clear at which heliocentric distances the gas (plasma and H atoms) parameters become indistinguishable from local interstellar parameters, or in other words, how far signatures of the solar system are noticeable in the interstellar medium. It is still not clear which of the two competiting processes is the most important in the tail region - charge exchange or plasma transport across the heliopause due to different instabilities. Studies of Saturn’s and Earth’s magnetic tails show that such tails can be very extended , .

Finally, growing interest in heliospheric interface studies is connected with expectations that Voyager 1 will cross the termination shock soon. Many predictions of the time of the termination shock crossing by Voyager appeared in the literature. However, it seems that much more work should be done to explain and reconcile all available indirect observations of the heliospheric interface based on the unique model of the heliospheric interface. This work should be done especially because NASA plans to send a spacecraft to a heliocentric distance of at least 200 AU with a flight-time of only 10 or 15 years. Intensive theoretical study will help to optimize goals, instrumentation, and, finally, the scientific profit of this ”interstellar” mission.

Acknowledgements. This work was supported in part by CRDF Award RP1-2248, INTAS Award 2001-0270, YSF 00-163, RFBR grants 01-02-17551, 02-02-06011, 01-01-00759, and the International Space Science Institute in Bern. I thank V. B. Baranov and S. V. Chalov for useful discussions.


format: text

text:
# RXJ1856-3754: Evidence for a Stiff EOS

## 1 Introduction

$`\mathrm{RX}\mathrm{J1856}3754`$, discovered by Walter, Wolk, & Neuhäuser (1996), is the nearest and brightest known neutron star not showing emission dominated by non-thermal magnetospheric processes. As such it offers a unique opportunity to study the bare thermal surface emission. Measurements of the spectrum can probe the neutron star mass ($`\mathrm{M}_{\mathrm{NS}}`$) and radius ($`\mathrm{R}_{\mathrm{NS}}`$), constraining the high density equation of state (EOS). Since the discovery, there have been several intensive observing campaigns covering the optical-UV (HST) and most recently the detailed soft X-ray spectrum (CXO). An initial $`50\mathrm{ks}`$ Low-Energy Transmission Grating (LETG) spectrum showed a broad band spectrum remarkably consistent with a simple blackbody (Burwitz et al., 2001), although hints of spectral features were suggested (van Kerkwijk, 2002). A deeper observation using $`450\mathrm{ks}`$ of Director’s Discretionary Time (DDT) was made. This unique data set has been the subject of prompt study; several authors show that lines in the spectrum are undetectable, while pulse searches have placed stringent limits on the observed soft X-ray pulse fraction (Ransom, Gaensler, & Slane, 2002; Drake et al., 2002). These data have been variously interpreted, including the widely reported speculation (based on the X-ray spectrum alone) that $`\mathrm{RX}\mathrm{J1856}3754`$ might be a bare quark star (esp. Drake et al., 2002).

Despite the very stringent constraints placed on the X-ray pulse fraction, ($`<4.5\%`$ at $`99\%`$ confidence, including a $`\dot{P}`$ search Ransom, Gaensler, & Slane, 2002), there is strong evidence that $`\mathrm{RX}\mathrm{J1856}3754`$ is a rotation-powered pulsar. van Kerkwijk & Kulkarni (2001b) discovered an H$`\alpha `$ nebula surrounding the neutron star, concluding that it could be best interpreted as a bow-shock nebula powered by a relativistic wind of $`e^\pm `$ generated by pulsar spindown. The bow shock geometry then provides an estimate of the spindown power $`\dot{E}=I\mathrm{\Omega }\dot{\mathrm{\Omega }}=8\times 10^{32}\mathrm{erg}/\mathrm{s}d_{140}^3`$ (van Kerkwijk & Kulkarni, 2001b). Adopting magnetic dipole braking at constant $`B`$, this gives $`\dot{E}=10^{34}(B_{12}\tau _6)^2\mathrm{erg}/\mathrm{s}`$ for a surface dipole field $`10^{12}B_{12}\mathrm{G}`$ and characteristic age $`10^6\tau _6\mathrm{y}`$, suggesting $`B_{12}\tau _63`$.

A critical parameter in the discussion of this source is the distance, which has been the subject of some controversy. Initial estimates from HST astrometry gave a parallax distance of $`61\mathrm{pc}`$ (Walter, 2001). Kaplan, van Kerkwijk, & Anderson (2002), however re-analyzed these data, deriving $`d=143\mathrm{pc}`$. A fourth HST observation appears to have resolved this discrepancy, giving an overall measurement of $`d=117\pm 12\mathrm{pc}`$ (Walter & Lattimer, 2002); we adopt this value.

## 2 Spectral Fits

For some time now, it has been clear that the broad band spectrum of $`\mathrm{RX}\mathrm{J1856}3754`$ from HST, ROSAT, and EUVE data (for a detailed discussion, see Pons et al., 2002) is inconsistent with a light element ($``$ Kramer’s law opacity) H or He atmospheres. H models, for example, overpredict the optical/X-ray flux ratio by a factor $`100`$ (Pavlov et al., 1996). Blackbody, heavy element or composite models gave acceptable fits. To produce such $``$ Planckian spectra, one must have nearly isothermal conditions at optical depth $`\tau _\nu 1`$ across the observed band. One possibility is that the surface is in a solid or liquid state, precluding large temperature gradients through the photosphere. Theoretical studies to date, while limited, suggest that this is not the case unless H is present and $`B>10^{13}\mathrm{G}`$ (Lai & Salpeter, 1997). In an atmosphere, the easiest way to form the spectrum in a small depth range is to invoke a rich line spectrum; thus, at first sight, blackbody-like spectra should require a high spectral density of opacity features (lines and edges). This led to the expectation that the blackbody-like spectrum of $`\mathrm{RX}\mathrm{J1856}3754`$ would show many spectral features when examined with good S/N LETG spectra. Unfortunately, the initial $`50\mathrm{ks}`$ exposure already showed that the line features were substantially weaker than expected in a simple low-$`B`$ single temperature atmosphere model dominated by heavy elements. The new DDT exposure only enhances this conclusion, placing very strong limits (typical equivalent width $`0.02`$Å Drake et al., 2002) on spectral features.

We will be concerned with the quality of statistical fits to various atmosphere models. It is important to note here that blackbody spectral fits to the CXO $`\mathrm{RX}\mathrm{J1856}3754`$ LETG data are, contrary to early reports, *not* statistically perfect fits to a simple Planck spectrum. This conclusion was drawn from fits to basic CIAO extractions, which appreciably underbin the spectrum. Instead we find that at a more modest binning (equally spaced $`0.7`$Åbins) we obtain $`\chi ^2/`$DOF=1.6. As the bin size is increased, the $`\chi ^2/`$DOF grows to $`4.8`$, until the number of degrees of freedom becomes small. This is a clear signature of spectral departures on resolved energy scales, and with appropriate binning one indeed finds systematic, grouped residuals to the Planck function fit at the $`10`$% level. We believe that these represent the limit of accuracy in calibration of the response matrix, as the broad band spectral shape is an excellent fit to the Planck function. Drake et al. (2002) have reached similar conclusions. Recognizing that very subtle departures from a pure blackbody may be present in these data, we adopt the conservative assumption that these departures are fully accounted for by response matrix systematics.

To accommodate an extended atmosphere, one must suppress the spectral features. One possibility is that external heat sources (such as precipitating magnetosphereic $`e^\pm `$) drive the atmosphere towards isothermality. Sample atmospheres showing this effect have been computed in Gänsicke, Braje, & Romani (2002). A second possibility is that the line energies for a given species vary strongly across the neutron star surface. For normal pulsar fields, $`B10^{12}\mathrm{G}`$ or higher, the strong dependence of the transition energies on the local $`B`$ (eg. Rajagopal, Romani, & Miller, 1997; Pavlov et al., 1995), coupled with substantial $`2\times `$ variation of $`B`$ across the surface, even for the simplest dipole models, ensures that such ‘magnetic smearing’ will strongly suppress the phase-averaged line width (Romani, in prep.). We discuss briefly here a third possibility, that the lines experience variable shifts as the pulsar rotates due to Doppler and other dynamic effects. Again, the phase-averaged spectrum observed for $`\mathrm{RX}\mathrm{J1856}3754`$ would be expected to show broadened and blended lines, driving the spectrum towards a Planck curve.

If, as is required for a normal neutron star, the soft X-ray emission of $`\mathrm{RX}\mathrm{J1856}3754`$ is dominated by hot polar caps, the rotation of these past the line of sight produces phase dependent Doppler shifts (Braje, Romani, & Rauch, 2000). These are only significant for $`v_{surf}2\pi \mathrm{R}_{\mathrm{NS}}/\mathrm{P}_{\mathrm{NS}}\mathrm{c}`$, i.e. $`\mathrm{P}_{\mathrm{NS}}`$ a few ms. Such small $`\mathrm{P}_{\mathrm{NS}}`$ are not excluded, since the HRC-S wiring error limits arrival time accuracy and precludes sensitive searches for $`\mathrm{P}_{\mathrm{NS}}10\mathrm{ms}`$. Moreover, the $`\dot{E}`$ from the bow shock standoff gives $`\mathrm{P}_{\mathrm{NS}}=4.6(\mathrm{B}/10^8\mathrm{G})^{1/2}\mathrm{ms}`$, so a low field star, having a non-magnetic atmosphere would have a $`\mathrm{ms}`$ spin period.

For a concrete example, we assume $`\mathrm{M}_{\mathrm{NS}}=1.4\mathrm{M}_{}`$, $`\mathrm{R}_{\mathrm{NS}}=10\mathrm{km}`$, and $`\mathrm{P}_{\mathrm{NS}}=1.5\mathrm{ms}`$ (allowed, as the existence of $`PSR1937+21`$ shows). We have tried both solar abundance and iron model atmospheres, tested a range of magnetic inclinations $`\alpha `$, and computed phase averaged spectra using an extension of the Monte Carlo simulation code described in Braje, Romani, & Rauch (2000). We then tested these models, fitting the most recent CXO data, allowing the temperature, observer viewing angle $`\zeta `$, and interstellar absorption $`\mathrm{n}_\mathrm{H}`$ to vary. The large L and M shell edges in the iron models ensure that these are always poor fits. The solar abundance atmospheres have a much richer line structure which is more easily blurred into a pseudo-continuum. In 2, we display the solar abundance millisecond pulsar model, overlaid on the CXO data. For comparison the best fit non-rotating model is shown in the upper panel. Doppler boosting produces a qualitatively acceptable fit below $`0.5\mathrm{keV}`$, but the simple blackbody remains statistically an appreciably better model ($`\chi ^2`$/DOF=1.6 vs. 3.7 for the binning chosen). We must conclude that a simple blackbody fit remains the best available, although not yet physically explained. Planck emission from a physical neutron star can be compared with the data to obtain significant constraints on the neutron star parameters; we pursue this in the remainder of the paper.

## 3 Two Temperature Blackbody Fits

A blackbody fit to the CXO data alone results in a temperature of $`T_{\mathrm{}}61\mathrm{eV}`$ with very small statistical uncertainty (Drake et al., 2002, and our own analysis). We find, as also reported by Drake et al. (2002), that systematics (likely in the the effective area, as noted above) provide the dominant error. Drake et al. (2002) quote $`61.2\pm 1.0\mathrm{eV}`$. Taken at face value this $`T_{\mathrm{}}`$ with the parallax distance gives a radius as measured at infinity of $`\mathrm{R}_{\mathrm{}}=(1+\mathrm{z})\mathrm{R}_{\mathrm{NS}}=3.88.2\mathrm{km}`$. If interpreted as the full star radius, this demands exotic equations of state (i.e. quark stars). Of course, this is only a lower limit to the stellar radius. The optical/UV data points, which closely follow a Rayleigh-Jeans spectrum (eg. Pons et al., 2002; van Kerkwijk & Kulkarni, 2001a), are most easily interpreted as a second cooler Planck spectrum representing flux from the full surface. This has been previously recognized, but Drake et al. (2002) argue against this interpretation, citing the absence of the X-ray pulse expected from such a hot polar cap/cool surface combination. We have addressed this concern quantitatively, computing detailed light curves and spectra.

### 3.1 Analytic Two-Temperature Models

A simple analytic two-temperature blackbody fit delineates the basic model parameters. For a range of effective (cold surface) radii, we fit $`\mathrm{T}_{\mathrm{hot}}`$, $`\mathrm{T}_{\mathrm{cold}}`$, the hot area, and the absorption column density. In Figure 2, we plot the minimum $`\chi ^2`$ as a function of cold radius. The optical-UV data points fix $`\mathrm{R}_{\mathrm{NS}}^2\mathrm{T}_{\mathrm{cold}}`$. The fit becomes poor when $`\mathrm{T}_{\mathrm{cold}}`$ starts to allow significant Wien peak contribution to the CXO X-ray band; this sets the minimum stellar radius. Formally, there is a maximum acceptable radius beyond which low $`\mathrm{T}_{\mathrm{cold}}`$ predicts Wien peak curvature inconsistent with the $``$Rayleigh-Jeans UV data points.

### 3.2 More Realistic Two-Temperature Models

In addition to the phase averaged spectrum, we have a limit on the CXO-band pulse fraction. A more realistic model is required to address the detailed spectrum and pulsations. We adopt a two-temperature model with two opposing hot spots (polar caps) at $`\mathrm{T}_{\mathrm{hot}}`$ and the remainder of the surface at $`\mathrm{T}_{\mathrm{cold}}`$. The caps’ orientation ($`\alpha `$ and $`\zeta `$) are free parameters. We radiate from these surfaces, tracing the photons to infinity to form phase resolved spectra and light curves. For details of these Monte-Carlo sums see Braje, Romani, & Rauch (2000).

The analytic model results allow some useful simplifications. Since $`\mathrm{T}_{\mathrm{hot}}`$ is virtually constant over the full acceptable $`\mathrm{R}_{\mathrm{NS}}`$ range, we fix this value in exploring the rest of the parameter space. Further, the X-ray flux amplitude allows an initial estimate of the cap half angle $`\mathrm{\Delta }`$ (which depends on $`\alpha `$ and $`\zeta `$) that ensures that the pulse formation is accurate and results in quick convergence to the true minimum.

The fit parameters of interest are $`\mathrm{R}_{\mathrm{NS}}`$, $`\mathrm{\Delta }`$, $`\mathrm{T}_{\mathrm{cold}}`$, $`\zeta `$, $`\alpha `$, and $`\mathrm{n}_\mathrm{H}`$. To explore the sensitivity to several parameters, we have computed a model grid. We calculate all models for $`\alpha =5\mathrm{°}`$ to $`\alpha =90\mathrm{°}`$ in five degree steps; $`\zeta =0\mathrm{°}`$ to $`\zeta =90\mathrm{°}`$ in five degree steps; and $`\mathrm{R}_{\mathrm{NS}}=12\mathrm{km}`$ to $`\mathrm{R}_{\mathrm{NS}}=20\mathrm{km}`$ in one kilometer steps. While not strictly a fit parameter, we also vary the stellar mass $`\mathrm{M}_{\mathrm{NS}}`$. The quality of the spectral fit turns out to be quite insensitive to the choice of $`\alpha `$ and $`\zeta `$. In Figure 3, we display a typical spectrum for a model with $`\mathrm{M}_{\mathrm{NS}}=1.4\mathrm{M}_{}`$, $`\mathrm{R}_{\mathrm{NS}}14\mathrm{km}`$. As the radius becomes smaller, the optical-UV flux is under-predicted and the fit becomes unacceptable.

For each $`\alpha `$ and $`\zeta `$ we compute $`\chi ^2`$ as a function of $`\mathrm{R}_{\mathrm{NS}}`$. The average over angles is shown by the points in Figure 2. Note that the minimum value is not at $`\chi _\nu ^2=1`$, a consequence of the aforementioned systematic errors. This means that the errors are not Poisson distributed according to the bin counts. To establish confidence levels (CL), we must Monte Carlo according to the observed error distribution. We have computed the $`\chi ^2`$ distribution about the best fit model for each ($`\alpha `$, $`\zeta `$) combination, obtaining a histogram of $`\chi ^2`$ values. These were almost completely insensitive to the angles, so we combined all the $`\chi ^2`$ distributions to obtain confidence levels to the $`\chi ^2`$ increases associated with variations in $`\mathrm{R}_{\mathrm{NS}}`$. Examinations of the differences between the different angles and between independent Monte Carlo runs show that the estimates of 90% and 99% confidence level limits on the radius are uncertain by no more than $`\pm 0.5\mathrm{km}`$ from systematic and computational errors in this procedure.

### 3.3 Cap Shape Constraints

One might question whether a simple, circular uniform $`\mathrm{T}_{\mathrm{hot}}`$ polar cap is merely an adequate approximation to the data. We have fit some alternative surface temperature distributions $`T(\eta )`$, where $`\eta `$ is the magnetic co-latitude, comparing with the best fit circular cap which had $`k\mathrm{T}_{\mathrm{hot}}=62.8\mathrm{eV}`$, $`\mathrm{\Delta }=21\mathrm{°}`$. The best fit Gaussian $`T(\eta )`$ had a peak temperature of $`k\mathrm{T}_{\mathrm{hot}}=74.8\mathrm{eV}`$ and width $`\sigma =19\mathrm{°}`$, but showed an increase of $`\chi ^280`$ over the simple cap model, sufficient to exclude at the $`99`$% CL. Adding a simple linear $`\mathrm{T}_{\mathrm{hot}}`$ to $`\mathrm{T}_{\mathrm{cold}}`$ ramp to a uniform cap makes no discernible difference until the ramp width is twice that of the cap. At this point the best fit model has a cap with $`k\mathrm{T}_{\mathrm{hot}}=69\mathrm{eV}`$ and $`\mathrm{\Delta }12\mathrm{°}`$, but the model is excluded at the $`90`$% CL. Finally, we fit both a simple $`\mathrm{T}_{\mathrm{hot}}\mathrm{cos}(\eta )`$ model and the surface $`T(\eta )`$ distribution of Greenstein & Hartke (1983) which is motivated by magnetic anisotropy in the thermal conductivity. Both were excluded by the CXO data at very high confidence. All fits were at $`1.4\mathrm{M}_{}`$, best fit stellar radii and ($`\alpha `$, $`\zeta `$) chosen so that the pulse fraction is $`5\%`$ for the default cap model.

Evidently, the CXO data require a quite uniform distribution of the high temperature excess, and suggests that it is induced by exterior heating rather than interior conductivity. This cap size is substantially larger than the $`\mathrm{\Delta }3\mathrm{°}(100\mathrm{ms}/\mathrm{P}_{\mathrm{NS}})^{1/2}`$ expected for a dipole surface cap, unless the period is very small. Higher magnetic multipoles would generally have even smaller open zone caps. The large, uniformly heated area is puzzling in the context of pulsar surface acceleration models, but might be most easily accommodated in the more modern pictures of a GR-induced potential that is relatively uniform across the polar cap and forms a pair formation front at relatively high altitudes (Harding & Muslimov, 1998). Perhaps a more plausible interpretation invokes a high-altitude acceleration zone, with the inward-directed $`\gamma `$-rays pair converting in the closed zone above the polar cap (Wang et al., 1998), which should give a $``$ large zone of uniform surface heating.

## 4 Pulse Fraction Constraints

We have seen that the spectral fits are quite insensitive to the cap viewing angles. They place a firm minimum on the allowed $`\mathrm{R}_{\mathrm{NS}}`$ but allow radii that are implausibly large. However, the observed pulse fraction depends strongly on viewing geometry and, through gravitational focusing, the value of $`\mathrm{M}_{\mathrm{NS}}/\mathrm{R}_{\mathrm{NS}}`$. In general, we expect a strong X-ray pulse when the magnetic axis passes close to Earth’s line-of-sight (small $`|\alpha \zeta |`$). Of course, for an aligned rotator $`\alpha 0`$ the pulse can be very small and for $`\alpha \pi /2`$ there is an appreciable region where the pulse is weak as viewed from Earth. We have examined the CXO-band pulse profile of our model grid to find the allowed region of $`(\alpha ,\zeta )`$ parameter space in which the model pulse fraction is weaker than that observed ($`4.5\%`$); this is shown in the inset of 4.

This exercise is repeated for each $`\mathrm{M}_{\mathrm{NS}},\mathrm{R}_{\mathrm{NS}}`$; the allowed phase space is larger for small $`\mathrm{R}_{\mathrm{NS}}/\mathrm{M}_{\mathrm{NS}}`$ as gravitational bending dilutes the observed pulse. By demanding a certain minimum probability that a pulse fraction as low as that observed is seen at Earth, we obtain a maximum acceptable radius for the neutron star. One important additional constraint can be invoked. The non-detection of this neutron star in the radio band (Brazier & Johnston, 1999) puts a very strong bound on pulsar emission directed towards Earth. Assuming that this is a normal radio pulsar (consistent with $`\dot{E}`$ inferred from the bow shock), we must conclude a priori (independent of the X-ray data) that our line of sight lies outside of the pulsar radio beam and that $`|\alpha \zeta |`$ is not small. The region excluded depends on the size of the radio beam ($`\mathrm{\Theta }_{\mathrm{rad}}`$), which in turn depends on the spin period. A typical estimate is (eg. Rankin, 1993)

$$\mathrm{\Theta }_{\mathrm{rad}}=5.8\mathrm{°}(\mathrm{P}_{\mathrm{NS}}/1\mathrm{s})^{1/2}$$
(1)
at the radio frequency $`\nu =1`$ GHz. With $`\mathrm{P}_{\mathrm{NS}}=0.3\mathrm{s}`$, this gives $`\mathrm{\Theta }_{\mathrm{rad}}10.6\mathrm{°}`$. Radio limits on $`\mathrm{RX}\mathrm{J1856}3754`$ are strong at even lower $`\nu `$ where radius-to-frequency mapping produces larger $`\mathrm{\Theta }_{\mathrm{rad}}`$, and for this nearby object, the radio luminosity constraints are so severe that we are unlikely to intersect even the faint fringe of the radio beam. Both effects argue for $`\mathrm{\Theta }_{\mathrm{rad}}`$ larger than that above. Accordingly, in the 4 inset we show by diagonal lines the regions near $`\alpha =\zeta `$ excluded by the lack of radio detection for $`1\times `$ and $`2\times `$ this fiducial beam size. Shorter periods allow an a priori exclusion of even more phase space. The point is that from the lack of radio detection we should have expected a low X-ray pulse fraction. The probability for obtaining pulses as weak as observed is then set by the fraction of the remaining allowed solid angle. These fractions are plotted in 4 with and without the radio prior. The allowed radius range depends on mass and in Table 1, we give the 90% and 99% CL upper bound on the neutron star radius for these priors and several neutron star masses.

## 5 Equation of State Constraint

We see that the X-ray/optical data, using the spectral and pulse fraction arguments, give a range of allowed stellar radii for each mass. Strictly speaking, the minimum and maximum radius CL have somewhat different interpretations, but it is interesting to place these bounds in the $`\mathrm{M}_{\mathrm{NS}}\mathrm{R}_{\mathrm{NS}}`$ plane to compare with the predictions of various EOS.

In 5, we show the combined constraints, assuming $`\mathrm{\Theta }_{\mathrm{rad}}=21.2\mathrm{°}`$. The spectral radius lower bounds at 90% and 99% CL approximate curves of constant $`\mathrm{M}_{\mathrm{NS}}/\mathrm{R}_{\mathrm{NS}}`$. The pulse fraction upper bounds (90%, 95%) rapidly drive one to small radii at low masses. For $`\mathrm{M}_{\mathrm{NS}}1.3\mathrm{M}_{}`$ no simultaneous solutions are allowed consistent with the 90% bounds. At $`\mathrm{M}_{\mathrm{NS}}=1.5\mathrm{M}_{}`$ the allowed range from the fit is quite small, $`\mathrm{R}_{\mathrm{NS}}=13.7\pm 0.6\mathrm{km}`$. The additional uncertainty in the distance actually dominates the errors (arrowed bar).

For comparison, several EOS curves (after Lattimer & Prakash, 2001) are plotted. We see that large radius (stiff at nuclear density) EOS are preferred. Formally, the relativistic field theoretical model by Müller & Serot (1996) and the model GS2 by Glendenning & Schaffner-Bielich (1999) are the only modern models allowed (the original PS model of Pandharipande & Smith (1975) is also allowed). We note that the GS models are very sensitive to the K meson potential; for example, GS1 is strongly excluded. Interestingly, no potential or variational method computations agree with the formal overlap. If one includes the distance uncertainty, a few more intermediate radius models are not excluded at the 90% CL. However, even including the distance uncertainties, all quark star models are excluded at the $`95`$% level for $`\mathrm{M}_{\mathrm{NS}}1.5\mathrm{M}_{}`$ and are only barely consistent at the highest allowed masses. Improved HST parallax measurements could boost this exclusion to the 3$`\sigma `$ level.

The model GS2 differs from the MS models in that it includes a kaon condensate in the core. One principal effect of exotic interior condensates is to enhance neutrino cooling after $`100\mathrm{y}`$. This can also be achieved without exotica for a proton fraction $`11\%`$ through the direct URCA process (Lattimer et al., 1991). Our inferred $`\mathrm{T}_{\mathrm{cold}}`$, interpreted as the signature of the cooling from the initial heat of formation, corresponds to $`L3.2\times 10^{31}\mathrm{erg}/\mathrm{s}`$ for our best fit radius. Such luminosities are reached in $`5\times 10^5\mathrm{y}`$ (the preferred age of Walter & Lattimer, 2002) when enhanced neutrino cooling can occur, but are achieved after $`1.3\times 10^6\mathrm{y}`$ for stars with low density (stiff) cores Tsuruta et al. (2002). Thus, depending on the actual stellar age, this cool surface may be seen as weak evidence for an exotic composition and/or significant softening of the EOS at very high densities.

The prospects for further tests of the ideas in this paper hinge on the detection of a pulse from $`\mathrm{RX}\mathrm{J1856}3754`$ along with measurement of the pulse fraction and the period derivative. Given that significant allowed pulse fraction parameter space lies just slightly below the CXO LETG detection threshold, the prospects for a pulse measurement with the recently completed $`58\mathrm{ks}`$ XMM observation are quite good. Even if only thermal, the phase resolved spectrum should provide important constraints on the cap temperature, size and orientation. In particular our thermal model predicts that the pulse fraction should increase by 70% from $`0.15\mathrm{keV}`$ to $`0.5\mathrm{keV}`$, where the XMM data should still deliver $`0.8`$ PN camera cps, allowing detection of pulse fractions well under 1%.

We are grateful to Boris Gänsicke for collaboration on the atmosphere models used in §2; and to Herman Marshall for sharing an independent LETG response matrix. This work was supported in part by NASA grant SP2-2002X.