text:
# Surface Emission Properties of Strongly Magnetic Neutron Stars

## 1 INTRODUCTION

The magnetized surface layers of young neutron stars strongly affect their observable characteristics, as well as many physical processes that take place in their interiors. Understanding the properties of these layers becomes most important when the questions under study are related to the transport of thermal energy from an internal source to the surface of neutron stars, as the magnetized envelopes determine the global thermodynamic properties of the neutron star in addition to generating the thermal emission spectra. Strong magnetic fields ($`B10^{12}`$ G) significantly modify the plasma properties and the transport of radiation through the neutron star atmospheres. Thermal emission from neutron stars with such magnetic fields has been observed as optical and soft X-ray emission from young radio pulsars (e.g., Ögelman 1995) and has also recently been suggested as the origin of the X-ray emission from two new classes of X-ray pulsars.

Neutron stars with ultrastrong magnetic fields ($`B10^{14}`$ G), or magnetars (Thompson & Duncan 1993), currently provide a very promising explanation for these two classes of sources, the Anomalous X-ray Pulsars (AXPs; Mereghetti & Stella 1995; see Mereghetti 2000 for a review) and the Soft Gamma-ray Repeaters (SGRs; e.g., Hurley 2000). According to the magnetar models, the quiescent X-ray emission from SGRs and the persistent pulsed X-ray emission of the AXPs may be powered by internal (crustal) heating arising from the decay of such strong magnetic fields, thought to be generated by dynamo processes in proto-neutron stars (Usov 1992; Thompson & Duncan 1993). Alternatively, strong magnetic fields direct the internal heat of the young neutron stars preferentially to the magnetic poles generating an anisotropic surface temperature distribution and may therefore give rise to the pulsed emission of SGRs (Usov 1997) and AXPs (Heyl & Hernquist 1997). In addition, magnetospheric processes may contribute significantly to the total emission from magnetars (e.g., Thompson & Duncan 1996) and have often been suggested as the explanation for the hard spectral tails observed in these sources.

Recent observations have provided support for the applicability of the magnetar models to AXPs and SGRs. The isolated nature of these sources, suggested by the lack of observable companions (Mereghetti, Israel, & Stella 1998; Hulleman, van Kerkwijk, & Kulkarni 2001), as well as their apparent youth, indicated by the potential supernova-remnant associations (Kaspi 2001) and their proximity to the Galactic plane, are naturally accounted for by magnetar models. Some of the timing properties, such as the large period derivatives and noisy spin-down (Woods et al. 1999; Thompson et al. 2000; Kaspi, Lackey, & Chakrabarty 2000) may also be related to the presence of an ultrastrong magnetic field.

In the absence of calculations addressing the angle- and energy-dependent emission from the surface of ultramagnetized neutron stars, however, it has not been possible to date to perform similar detailed comparisons of models to observable properties such as the pulse profiles and spectra. Analysis of folded lightcurves and spectral modeling so far has had to rely on simple estimates and parameterizations. Significant qualitative changes in the transport of radiation in the presence of strong magnetic fields prohibit an extrapolation of results from model atmospheres obtained for nonmagnetic or weakly-magnetic neutron stars. Furthermore, the magnetic opacities have a very strong dependence on photon energy and on the angle with respect to the magnetic field direction; thus, angle- and energy-integrated approaches to this problem are insufficient (see also §3). Motivated by this, we have carried out radiative transfer calculations in magnetar atmospheres, for magnetic fields in the range $`10^{13}B10^{15}`$ G, which is a previously unexplored regime.

Significant progress has been made recently by several authors in modeling neutron star atmospheres with typical ($`10^{12}`$ G) or negligible magnetic fields. These calculations have been carried out in the context of weakly-magnetic neutron stars with H atmospheres (Shibanov et al. 1992; Pavlov et al. 1994; Zavlin, Pavlov, & Shibanov 1996), as well as for magnetized heavy-element atmospheres (Rajagopal & Romani 1996; Rajagopal, Romani, & Miller 1997) with varying degrees of approximation. The coherent scattering of radiation has so far been included only in the case of weakly magnetized atmospheres with isotropic opacities (Zavlin et al. 1996). In the case of magnetized atmospheres with anisotropic opacities, angle-averaged solutions have been obtained using the moment equations (Rajagopal & Romani 1996, Shibanov et al. 1992, Pavlov et al. 1994) or scattering has been neglected altogether (Rajagopal et al. 1997).

The current work represents, to our knowledge, the first solution of the strongly magnetized radiative equilibrium atmosphere problem with angle-dependent scattering and the inclusion of vacuum polarization effects. As the increasing magnetic field strength renders the angle and energy dependence of the photon-electron interaction cross sections more severe, as well as qualitatively altering the radiation transport properties because of vacuum polarization effects, this new regime requires a new set of tools for the radiative transfer problem. The description of these methods as well as the results from the model atmospheres is the subject of this paper.

The important effects of the resonance due to the polarization of the magnetic vacuum on radiation spectra has previously been addressed by Bulik & Miller (1997). They show that broad-band absorption features appear due to the enhanced interaction cross sections when high energy photons propagate through a strongly magnetized plasma. Zane et al. (2000) also present a solution of the ultramagnetized problem with scattering but suffer from the use of ordinary lambda iteration method for solution of the scattering which is known to have convergence problems (see, e.g., Mihalas 1978). After the submission of the present paper, two independent treatments of ultramagnetized atmospheres have appeared. Ho & Lai (2001) and Zane et al. (2001) include the proton cyclotron line but neglect the vacuum polarization resonance. It should also be noted that Ho & Lai (2001) do not retain the full angle dependence, while Zane et al. (2001) employ a local temperature correction scheme for the radiative equilibrium atmosphere, which does not yield radiative equilibrium solutions.

The thermodynamic properties of magnetized crusts and envelopes have been addressed previously (e.g., Heyl & Hernquist 1998; Potekhin & Yakovlev 2001). It has been shown that heat transport is affected not only by the surface magnetic-field strengths, but by the chemical composition of the neutron star atmospheres as well. In particular, large fluxes through the neutron staratmosphere can be obtained only in the presence of a light-element (H and He) atmosphere, because the heat conductivity is inversely proportional to the atomic charge in strong magnetic fields. Furthermore, in the case of observations of soft X-ray emission from neutron stars and thermal optical emission from young radio pulsars, spectral fits have required the presence of a H layer on the neutron star (Pavlov et al. 2001). It is for these reasons that, in the current work, we restrict our attention to the transport of radiation through a H plasma.

In §2 we outline the physical processes in a strongly magnetized plasma and the equations that describe a magnetized atmosphere in radiative equilibrium. In §3 we present the method of solution, in §4 we discuss the energy and angle dependence of the radiation emerging from the neutron star surface for a range of plasma parameters and, in §5, we discuss the various implications of our results.

## 2 The Problem

In this paper, we consider the transport of radiation through a strongly magnetized ($`B10^{13}`$ G) atmosphere of a neutron star and construct radiative equilibrium models to calculate the spectrum of the emerging radiation. We assume that the atmosphere is a fully ionized H plasma with an ideal-gas equation of state in local thermodynamic equilibrium (LTE). Owing to the fact that the thickness of the atmosphere ($`10`$ cm) is much smaller than the radius of the neutron star, we treat the atmosphere in plane-parallel geometry. In addition, we assume that the atmosphere has no lateral density variations and that the magnetic field is normal to the neutron star surface. This reduces the radiative transfer problem to one spatial dimension which may be specified, e.g., by the column density in the atmosphere measured from the surface or the equivalent Thomson scattering optical depth $`\tau _\mathrm{T}`$.

The problem, then, consists of solving the equation of radiative transfer with scattering in a magnetized plasma, subject to the constraint of radiative equilibrium, under the assumption of LTE, together with the equation of hydrostatic equilibrium. Below, we first present the processes occuring during the interaction of photons with a magnetized plasma, as well as the approximate ionization equilibrium of the plasma in the presence of a strong magnetic field. We then state the equations for the model-atmosphere problem for the specific case of an ultramagnetized neutron star. Note that, throughout this paper, we express the photon energy $`E`$ and the plasma temperature $`T`$ in keV, and report all quantities as measured by a static observer on the neutron star surface.

### 2.1 Properties of Strongly Magnetized Plasmas

In this section, we briefly summarize some of the properties of electrons and atoms in strong magnetic fields, $`Bm_e^2e^3c/\mathrm{}^3=2.35\times 10^9`$ G, where $`m_e`$ and $`e`$ are the electron mass and charge, c is the speed of light, and $`\mathrm{}`$ is Planck’s constant (for a thorough treatment of this subject, see Mészáros 1992, a recent review by Lai 2001, and references therein). We concentrate, in particular, on field strengths greater than $`B_{\mathrm{cr}}=4.3\times 10^{13}`$ G, at which the cyclotron energy of the electron $`E_\mathrm{b}=eB/m_ec`$ equals its rest mass, 511 keV. At these magnetic field strengths and at temperatures $`T1`$ keV, $`E_\mathrm{b}E_\mathrm{F},E_{\mathrm{th}}`$ where $`E_\mathrm{F}`$ and $`E_{\mathrm{th}}`$ are the Fermi energy and the thermal energy of the electrons, respectively; hence, the electrons are confined to the ground Landau level (in the so-called adiabatic approximation). This greatly limits the motion of the electrons in the transverse direction and the electron is nearly confined to motion in one dimension along the magnetic field. The photon-electron cross section can thus be greatly reduced depending on the direction of polarization of the photon with respect to the magnetic field.

Neutral atoms can also play an important role in plasmas at high magnetic field strengths. The binding energies of atoms are highly enhanced so that, for $`B10^{12}`$ G, the H atom ground-state ionization potential is $`300500`$ eV (Potekhin 1998). Therefore, the ionization equilibrium of a plasma can differ significantly from the nonmagnetic case and the neutral fraction even at a temperature of 1 keV can be large. A competing effect in neutron star atmospheres is that of pressure ionization, which dominates the ionization equilibrium at high densities, $`\rho 10^2\mathrm{g}\mathrm{cm}^3.`$ The most complete treatment to date of ionization equilibrium in strongly magnetized plasmas is by Potekhin, Chabrier, & Shibanov (1999), who used the numerically calculated energy levels of a H atom to compute the thermodynamic properties of a non-ideal gas and the neutral fractions, for a range of temperatures and densities, and for magnetic field strengths of $`10^{13}`$ G.

The present calculations span a range that is higher than the calculations of Potekhin et al. (1999) in all three quantities, i.e., temperature, density, and magnetic field strength. Specifically, most of the photons which carry the flux in the model atmospheres originate at a depth with $`T1`$ keV and $`\rho 10^3\mathrm{g}\mathrm{cm}^3.`$ Extrapolating the binding energies for neutral H at the highest magnetic field strengths calculated by Potekhin (1998), assuming that only the centered (lowest quantum-number) states are populated due to the excluded-volume effects at these high densities, and using the magnetic Saha equation (e.g., Lai 2001), we estimate that the neutral fraction is at most $`f_H<10^5`$. Therefore, in the present work, we make the assumption of complete ionization in the plasma and defer a more complete treatment to a later paper. Note that we also employ in our calculations an ideal gas equation of state which provides a good approximation to the state of the plasma at the temperatures and densities where $`E=`$ few keV photons originate in the atmosphere.

### 2.2 Transport of Photons in Strongly Magnetized Plasmas and Interaction Cross Sections

At magnetic field strengths in the range $`10^{13}\mathrm{G}B10^{15}`$ G which we consider here, the electron cyclotron energy $`E_b`$ is much greater than the plasma temperature $`kT1`$ keV so that the thermal motions of the electrons can be neglected (the cold plasma approximation). At the same time, for these magnetic field strengths, the photon energies in the range of interest ($`0.1E10`$ keV) are $`E_b`$ and therefore lie sufficiently far from cyclotron resonances. Under these conditions, the transport of photons through a plasma can be described by two orthogonal normal modes which correspond to different polarization states of the photons (e.g., Gnedin & Pavlov 1974; Pavlov & Shibanov 1979). Here we adopt $`i=1`$ for the extraordinary mode and $`i=2`$ for the ordinary mode and denote by $`I_E^i`$ the monochromatic intensity of radiation emerging from the plasma in each mode at photon energy $`E`$. Note that in a pure plasma, the normal modes correspond to right- and left-hand circularly polarized photons for propagation angles $`\theta 0`$, while for $`\theta \pi /2`$, they are linearly polarized orthogonal to and along the magnetic field. The normal modes for all other directions of propagation are elliptically polarized. (See §2.3 for the effects of vacuum polarization on the normal modes.)

In our calculations, we take into account scattering as well as emission and absorption due to magnetic bremsstrahlung processes. Because the conditions $`Em_ec^2`$ and $`kTm_ec^2`$ are satisfied, we neglect the change in photon energy in the scattering process (see Gonthier et al. 2000 for these cross sections). The scattering coefficient $`\chi _{ij}(\mu ,\mu ^{})`$ gives the probability of scattering for a photon with direction of propagation $`\mu =\mathrm{cos}\theta `$ with respect to the magnetic field and polarization mode $`i`$ to a direction $`\mu ^{}`$ and polarization mode $`j`$. We follow Kaminker, Pavlov, & Shibanov (1982; eqs. –) to calculate the scattering coefficients. Note that, as in Kaminker et al. (1982), the scattering and absorption coefficients throughout this paper are given in units of $`N_e\sigma _T`$, where $`N_e`$ is the electron density and $`\sigma _T`$ is the angle-averaged Thomson cross section. The total scattering coefficient is denoted by $`\chi _i^s`$ and is given by $`\chi _i^s=_j\chi _{ij}(\mu ,\mu ^{})d\mu ^{}`$. The absorption coefficient $`\chi _i^a`$ is related to the scattering by the generalized Coulomb logarithms and can be found in Kaminker, Pavlov, & Shibanov (1983; eqs. \[A1\]–\[A9\]).

The magnetic scattering and absorption coefficients have strong dependences on the direction of propagation and photon energy, and the two modes show increasingly different behavior with increasing magnetic field strength. To facilitate the discussion in $`\mathrm{\S }4`$ and make our results more intuitive, we plot in Figure 1 the scattering coefficients $`\chi _i^s`$ as a function of $`\theta `$ for $`E/E_b=0.1`$, $`0.01`$, and $`0.001`$. The solid lines show the extraordinary and the dashed lines show the ordinary-mode opacities. In $`\mathrm{\S }4.3`$, we will show that the beaming of emerging radiation for different values of $`E/E_b`$ is closely related to the behavior of the opacities shown in Figure 1.

### 2.3 Vacuum Polarization

In strong magnetic fields, in addition to the response of the plasma electrons, the vacuum polarization effects due to the virtual electron-positron pairs also become significant (Adler 1971; Gnedin et al. 1978; Mészáros & Ventura 1978, 1979) and contribute to the dielectric tensor and the magnetic permeability tensor. Thus, the vacuum polarization alters the polarizations of the normal modes and affects the interaction cross sections in the plasma (see Pavlov & Shibanov 1979; Mészáros & Ventura 1979; Mészáros 1992). The effect of the magnetic vacuum can be quantified by the parameter

$$W=\left(\frac{3\times 10^{28}\mathrm{cm}^3}{N_e}\right)\left(\frac{B}{B_{cr}}\right)^4.$$
(1)
which enters the ellipticity parameter $`q`$ of the normal modes through

$$q=\frac{\mathrm{sin}^2\theta }{2\mathrm{cos}\theta }\sqrt{u}\left(1W\frac{u1}{u^2}\right),$$
(2)
where $`u=E_b^2/E^2`$. The resulting normal modes of the magnetic vacuum, in the absence of plasma electrons, are linearly polarized at $`\theta =0`$ and $`\theta =\pi /2`$, and elliptically polarized at all other angles of propagation. Note that the ellipticity parameter $`q`$ has a slightly modified form for $`B>B_{cr}`$ but this does not affect qualitatively the emerging spectrum (Özel & Narayan 2001).

When both plasma electrons and virtual pairs are present, the normal modes of propagation at very high plasma densities (plasma-dominated region) are different than those at very low plasma densities (vacuum-dominated region) and go through a transition at a critical density where plasma electron and virtual pair contributions identically cancel out (see Fig. 2). At this critical density, a resonance appears at the point where $`q`$ changes sign and the linear polarizations of the normal modes vanish, leaving the modes purely circularly polarized. (Strictly at the critical frequencies, the orthogonality of the normal modes also vanishes, rendering the normal mode treatment invalid; see Mészáros & Ventura 1979 for a discussion). For $`W>4`$, the $`q=0`$ condition can be written as

$$E^2=0.5E_b^2\left[1\pm (14/W)^{1/2}\right],$$
(3)
which gives the critical photon energies for each density where this resonance appears. Alternatively, the same equation can be used to determine the critical density at which the resonance occurs for each photon energy $`EE_b`$.

The most important effect of the vacuum resonance is that, for each photon energy, the opacities of the two modes have a close approach or become equal at the critical density depending on the value of $`\mu `$, with a strong enhancement of both the scattering and absorption coefficients of the extraordinary mode, including the off-diagonal, mode-changing terms of the scattering matrix. Figure 2 shows the opacities with vacuum resonance as a function of $`\tau _T=N_e\sigma _T𝑑z`$ in the atmosphere (see also §3.1). The resonance occurs over a narrow range of densities, which leads to broadband absorption-like features in the spectrum (see §4.3 and also Bulik & Miller 1996). Finally, the eigenmodes of the pure vacuum are similar to those of the electron plasma for most angles $`\theta `$, but differ significantly at a very small range of angles near $`\theta =0`$. Therefore, even though an extraordinary-mode photon can convert into an ordinary-mode photon across the resonance through scattering, the definitions of the “extraordinary” and the “ordinary” modes essentially remain the same above and below the resonant density (see also Fig. 2 and Mészáros 1992).

The vacuum resonance has a finite narrow width in photon energy and electron density. For a given photon energy and direction of propagation, the width of the resonance corresponds to a change in $`q`$, of

$$\mathrm{\Delta }q\frac{\mathrm{sin}^2\theta }{2\mathrm{cos}\theta }\frac{E_b}{E}\frac{\mathrm{\Delta }N}{N},$$
(4)
and occurs over a fractional change in electron density of $`\mathrm{\Delta }N/N1\%`$. We include in our calculations the effects of vacuum polarization through the normal mode ellipticity parameter, as given in Kaminker et al. (1982). However, due to the numerical resolution which prevents us from sampling the change in electron density to the required accuracy even with 400 depth points, we force the resonance to saturate by choosing for $`|q|`$ a minimum allowed value across the resonance of

$$|q|\xi \frac{\mathrm{sin}^2\theta }{2\mathrm{cos}\theta }\frac{E_b}{E}.$$
(5)
We have experimented with the factor $`\xi `$ and found that for our resolution, $`\xi 1\%`$ results in smooth spectra with negligible change in the total and scattering optical depth. This is equivalent to assuming that at a given $`\tau _\mathrm{T}`$, the density fluctuates to within $`1\%`$.

Our results show that specifically, for $`B>10^{14}`$ G and at densities $`N_e10^{27}\mathrm{cm}^3`$, $`W1`$ and therefore the temperature profiles of the atmospheres as well as the emerging spectra at energies above 1 keV are significantly modified by vacuum polarization (see §4.1 and 4.3). Already based on Figure 2, one can see that this resonance is expected to attenuate the emerging flux and lead to a heating of the atmosphere at the $`1\tau _\mathrm{T}10^4`$, which corresponds to the density range that includes the critical density for most flux carrying ($`0.55`$ keV) photons.

### 2.4 Model Atmospheres

We choose as the dimensionless depth variable the Thomson optical depth, $`d\tau _\mathrm{T}=N_e\sigma _Tdz`$, and define a logarithmic grid of typically 100-400 equally spaced points between $`10^5`$ and $`10^5`$. In terms of this stratification parameter, the equation of transfer in the two modes can be written as

$`y_G\mu {\displaystyle \frac{dI_E^i}{d\tau _\mathrm{T}}}`$ $`=`$ $`\chi _i^a\left[I_E^i{\displaystyle \frac{B_E(T)}{2}}\right]+\chi _i^sI_E^i`$ (6)
$`{\displaystyle \underset{j}{}}{\displaystyle \chi _{ji}^s(\mu ,\mu ^{})I_E^j(\mu ^{})𝑑\mu ^{}},`$
where $`B_E(T)`$ is the Planck function, $`\chi _i^a`$ and $`\chi _i^s`$ denote the total absorption and scattering coefficients in Thomson units ($`N_e\sigma _T`$), respectively, and $`\chi _{ij}^s(\mu ,\mu ^{})`$ is the angle- and polarization-dependent scattering coefficient. The parameter $`y_\mathrm{G}`$ takes into account the effects of general relativity and is defined as

$$y_\mathrm{G}\sqrt{1\frac{2GM_{\mathrm{NS}}}{R_{\mathrm{NS}}c^2}},$$
(7)
where $`G`$ is the gravitational constant, $`c`$ is the speed of light, and $`M_{\mathrm{NS}}`$ and $`R_{\mathrm{NS}}`$ are the neutron-star mass and radius, respectively. In the following calculations, we set $`M_{\mathrm{NS}}=1.4M_{}`$ and $`R_{\mathrm{NS}}=10^6`$ cm. Since there are no sources or sinks of heat, the atmospheres are in radiative equilibrium, specified by

$$H_{\mathrm{eff}}=\frac{1}{2}\underset{i}{}I_E^i(\mu )\mu 𝑑\mu 𝑑E\frac{1}{4\pi }\sigma T_{\mathrm{eff}}^4,$$
(8)
where $`4\pi H_{\mathrm{eff}}`$ is the (constant) radiation flux and $`T_{\mathrm{eff}}`$ is the effective temperature.

The momentum balance in the atmosphere is specified by the hydrostatic equilibrium condition

$$\frac{dP}{d\tau _\mathrm{T}}=\frac{g\rho }{y_G^2N_e\sigma _T},$$
(9)
where $`P`$ is the pressure, $`\rho `$ is the density, and $`g=GM_{\mathrm{NS}}/R_{\mathrm{NS}}^2`$. This system of equations describe the state of the plasma and can be closed assuming an equation of state, which to first approximation can be taken to be that of an ideal gas (see §2.1),

$$P=2N_ekT,$$
(10)
where the factor of 2 is valid for completely ionized H gas. Equations (9) and (10) can be easily integrated from the surface inwards, to give

$$N_\mathrm{e}=\frac{gm_\mathrm{p}\tau _\mathrm{T}}{2y_G^2\sigma _TkT},$$
(11)
where $`m_\mathrm{p}`$ is the proton mass. The temperature and density profile of the atmosphere is obtained by solving these equations coupled with the radiation field. The iterative method of solution for this set of equations is described in the next section.

## 3 The Method

We compute the radiative equilibrium model atmospheres using an iteration procedure that involves the solution of the equation of radiative transfer through a Feautrier method and makes use of the linearized moment equations to determine the temperature profile.

### 3.1 Radiative Transfer

We follow a modified Feautrier method (Mihalas 1978) for the solution of the radiative transfer problem. This method allows for a complete, non-iterative solution (of the transfer equation), including the effects of angle-dependent scattering and absorption cross sections. Because we need only to consider conservative scattering, we solve the transfer equation at each photon energy separately.

We begin by writing the equation of radiative transfer for the two modes as

$$y_\mathrm{G}\mu \frac{I_E^i}{\tau _i}=I_E^iS^i,$$
(12)
where

$$d\tau ^i(\chi _i^a+\chi _i^s)d\tau _\mathrm{T}$$
(13)
is the total optical depth at each mode and we have defined the energy-dependent source function as

$$S^i\frac{\chi _i^a}{\chi _i^a+\chi _i^s}\frac{B}{2}+\frac{1}{\chi _i^a+\chi _i^s}\underset{j}{}\chi _{ij}^s(\mu ,\mu ^{})I_E^j(\mu ^{})𝑑\mu ^{}.$$
(14)
In the following equations we drop the energy subscript and explicitly specify any energy-integrated quantities.

We define the Feautrier variables

$$u^i\frac{I_E^{i+}+I_E^i}{2}$$
(15)
and

$$v^i\frac{I_E^{i+}I_E^i}{2}.$$
(16)
We then write the transfer equation  as two first-order equations

$$y_\mathrm{G}\mu \frac{u^i}{\tau _i}=v^i$$
(17)
and

$$y_\mathrm{G}\mu \frac{v^i}{\tau _i}=u^iS^i,$$
(18)
which can be combined to give one second-order equation for the variable $`u`$

$$y_\mathrm{G}^2\mu ^2\frac{^2u^i}{\tau _{i}^{}{}_{}{}^{2}}=u^iS^i.$$
(19)
Note here that, unlike the usual Feautrier equations, in which the source function involves an integral over only photon momentum, the source function in equation  involves also a summation over the photon polarization mode, introducing a coupling between the two equations. However, the off-diagonal (mode-changing) elements of the scattering coefficients, $`\chi _{ij}^s(\mu ,\mu ^{})`$ ($`ij`$), are in general much smaller than the diagonal elements away from the cyclotron resonance, and can be treated as a perturbation. Making use of this fact, we solve iteratively the Feautrier equation for each mode separately, assuming that the coupling terms are known from the previous iteration, until the local specific intensity converges to one part in $`10^9`$.

The transfer problem is highly coupled in the three variables $`\tau _\mathrm{T},E`$, and $`\mu `$ so that choosing a discretization grid is one of the most delicate parts of the calculation. We discretize the second-order equation (19) on the predefined grid over the variable $`\tau _\mathrm{T}`$ (and not over $`\tau _i`$, which depends on $`\mu `$ and $`E`$), according to Mihalas (1978; eq. \[6.26\]-\[6.30\]). We use

$$\mathrm{\Delta }\tau _{i,d\pm 1/2}\frac{1}{2}(\chi _{i,d}^{\mathrm{tot}}+\chi _{i,d\pm 1}^{\mathrm{tot}})|\tau _{\mathrm{T},\mathrm{d}}\tau _{\mathrm{T},\mathrm{d}\pm 1}|$$
(20)
where $`\chi _i^{\mathrm{tot}}\chi _i^a+\chi _i^s`$. The chosen range for $`\tau _\mathrm{T}`$ (see §2.3) corresponds to a range in the optical depth in the extraordinary mode ($`\tau _1`$) of $`0.000110`$ (depending on the photon energy and direction of propagation) and in the ordinary mode ($`\tau _2`$) of $`0.00110^4`$. We chose this range so that all photon energies of interest turn optically thick in the deep interior and decouple from the atmosphere at small optical depths.

Because of the very strong dependence of magnetic opacities on angle near $`\mu 1`$, we choose a logarithmic grid (typically of 30–40 points) in the quantity $`1\mu `$, with $`\mathrm{log}(1\mu )=50`$. We also choose a logarithmic grid (typically of 40 points) in photon energy, with $`\mathrm{log}E=31.5`$. This range spans the photon energies that carry most of the radiative flux for the effective temperatures ($`0.10.5`$ keV) of our atmosphere. Note that the highest number of points in depth, photon energy, and angle are required in those calculations where the vacuum polarization effects are significant.

Finally, we follow the procedure described in Mihalas (1978) to evaluate and invert the Feautrier matrices.

### 3.2 Temperature Correction Scheme

We start with an initial guess for the temperature profile

$$T=\left[\frac{3}{4}\left(\tau _\mathrm{R}+\frac{2}{3}\right)\right]^{1/4}T_{\mathrm{eff}}$$
(21)
where $`\tau _\mathrm{R}`$ is the optical depth that corresponds to the Rosseland mean opacity $`\chi _\mathrm{R}`$, defined as (Pavlov et al. 1992)

$$\frac{1}{\chi _\mathrm{R}}=\frac{3\pi }{8\sigma T^3}_1^1𝑑\mu 𝑑E\left[\frac{1}{\chi _1^{\mathrm{tot}}}+\frac{1}{\chi _2^{\mathrm{tot}}}\right]\frac{B_E(T)}{T}$$
(22)
for two polarization modes and calculate the corresponding density $`N_e(\tau _\mathrm{T})`$ using equation (11). We then proceed by calculating the magnetic opacities with this guess temperature profile on our predetermined grid and solve the radiative transfer problem.

In order to determine the temperature profile of the model neutron star atmosphere in radiative equilibrium, we use an iterative temperature correction scheme similar to the Unsöld-Lucy algorithm, which utilizes the energy-integrated moment equations in order to obtain a linear perturbation equation for $`\mathrm{\Delta }T(\tau _\mathrm{T})`$ (see Mihalas 1978). However, despite the similarity of the approach, the scheme we describe here differs from the Unsöld-Lucy scheme in two important respects:

(i) We do not make any assumptions about the angular distribution of the radiation field, which is the case when the Eddington approximation is used. Instead, we calculate at each iteration the Eddington factor $`f=K/J`$ as a function of $`\tau _\mathrm{T}`$, as well as the ratio $`K/H`$ at the boundary $`\tau _\mathrm{T}=0`$, where $`J,H`$, and $`K`$ denote the zeroth, first, and second moments of the specific intensity, respectively. This takes into account the deviations from anisotropy in the radiation field which are likely to be present in strongly magnetized plasmas.

(ii) We include first-order corrections in all temperature-dependent quantities due to a temperature change $`\mathrm{\Delta }T`$, except for the radiation intensity itself. This is in contrast to the assumption in the Unsöld-Lucy procedure, that only the local blackbody source function changes with temperature.

With these modifications, the correction scheme proceeds as follows. We start with the energy-integrated zeroth and first moment equations

$$y_G\frac{H}{\tau _\mathrm{T}}=\chi _JJ\chi _PB$$
(23)
and

$$y_G\frac{K}{\tau _\mathrm{T}}=\chi _FH,$$
(24)
where $`\chi _J,\chi _P`$, and $`\chi _F`$ are similar to (Mihalas 1978) the absorption mean, Planck mean, and flux mean opacities, respectively, appropriately defined for a magnetized atmosphere. In particular

$$\chi _P\frac{\pi }{2}\left[\frac{𝑑\mu 𝑑EB_E(T)(\chi _1^a+\chi _2^a)/2}{\sigma _TT_{\mathrm{eff}}^4}\right],$$
(25)
$$\chi _J\frac{𝑑\mu 𝑑E(\chi _1^au^1+\chi _2^au^2)}{𝑑\mu 𝑑E(u^1+u^2)},$$
(26)
and

$$\chi _F\frac{𝑑\mu 𝑑E\mu (\chi _1^{\mathrm{tot}}v^1+\chi _2^{\mathrm{tot}}v^2)}{𝑑\mu 𝑑E\mu (v^1+v^2)}.$$
(27)
In the above definitions, we implicitly assume that the moments and the mean opacities depend on $`\tau _\mathrm{T}`$ and omit any explicit references.

Integrating equation (24) and substituting into equation (23) by making use of the Eddington factors, which are known from the current iteration, yields

$$\chi _PB=\chi _J\left(\frac{J}{K}\right)\left[\left(\frac{K}{H}\right)_0H(0)+\frac{1}{y_G}\chi _F(\tau _\mathrm{T}^{})H(\tau _\mathrm{T}^{})𝑑\tau _\mathrm{T}^{}\right]y_G\left(\frac{dH}{d\tau _\mathrm{T}}\right),$$
(28)
where the subscript 0 corresponds to the outer boundary. We assume that the constant flux solution can be obtained by making a small correction $`\mathrm{\Delta }T(\tau _\mathrm{T})`$ everywhere in the atmosphere and allow for the local black body and the mean opacities to vary linearly with $`\mathrm{\Delta }T`$, which gives $`\mathrm{\Delta }B(\tau _\mathrm{T})=(B/T)\mathrm{\Delta }TB^{}(\tau _\mathrm{T})\mathrm{\Delta }T,\mathrm{\Delta }\chi _P=\chi _P^{}\mathrm{\Delta }T`$, and similarly for the other mean opacities. Note that all the derivatives are computed at the $`\tau _\mathrm{T}`$ grid points. Now proceeding as in the Unsöld-Lucy scheme, we note that in the desired radiative equilibrium solution, the term $`H/\tau _\mathrm{T}`$ vanishes and the radiation flux becomes $`H(\tau _\mathrm{T})=H_{\mathrm{eff}}=\frac{1}{4\pi }\sigma _TT_{\mathrm{eff}}^4`$. Thus, subtracting the current equation from the desired constant flux solution and solving for $`\mathrm{\Delta }T`$ we obtain

$$\mathrm{\Delta }T(\tau _\mathrm{T})=\frac{\chi _J(J/K)_\tau \left[(K/H)_0\mathrm{\Delta }H(0)+I_1\right]+y_G(\mathrm{d}H/\mathrm{d}\tau )_\tau }{\chi _P^{}B+\chi _PB^{}\chi _J^{}(J/K)_\tau (K/H)_0H_{\mathrm{eff}}+I_2]},$$
(29)
where

$$I_1=\frac{1}{y_G}_0^\tau [\chi _F(\tau _\mathrm{T}^{})\mathrm{\Delta }H(\tau _\mathrm{T}^{})+\chi _F^{}(\tau _\mathrm{T}^{})\mathrm{\Delta }T(\tau _\mathrm{T}^{})H_{\mathrm{eff}}]d\tau _\mathrm{T}^{},$$
(30)
$$I_2=\frac{1}{y_G}_0^\tau [\chi _F(\tau _\mathrm{T}^{})+\chi _F^{}(\tau _\mathrm{T}^{})\mathrm{\Delta }T(\tau _\mathrm{T}^{})]H_{\mathrm{eff}}d\tau _\mathrm{T}^{},$$
(31)
and $`\mathrm{\Delta }H(\tau _\mathrm{T})=H_{\mathrm{eff}}H(\tau _\mathrm{T})`$ refers to the deviation of the current solution from the constant flux solution at any depth in the atmosphere. To quantify this deviation, we define the error $`e=[H_{\mathrm{eff}}H(\tau _\mathrm{T})]/H_{\mathrm{eff}}`$ and quote the convergence of the obtained solutions in terms of this quantity.

Due to the assumption of linear changes in various quantities in response to a temperature correction $`\mathrm{\Delta }T`$, care must be taken when applying large temperature corrections which may violate this assumption. To this end, we have introduced a damping factor to control the actual applied change everywhere in the atmosphere. We experimented with damping factors that depend on the magnitude of $`\mathrm{\Delta }T(\tau _\mathrm{T})`$ and found that reducing the $`\mathrm{\Delta }T`$ obtained from the above scheme typically by $`0.30.5`$, even though somewhat slowing down the convergence rate, results in a marked increase in the stability of the algorithm. With this method, we achieve flux convergence down to a maximum error of $`e5\%`$ in $`O(10)`$ iterations, improving it to everywhere better than $`1\%`$ in $`O(100)`$ iterations.

Comparing the scheme described here with the Unsöld-Lucy algorithm, we note that the correction in temperature at any point in the atmosphere, $`\mathrm{\Delta }T(\tau _\mathrm{T})`$, depends not only on the errors in the flux everywhere in the atmosphere outward of that location (as in the Unsöld-Lucy scheme) but also on the temperature corrections applied at all other points. This allows for stronger coupling between all points in the atmosphere and leads to a more uniform convergence. At the same time, convergence is much more stable, and artificial temperature inversions that commonly arise when applying the Unsöld-Lucy scheme to magnetized atmospheres are not present.

Finally, note that any temperature correction scheme that allows for a change in the intensity field with temperature, such as the complete linearization schemes (Mihalas 1978), couples the different photon energies in the solution of the radiative transfer problem and thus requires $`O(N_E^3)`$ operations in the Feautrier method, where $`N_E`$ is the number of grid points in photon energy. As a result, the total computational time required for solving the radiative transfer problem scales as $`O(M^3N_E^3)`$, with M denoting the number of points on the $`\mu =\mathrm{cos}(\theta )`$ grid, thus limiting the computational resolution in one or both of these variables. In such a highly angle- and photon energy-dependent problem as the magnetized radiative transfer, we have opted to work with the largest number of grid points in a feasible amount of computational time and therefore chose to work with the temperature correction scheme discussed above.

We apply the method described above to construct magnetized, radiative equilibrium atmospheres with a flux $`H(\tau )`$ that is converged to within $`1\%`$ of $`H_{\mathrm{eff}}`$, everywhere in the atmosphere. Figure 3 shows the magnitude of typical errors that result from such a calculation. Throughout most of the atmosphere, the fractional error is $`e1\%`$, reaching $`10^8`$ in the deep interior. We have found that the slowest convergence, and thus the maximum error, occurs typically in the transition region, where the photons carrying the largest flux become optically thin; we will discuss this transition region below. We also note that the convergence criterion of $`1\%`$ does not arise from an intrinsic limitation of the algorithms, but is chosen in view of the required computational time.

## 4 Results

We have applied the method described above to construct magnetized, radiative equilibrium atmospheres of neutron stars with a range of magnetic field strengths ($`10^{13}`$ G$`B10^{15}`$ G) and effective temperatures ($`0.1`$ keV$`<T_{\mathrm{eff}}<0.5`$ keV). We show below the resulting temperature profiles, spectra, and beaming of emerging radiation and identify the effects of magnetic field strength, effective temperature, scattering, and vacuum polarization on each of these properties.

### 4.1 Temperature Profiles

Figure 4 shows the temperature profiles for three model atmospheres with different magnetic field strengths, at $`T_{\mathrm{eff}}=0.3`$ keV (top panel). These profiles have distinct features arising from the strong angle- and energy-dependence of the magnetic opacities and the presence of two polarization modes and are significantly different from the temperature profile of a grey atmosphere (see, e.g., Mihalas 1978) and from an atmosphere calculated using Rosseland mean opacities (dashed lines). In particular, at low optical depths, the temperature is determined primarily by the ordinary-mode opacities, which depend very weakly on the magnetic field strength, for most directions of propagation. Therefore, the temperature profiles in this region are similar in all three model atmospheres (but see below for the role of vacuum polarization in this region). This is also the reason why the Rosseland mean opacities, which are dominated by the small opacity (extraordinary) mode, fail to describe the temperature of the outer atmosphere correctly. At the limit of high optical depth, the diffusion limit is reached and the temperatures calculated using the Rosseland mean opacities provide an accurate description of the atmospheres. At these high optical depths, it is the extraordinary-mode opacities, which are low and depend strongly on the field strength, that determine exclusively the temperature profiles. As a result, with increasing magnetic-field strength, the same plasma temperature is reached deeper in the atmosphere.

The presence of multiple plateaus in the temperature profiles correspond to optically thick/thin transitions of photons with different energies, directions of propagation, and polarizations. Specifically, the plateaus at $`\tau _\mathrm{T}10^010^2`$ arise from such transitions of extraordinary-mode photons that carry most of the flux, while the plateaus at $`\tau _\mathrm{T}10^4`$ are the result of such transitions for the ordinary-mode photons with low energies or large angles of propagation.

One other significant difference between the profiles of Figure 4 and the temperature profile of a grey/Rosseland atmosphere is the relation between the effective temperature and the plasma temperature $`T_0`$ at the outer boundary. For a grey atmosphere, $`T_00.8T_{\mathrm{eff}}`$, whereas for the calculations shown in Figure 4, $`T_00.1T_{\mathrm{eff}}`$. The same conclusion is reached by Shibanov et al. (1992), but the agreement is only qualitative due to the insufficiency of the diffusion approximation employed in that work in capturing the details of the temperature profiles.

The temperature profiles of model atmospheres with different effective temperatures, as shown in Figure 4 (bottom panel), have the same qualitative characteristics, but show an overall displacement. In addition, the plateaus occur further out in the atmosphere with decreasing effective temperature, as the absorption cross sections increase with decreasing temperature, moving the optically thick/thin transitions to smaller column densities.

At $`B10^{15}`$ G, the temperature profiles are significantly affected by the vacuum polarization resonance. To demonstrate the vacuum resonance effects on the temperature profile, we show in Figure 5 the results of two calculations at $`B=10^{15}`$ G and $`T_{\mathrm{eff}}=0.3`$ keV, where we include (solid lines) or neglect (dashed-lines) the vacuum polarization contribution to the magnetic opacities. The temperature is higher for a range of optical depths $`1\tau _\mathrm{T}10^3`$ in the case where the vacuum polarization contribution to the opacities are included. Such an increase in temperature is typical of radiative equilibrium atmospheres with a resonant layer at low densities for reasons discussed below.

The net effect of the vacuum contribution to the magnetic opacities is to increase the extraordinary-mode opacity for a narrow, photon energy-dependent range of densities in the atmosphere (see §2.2 and Fig. 2). For $`B10^{15}`$ G, this resonance appears at $`1\tau _\mathrm{T}10^3`$. Because this range lies within the photosphere where the extraordinary mode photons of $`0.510`$ keV have decoupled from the atmosphere, the resonance leads to an attenuation of the flux at these photon energies, which appears as a broad-band absorption feature in the spectrum (see §4.3). At the same time, the additional absorption of radiation heats the atmosphere, thus leading to an increase in the temperature at this range of column densities (Fig. 5). This in turn increases the emission at these column depths and contributes to the flux at lower photon energies. We note that it is also possible to understand these results in terms of the radiative equilibrium condition in the atmosphere. In the presence of resonant “layers” of high opacity, the overall temperature in the atmosphere has to increase in order for the energy-integrated flux to be constant throughout the atmosphere, for a given effective temperature. This change in the temperature profile results in a corresponding change in the entire emerging spectrum at the highest magnetic field strengths, as we will show in the next section.

### 4.2 Spectra of Emerging Radiation

In Figure 6, we show the spectra emerging from neutron star atmospheres with different magnetic field strengths and effective temperatures. In the top panel, the flux is decomposed into the relative contributions of the two polarization modes while in the bottom panel the total flux is shown. Note that the fluxes shown are in accordance with the definition of $`H`$ given in $`\mathrm{\S }2.3`$ and the total radiative flux is given by $`4\pi H`$; the subscript $`E`$ as usual denotes the energy dependence. The detailed features of the spectra result mainly from the combined effects of the temperature profile, the scattering between polarization modes, and the polarization of the magnetic vacuum, which we discuss below. The magnetic field strength enters all of these factors and thus affects the emerging spectra in multiple ways.

Two results hold in all cases: (i) the extraordinary mode carries most of the radiative flux and (ii) the spectra are not Planckian, showing various amounts and shapes of hard excess. The first result arises from the fact that the opacity of the extraordinary mode is significantly smaller and hence it decouples from the plasma at greater depths, where the temperature is higher (Fig. 6). The second result is primarily, but not exclusively, due to a similar reason: as shown previously (e.g., Shibanov et al. 1992; Rajagopal & Romani 1997), the high-energy photons originate deeper in the atmosphere due to their smaller interaction cross section and hence the emerging spectra are bluer than the blackbody spectrum at the effective temperature. The spectral hardening because of the energy-dependence of the opacities is less pronounced in the case of high magnetic fields than non-magnetic neutron star atmospheres because the dependence of the magnetic free-free opacity on photon energy $`(KE^1)`$ is weaker than that of the non-magnetic one $`(KE^3)`$ (see also Zane, Turolla, & Treves 2000). However, it is important to note that in the case of ultrastrong magnetic fields, different photospheres for different energy photons is not the only origin of hard excess in the spectra and the hard tails do not decrease monotonically with increasing magnetic field strength. Instead, the spectral shape is determined by two other effects. First, as we will show in §4.3, the broad-band absorption due to vacuum resonance changes the spectral shape to reduce the color temperature $`T_c`$ of the best-fit blackbody but gives rise to a pronounced power-law tail. This is because of the energy- and density-dependence of the vacuum resonance as shown in Figure 2. Second, scattering of photons from the extraordinary mode to the ordinary mode, especially enhanced in the presence of vacuum polarization resonance, reduces the spectral hardening. We also show this below.

As the bottom panel in Figure 6 shows, the effective temperature controls the total emerging flux, as well as the peak and broadness of the spectrum. Although the scaling of the total flux with $`T_{\mathrm{eff}}`$ is straightforward, the flux at each photon energy does not scale in the same way. This is because the temperature profiles corresponding to different $`T_{\mathrm{eff}}`$ differ by more than a simple displacement (see §4.1). Specifically, the outer layers can be significantly hotter for higher effective temperatures, giving rise to more flux at low photon energies. For a similar reason, the spectra emerging from atmospheres with higher $`T_{\mathrm{eff}}`$ are more peaked, with an increasing fraction of the total flux emerging around the peak photon energy.

The effects of scattering between different polarization modes are strongest at the lowest magnetic-field strengths, because, for a given photon energy, the off-diagonal element of the scattering coefficient from the extraordinary mode to the ordinary mode increases with decreasing field strength. On the other hand, the effects of the polarization of the magnetic vacuum are strongest at the highest field strengths, as discussed in §2.2. In order to disentangle these two effects, we discuss the lowest ($`10^{13}`$ G) and the highest ($`10^{15}`$ G) magnetic field strengths in more detail, thus isolating each effect.

We use the temperature profile of the model atmosphere that corresponds to $`B=10^{13}`$ G and $`T_{\mathrm{eff}}=0.3`$ keV and calculate the emerging spectra, including or neglecting the scattering between polarization modes. As can be seen in Figure 7, the scattering between modes increases the flux of the ordinary mode at high energies. This occurs because the off-diagonal elements of the scattering coefficients increase, when $`E`$ approaches $`E_\mathrm{b}`$. An observable result of this cross-mode scattering is a small suppression of the hard excess at $`B10^{13}`$ G. The discussion of the vacuum polarization effects which are strong at $`B10^{15}`$ G are presented in the next section.

### 4.3 Effects of Vacuum Polarization on the Spectrum

Taking into account the contribution of vacuum polarization to the magnetic opacities leads to significant modification in the spectrum of radiation emerging from the neutron star atmosphere. To demonstrate the effects of vacuum resonance, we show in Figure 8 two model spectra where include (solid line) or neglect (dashed line) the vacuum contribution to the scattering and absorption cross sections, using model atmospheres that corresponds to $`B=10^{15}`$ G and $`T_{\mathrm{eff}}=0.3`$ keV. The result is a drastic change in the overall spectral shape: an increase of the flux at $`E2`$ keV, an attenuation of the flux at $`E2`$ keV and the emergence of a hard spectral tail at $`E3`$ keV due to broadband absorption across the vacuum resonance.

As Figure 2 shows, the main effect of the contribution of the vacuum polarization to the magnetic opacities for most angles of photon propagation is the emergence of a resonance at a photon-energy dependent critical density in the plasma. When the effects of vacuum polarization are included, at each photon energy, the opacity of the extraordinary mode becomes comparable to that of the ordinary mode and the scattering from the extraordinary mode into the ordinary mode is enhanced for a narrow range of densities. This rapid increase in opacity of the extraordinary mode photons typically occurs at low densities further out in the atmosphere (see Fig. 2), and thus at a density where these photons have already decoupled from the atmosphere. Therefore, the reduction in the flux arises from the fact that going outward in the atmosphere, the essentially free-streaming extraordinary mode photons encounter a resonance and either change polarization mode or are absorbed over a narrow range of $`\tau _T`$. In the latter case, they decouple again outside the resonant layer and stream outward. Note that in the absence of temperature adjustments in the atmosphere, flux in the extraordinary mode at nearly all photon energies would suffer such a reduction. On the other hand, at the densities where the vacuum effects become important, the ordinary mode is mostly optically thick and the reduction of the ordinary-mode opacity due to vacuum contribution is minimal. As a result, the primary effect of vacuum polarization is an attenuation of the flux carried by the extraordinary mode at high photon energies.

The radiative equilibrium condition introduces an additional effect and a further modification of the emerging spectrum. This is because the atmosphere adjusts to the increase in opacity, and the corresponding reduction in the flux, by an overall increase in the temperature in the resonant layers (see §4.1 and Fig. 5). This can also be understood as a heating of the atmosphere due to resonant absorption. This in turn leads to an increase of the flux carried by the low-energy photons (so that the energy-integrated flux is constant) and therefore to an overall shift of the spectrum towards lower energies and a further distortion away from a Planckian spectrum.

### 4.4 The Beaming of Emerging Radiation

Finally, we discuss the angular pattern of the radiation emerging from the different model atmospheres. In Figure 9 (top panel), we show the beaming of radiation at a photon energy of 1 keV and for three different magnetic-field strengths, while in the bottom panel, we show the beaming of radiation emerging from a model atmosphere with $`B=10^{14}`$ G, at three different photon energies. Note that the emerging intensity in the ordinary and extraordinary modes are added together in these figures.

Both figures show that the radiation pattern consists of two components: a narrow pencil beam at small angles ($`5^{}`$) and a broader fan beam at intermediate angles ($`2060^{}`$). This pattern is a direct consequence of the angular dependence of the cross sections in the two modes (see Fig. 1), combined with the temperature profiles in the atmospheres (Fig. 4). In particular, the extraordinary mode contributes to both the pencil and the fan beams, because its opacity is peaked at an angle away from $`\theta =0`$, decreasing towards both smaller and larger angles. On the other hand, the ordinary mode contributes (roughly equally as the extraordinary mode) to the narrow pencil beam of the pattern, because its opacity decreases significantly at small angles. Therefore while not contributing significantly to the total flux (see §4.2), the ordinary mode plays an important role in determining the beaming pattern with its large specific intensity near $`\theta =0`$.

The angular dependence of the magnetic cross sections, and hence the beaming of radiation, depend both on photon energy and on magnetic field strength, through the ratio $`E/E_\mathrm{b}`$. As Figure 9 shows, at photon energies much smaller than the electron cyclotron energy, the pencil beam is constrained to increasingly smaller angles, $`\mathrm{sin}\theta (E/E_\mathrm{b})^{1/2}`$ while the fan beam becomes more prominent, resulting in a flatter radiation pattern. Therefore, for a given magnetic field, the emerging radiation is more strongly and radially beamed at high photon energies, and for a given photon energy the radiation is more beamed at lower magnetic field strengths. Note that with increasing field strength or decreasing photon energy, the peak value of the pencil beam may in fact increase; however the total flux in the pencil beam decreases, thus making it less prominent. Numerically, the pencil beam in such cases, constrained to $`\theta 0.1^{}`$, becomes hard to resolve and has minimal observational consequences.

In addition, at $`B10^{15}`$ G, vacuum polarization affects the beam shape at photon energies $`E1`$ keV. Specifically, the scattering cross sections are enhanced due to the vacuum contribution, redistributing the photons between the two polarization modes and different angles of propagation more efficiently. This results in a beam shape with a broader pencil component and a broader plateau between the pencil and fan beams than that of the $`B=10^{14}`$ G case (Fig. 8) and the case when the vacuum contribution is ignored.

## 5 Discussion and Conclusions

We have constructed model atmospheres for strongly magnetic neutron stars in radiative equilibrium. We have considered the complete angle, polarization and energy dependence of the magnetic free-free absorption and scattering processes, taking into account the effects of polarizability of the magnetic vacuum. We find that the resulting spectra are broader and bluer than the black body spectrum at the effective temperature, and depend on the magnetic field strength.

Figure 10 shows the dependence of the emerging spectra on magnetic field strength for $`T_{\mathrm{eff}}=0.3`$ keV along with the best-fit black body spectra in the photon energy range $`0.510`$ keV, displaced vertically for clarity. In determining the best fit blackbody and the color temperature $`T_\mathrm{c}`$, we fit the peak of the spectra and allow for a hard excess, either in the form of a power-law tail or a broad excess. For $`B10^{14}`$ G, the spectra have distorted Planckian shapes and are generally broader than a black body. Therefore, when fitted with black body shapes, the best-fit (color) temperatures, $`T_\mathrm{c}`$ depend on the range of photon energies used (see Table 1). Typically, for these magnetic fields, $`T_\mathrm{c}/T_{\mathrm{eff}}1.51.8`$.

At $`B>10^{14}`$ G, the spectra are significantly modified by the vacuum polarization effects and are in general well-described by a black-body shape at low energies with a high energy power-law tail. The corresponding color temperatures are smaller (see Table 1), with $`T_\mathrm{c}/T_{\mathrm{eff}}1.2`$ at $`B=10^{15}`$ G. For this magnetic field strength, the power-law tail has an energy index of $`4`$.

At all magnetic field strengths, the angular dependence has two maxima: a narrow (pencil) beam at small angles ($`5^{}`$) with respect to the normal and a broad maximum (fan beam) at intermediate angles ($`2060^{}`$). The relative importance and the opening angle of the radial beam decreases strongly with increasing magnetic field strength and decreasing photon energy (see §4.3). As a result, the emerging radiation from a thermally-emitting neutron star has a prominent radial peak only at relatively low magnetic field strengths ($`B10^{13}`$ G) and at high photon energies ($`E1`$ keV) within the range of $`B`$ and $`E`$ values considered here.

Our results have significant implications for models of thermally emitting neutron stars with strong magnetic fields, and specifically for magnetar models of SGRs and AXPs. First, the observed power-law spectral tails at high photon energies need not arise from a non-thermal (possibly magnetospheric) emission mechanism but may be the result of vacuum polarization at very strong magnetic fields (see Fig. 9). Second, the beaming of the emerging radiation determines the pulse shapes and amplitudes of spinning neutron stars with anisotropic surface temperature distributions. These properties depend also on the bending of photon trajectories in the gravitational field of neutron stars and can be used in constraining models of SGRs and AXPs (DeDeo, Psaltis, & Narayan 2001; Psaltis, Özel, & DeDeo 2000). We will report the results of a detailed comparison of the observed spectral and timing properties of SGRs and AXPs with the predictions of the models discussed here in a forthcoming paper (Özel, Psaltis, & Kaspi 2001).

I am grateful to Ramesh Narayan for countless stimulating discussions and especially for his help in developing the temperature correction scheme. I am also grateful to Dimitrios Psaltis for his help in many aspects of constructing and understanding models of neutron star atmospheres throughout this work. I thank George Rybicki and Dimitar Sasselov for useful discussions on radiative transfer. I have benefited a lot from discussions with George Pavlov, Dong Lai, and Roger Romani on magnetic processes during the conference on “Spin and Magnetism in Young Neutron Stars” at ITP. I also thank Vicky Kaspi, Cole Miller, George Pavlov, and Peter Woods for useful comments on the manuscript. This work was supported in part by NSF Grant AST 9820686.


format: text

text:
# On a possible manifestation of 𝑓₁ trajectory in 𝐽/𝜓 photoproduction

## Abstract

We analyze a possible manifestation of $`f_1`$-trajectory in elastic $`J/\psi `$ photoproduction at high energy and large momentum transfer. Inspite of the small contribution of $`f_1`$-trajectory in total cross sections, it becomes significant in various spin observables. In particular, we show that the crucial test for $`f_1`$-exchange can be made by measuring the single beam- and double parity- and beam-target asymmetries at large momentum transfers, where a strong deviation from the exchange of conventional Pomerons is expected. This effect is caused by the interference of natural (Pomeron) and unnatural ($`f_1`$) parity exchange parts of amplitude in the region where their contributions become comparable to each other and might be interesting to observe in forthcoming experiments, if feasible.

preprint: KUHD/00-8
Hadronic diffraction has regained popularity in the recent years due to many new interesting experiments at HERA and TEVATRON. One of them is the vector meson photoproduction at HERA. This channel is exciting because here one can study any of the known vector mesons from the lightest $`\rho `$ to heavier $`J/\psi `$ and $`\mathrm{{\rm Y}}`$ and hope to see the transition from the “soft” to the “hard” regime. The most popular models for description of the diffractive processes at high energy are based on the Regge theory  where the corresponding multiparticle exchanges with definite quantum numbers are expressed by the effective exchange of one Regge pole whose propagator is given by $`(s/s_0)^{\alpha (t)}`$. The Regge trajectory $`\alpha (t)`$ has a simple linear form $`\alpha (t)=\alpha (0)+\alpha ^{}t`$ and the main contribution to the elastic forward photoproduction of light vector mesons comes from the Pomeron (so-called “soft” Pomeron) trajectory with intercept 1.08 and slope of 0.25 GeV-2 . This leads to rather weak energy dependence of the cross section gradually becoming almost constant at very high energies. However, the Pomeron, though amply used to successfully analyze the high energy total hadron cross sections, still remains mysterious in its partonic content, hence its manifestation in QCD also remains unknown.

Recent studies of elastic vector meson photoproduction data shows that the $`J/\psi `$ cross-section shows a steep rise with the energy $`W`$ . This behavior is quite different from that of other light vector mesons like $`\rho ,\omega `$ and $`\varphi `$ which are characterized by a weak dependence of their cross section on energy \[$`\sigma W^{0.22}`$\], while for $`J/\psi `$ it is parameterized as $`W^{0.8}`$. This new phenomenon drew a lot of attention in this field leading to several models being proposed to explain the data . However, this subject still remains controversial. Donnachie and Landshoff  claimed that perturbative QCD does not work in this energy region and suggested the use of a small admixture of the second “hard” Reggeon exchange whose trajectory is given by $`\alpha _h(t)=1.418+0.1t`$. The conventional mesonic trajectories arising from the fits of Donnachie and Landshoff (for instance, $`f_2`$-trajectory ) contribute to light vector meson productions at relatively low energies. Most of them have a slope of 0.9 GeV-2 and an intercept $`\alpha (0)0.5`$. These trajectories with positive and negative signature are important in the near-threshold region though their contribution decreases with energy. But for description of the $`J/\psi `$ photoproduction for all energy regions, one must take into account all of the above trajectories: Pomeron, hard Reggeon and mesonic (supersoft) trajectories .

The Pomeron and “hard” Reggeon which are dominant at high energy on the phenomenological level are described by the effective “$`C=+1`$ photon” exchange which leads to similar spin observables for all components typical for the amplitude with $`t`$ \- channel natural parity exchange. However, so far the study of the unnatural parity exchange trajectories at high energy has not been done widely. Recently, an interesting analysis for testing the new $`f_1`$-trajectory associated with the axial vector meson exchange was done by Kochelev et al. , in which the slope of the trajectory was taken to be almost zero. It is very interesting to test the effect of this trajectory even for heavy vector meson $`J/\psi `$ photoproduction. Since its slope is close to zero, it may be active in the high energy and large momentum transfer region.

In this work, we are concerned with the large momentum transfer region, where the conventional soft Pomeron, the hard Reggeon and $`f_1`$ trajectory contribute. We show that the spin observables are sensitive to $`f_1`$-trajectory and will help to differentiate the different components of the amplitude in this region. Let us start with the definition of the kinematical variables for the $`\gamma pJ/\psi p`$ reaction using standard notation. The four momenta of the incoming photon, outgoing $`J/\psi `$, initial (target) and final (recoil) proton are represented by $`k,q,p`$ and $`p^{}`$, respectively. Hereafter, $`\theta `$ denotes the $`J/\psi `$ production angle in c.m.s. and $`sW^2=(p+k)^2,M_N`$ is the nucleon mass, $`M_V`$ is the $`J/\psi `$ mass and $`m_{f_1}`$ is the $`f_1`$-meson mass. We use the convention of Bjorken and Drell to define the $`\gamma `$ matrices and Dirac spinors.

For simplicity we call the “soft” Pomeron and the “hard” Reggeon trajectories \[Donnachie-Landshoff(DL)\] as the Pomeron part . The corresponding invariant amplitude in the standard notation reads,

$`T_{fi}^P=i{\displaystyle \underset{i=s,h}{}}\overline{u}_{m_f}(p^{})M_{R_i}ϵ_{V\mu }^{}(\lambda _f)((\gamma k)g^{\mu \nu }\gamma ^\nu k^\mu )ϵ_{\gamma \nu }(\lambda _i)u_{m_i}(p),`$ (1)
The factor $`M_{R_i}`$ is given by the conventional Regge pole amplitude,

$`M_{R_i}=12a_is_0\beta _0^2\left({\displaystyle \frac{3\mathrm{\Gamma }_{Ve^+e^{}}}{\alpha _{\mathrm{em}}M_V}}\right)^{1/2}F_N(t)F_V(t)e^{i\frac{\pi }{2}(\alpha _i(t)1)}\left({\displaystyle \frac{s}{s_0}}\right)^{\alpha _i(t)},`$ (2)
where $`i=s,h`$ for soft Pomeron and hard Reggeon, respectively. $`VJ/\psi `$; $`s_0=4`$ GeV2, $`\beta _0=1.64`$ GeV-1; $`a_s=1`$, $`a_h=0.05`$; $`\alpha _{\mathrm{em}}=1/137`$, $`\mathrm{\Gamma }_{J/\psi e^+e^{}}`$ = 5.26 KeV . The trajectories of these Pomerons are given by,

$`\alpha _s(t)=1.08+0.25t,`$ (3)
$`\alpha _h(t)=1.418+0.1t.`$ (4)
$`F_N`$ is the isoscalar electromagnetic form factor of the nucleon and $`F_V`$ is the form factor for the vector-meson-photon-Pomeron coupling

$`F_N(t)`$ $`=`$ $`{\displaystyle \frac{(4M_N^22.8t)}{(4M_N^2t)(1t/0.7)^2}},`$ (5)
$`F_V(t)`$ $`=`$ $`{\displaystyle \frac{M_V^2}{(M_V^2t)^2}}{\displaystyle \frac{\mu ^2}{2\mu ^2+M_V^2t}},`$ (6)
with $`\mu =1.1`$ GeV2.

For the $`f_1`$-trajectory, we follow the consideration of Ref.  where it is assumed that the intercept of this trajectory is equal to 1, in agreement with the low $`x`$ behavior of the spin-dependent deep inelastic structure function $`g_1(x,Q^2)1/x^\alpha `$ with $`\alpha =0.9\pm 0.2`$. Moreover it is assumed that the slope of the trajectory is close to zero, which is in agreement with large $`t`$ behavior of elastic $`pp`$ scattering. This means that the excited meson states of this trajectory have very large masses and the total amplitude with unnatural parity is saturated by the lowest $`f_1(1285)`$ meson exchange which may be expressed in terms of the one boson exchange amplitude. The effective vertices of the axial-vector $`f_1`$-meson interaction with the nucleon and the $`f_1`$-photon-vector-meson interaction reads,

$`V_{f_1NN}`$ $`=`$ $`ig_{f_1NN}\overline{u}_f\gamma _\mu \gamma ^5\overline{u}_i\xi ^\mu ,`$ (7)
$`V_{f_1V\gamma }`$ $`=`$ $`g_{f_1V\gamma }ϵ_{\mu \nu \alpha \beta }\xi ^\beta ϵ_1^\nu ϵ_2^\alpha M_V^2k^\mu ,`$ (8)
where $`\xi `$ and $`ϵ`$ are the polarization vectors of the axial vector and vector mesons, respectively, and $`g_{f_1NN}=2.5`$ is fixed from the proton spin analysis. The coupling constant $`g_{f_1V\gamma }`$ is extracted from the radiative decay $`J/\psi f_1\gamma `$,

$$\mathrm{\Gamma }_{Vf_1\gamma }=g_{f_1V\gamma }^2\frac{(M_V^2M_{f_1}^2)^3(M_V^2+M_{f_1}^2)}{96\pi M_VM_{f_1}^2}.$$
(9)
Taking $`\mathrm{\Gamma }_{J/\psi f_1\gamma }=56.5`$ eV , one gets $`|g_{f_1J/\psi \gamma }|=1.245\times 10^4`$ GeV-2.
Using the above vertex functions and coupling constants, the corresponding matrix element of the vector meson production can be written as,

$`T_{fi}=ig_{f_1V\gamma }g_{f_1NN}F_{f_1NN}(t)F_{f_1V\gamma }(t){\displaystyle \frac{m_V^2}{tm_{f_1}^2}}ϵ_{\mu \nu \alpha \beta }q^\mu ϵ_V^\nu (\lambda _V)ϵ_\gamma ^\alpha (\lambda _\gamma )`$ (10)
$`\times \left(g^{\beta \delta }{\displaystyle \frac{(pp^{})^\beta (pp^{})^\delta }{m_{f_1}^2}}\right)\overline{u}_{m_f}(p^{})\gamma _\delta \gamma _5u_{m_i}(p),`$ (11)
where $`m_{f_1}=1.285`$ GeV. For the $`f_1NN`$ vertex, we use the flavor singlet axial vector form factor which is fixed from fitting to elastic $`pp`$ scattering at high $`\sqrt{s}`$ and large $`|t|`$ . (See also .)

$`F_{f_1NN}(t)=1/(1t/m_{f_1}^2)^2.`$ (12)
The form factor of the $`f_1V\gamma `$ vertex is taken as

$$F_{f_1V\gamma }(t)=\frac{\mathrm{\Lambda }_V^2m_{f_1}^2}{\mathrm{\Lambda }_V^2t}$$
(13)
where $`\mathrm{\Lambda }_V`$ for $`J/\psi `$ is estimated to be $`\mathrm{\Lambda }_{J/\psi }=4.2`$ GeV by extrapolating the $`\mathrm{\Lambda }_V`$ found for $`\rho `$ and $`\varphi `$ mesons to the one for the heavy $`J/\psi `$ meson assuming a linear dependence.

Note that the $`t`$-dependence of Eq.(11) is flatter than that of Eq.(1). One reason for this is that the conventional Pomerons have finite trajectory slopes ($`\alpha _i>0`$) whereas the slope of $`f_1`$ trajectory is almost zero. A slight difference in $`t`$-dependence of the corresponding form factors $`F_NF_V`$ and $`F_{f_1NN}F_{f_1V\gamma }`$ also play a part in determining the overall $`t`$-dependence.

Fig. 1 shows the differential cross sections of $`J/\psi `$ photoproduction for $`W=20`$ and $`94`$ GeV. One can see that although the $`f_1`$-contribution is not conspicuously seen at forward production angles (and in total cross sections), it becomes comparable with (and dominant over) the conventional Pomerons in the differential cross sections at large $`|t|`$. Since its contribution does not qualitatively change the shape of the cross section, it is difficult to extract the $`f_1`$-exchange from the cross section alone. Therefore, we should study the polarization observables for its possible manifestation.

Our analysis of single spin observables shows that the one which is most sensitive to the $`f_1`$-trajectory is the beam asymmetry

$`\mathrm{\Sigma }_x{\displaystyle \frac{\mathrm{Tr}[T\sigma _\gamma ^xT^{}]}{\mathrm{Tr}[TT^{}]}}={\displaystyle \frac{\sigma _{}\sigma _{}}{\sigma _{}+\sigma _{}}},`$ (14)
where $`()`$ corresponds to a photon linearly polarized beam along (perpendicular) to the $`J/\psi `$ production plane. The prediction for $`\mathrm{\Sigma }_x`$ is shown in Fig.2. All the calculations shown below are done for $`W=94`$ GeV. The results shown are for both with (solid curve) and without (dot-dashed curve) $`f_1`$ contribution. The pure Pomeron part results in monotonic decreasing of $`\mathrm{\Sigma }_x`$ from 0 to some negative value depending on $`t`$. Solving Eq.(12) using Eq.(9), for pure $`f_1`$ channel at high energy and small vector meson production angle $`\theta `$, we get,

$`\mathrm{\Sigma }_x^{f_1}={\displaystyle \frac{s\mathrm{sin}^2\frac{\theta }{2}}{s\mathrm{sin}^2\frac{\theta }{2}+M_V^2}}>0,`$ (15)
which results in strong increase of $`\mathrm{\Sigma }_x`$ up to the large positive values exhibiting non-monotonic $`t`$-dependence. To further facilitate the experimental test of our predictions, we have also investigated two other spin observables. The first observable is the parity asymmetry defined as

$`P_\sigma {\displaystyle \frac{\sigma ^N\sigma ^U}{\sigma ^N+\sigma ^U}}=2\rho _{11}^1\rho _{00}^1,`$ (16)
where $`\sigma ^N`$ and $`\sigma ^U`$ are the cross sections due to the natural and unnatural parity exchanges and $`\rho _{\lambda \lambda ^{}}^i`$ are the vector meson spin density matrices. In the region where the natural parity exchange Pomeron part is dominant, one expects $`P_\sigma =1`$. Thus, any deviation from this value will be due to the unnatural parity $`f_1`$-exchange. Fig. 3 shows the corresponding prediction. One can see strong $`t`$ dependence of $`P_\sigma `$ which varies from 1 to $`1`$ in a relatively small interval of $`t`$.

The second one is the double beam-target asymmetry which is defined as ,

$$C_{BT}A_{LL}=\frac{d\sigma ()d\sigma ()}{d\sigma ()+d\sigma ()},$$
(17)
where the arrows denote the relative orientations of the proton and photon helicity, respectively. For both pure natural and unnatural parity exchange amplitudes, $`C_{BT}`$ is exactly zero. The deviation from this value is due to the interference between natural (Pomeron) and unnatural ($`f_1`$-exchange) parts of the total amplitude. The corresponding calculation is shown in Fig. 4. It can be seen clearly that at some region of $`t`$, the beam-target asymmetry deviates from zero considerably. This deviation is proportional to $`\mathrm{Re}(T_PT_{f_1})`$ and can be used to disentangle and identify the relative strengths of the two parts of the amplitude in the data.

In summary, we have analyzed the possible manifestation of the $`f_1`$ trajectory in $`J/\psi `$ elastic photo-production at high energy and large momentum transfers. The trajectory parameters are taken from the prediction of Ref.  and strength $`g_{Vf_1\gamma }`$ is estimated using . It is found that the $`f_1`$ exchange contribution can significantly influence the differential cross section at large $`|t|`$. We have presented predictions showing the $`f_1`$ effect on several spin observables. In particular, we have shown notable effects in single beam and double beam-target and parity asymmetries. It should be noted that our prediction is based on the assumption of the zero slope of the $`f_1`$ trajectory discussed in Ref. , in which two physical grounds for almost zero slope of the $`f_1`$ trajectory were given as (i) there is a nonvanishing contribution to the elastic $`pp`$ scattering at high energy due to such $`f_1`$ trajectory and (ii) $`f_1`$-meson exchange is deeply related to the axial anomaly which is associated with the instanton fluctuations whose space scale is much smaller than the scale of confinement. On the contrary, if the slope of $`f_1`$ trajectory is finite, we expect strong decreasing of its contribution both in $`s`$-dependence of the $`\sigma _{tot},d\sigma /dt(t=0)`$ and in $`t`$-dependence of $`d\sigma /dt(t=0)`$ at large $`s`$ and in consequence, any deviation from the standard Pomeron model for all the polarization observables discussed in our paper would disappear. Experimental test of our predictions will be a useful step towards understanding the structure of the vector meson photoproduction amplitude at large $`|t|`$. To test the effect of the $`f_1`$ trajectory contribution alone, it might be advantageous to use the other light vector meson $`(\rho ,\varphi )`$ photoproductions discussed in Ref. . However, if non-trivial interference effects in spin variables predicted in this work could be observed in experiment, it will not only be just a supplementary support for the $`f_1`$ trajectory but also a good test of the hard Pomeron contribution as well, because the hard Pomeron is important only in the $`J/\psi `$ photoproductions.

We should emphasize that the present investigation is the first step from the point of view of a dynamical treatment of the problem, as, for example, has been done for effective two-gluon model of Pomeron . Here, we do not know the sign of $`g_{Vf_1\gamma }`$ vertex, and have predicted only the absolute value of beam-target asymmetry. Moreover, it would be desired to generate the the corresponding form factor independently.

Finally, the experimental facilities like SPring-8, RHIC and HERA have plans to work more on polarized vector meson production. Our aim here is to find a clear non-trivial effect for the manifestation of the $`f_1`$ trajectory which reflects some new physics. Without the data on the experimental/accelerator resources which are needed to perform this measurement, it is rather difficult to give a detailed and complete analysis on the feasibility at present. We can say that very high luminosity, hopefully larger than $`10^{32}cm^2sec^1`$, might be necessary to actually test our prediction since the expected cross section at large $`|t|`$ regions where the predicted effect becomes evident, is quite small. We hope that such high luminosity could be realized in the forthcoming experiments so that we can have some good results to support this trajectory and also to establish its manifestation in QCD.

Acknowledgment

We gratefully acknowledge fruitful discussion with N.I. Kochelev and Y. Oh. One of us (T.M) would like to thank the Grant-in-Aid for Scientific Research, Ministry of Education, Science and Culture, Japan (No. 11694081), for financial support.