# Towards constraining axions with polarimetric observations of the isolated
neutron star RX J1856.5-3754

Aleksei Zhuravlev  Faculty of Physics, Lomonosov Moscow State University,
Vorobjevy Gory 1, Moscow, 119234, Russia  Roberto Taverna  Dipartimento di
Fisica e Astronomia, Università di Padova, via Marzolo 8, I-35131 Padova,
Italy  Roberto Turolla  Dipartimento di Fisica e Astronomia, Università di
Padova, via Marzolo 8, I-35131 Padova, Italy  MSSL-UCL, Holmbury St. Mary,
Dorking, Surrey RH5 6NT, UK<EMAIL_ADDRESS>
###### Abstract

Photon-axion mixing can create observable signatures in thermal spectra of
isolated, cooling neutron stars. Their shape depends on the polarization
properties of the radiation, which, in turn, are determined by the structure
of stellar outermost layers. Here we investigate the effect of mixing on the
spectrum and polarimetric observables, polarization fraction and polarization
angle, using realistic models of surface emission. We focus on RX
J1856.5-3754, the only source among the X-ray dim isolated neutron stars for
which polarimetric measurements in the optical band were performed. Our
results show that in the case of a condensed surface in both fixed and free-
ion limits, the mixing can significantly limit the geometric configurations
which reproduce the observed linear polarization fraction of 16.43%. In the
case of an atmosphere, the mixing does not create any noticeable signatures.
Complementing our approach with the data from upcoming soft X-ray polarimetry
missions will allow to obtain constraints on $g_{\gamma a}\sim 10^{-11}$ GeV-1
and $m_{a}\lesssim 10^{-6}$ eV, improving the present experimental and
astrophysical limits.

##  1 Introduction

Since their first introduction in the late 1970s, axion-like particles (or
simply axions) played a major role in modern extensions of the Standard Model
of particle physics. Apart from explaining the absence of CP violation in
strong interactions [Peccei & Quinn, 1977, Wilczek, 1978], they can provide a
solution to various problems in physics and astrophysics (see reviews by
Raffelt 1996, Marsh 2016, Di Luzio et al. 2020, Choi et al. 2020).

Primarily, axions manifest themselves through conversion into ordinary photons
(and vice versa) in external magnetic fields. The strength of the interaction
is mainly determined by the field geometry, coupling constant $g_{\gamma a}$,
and axion mass $m_{a}$. Currently, a significant portion of the axion
parameter space has been excluded by numerous experimental and astrophysical
studies. Among laboratory-based searches, the more stringent bounds on the
coupling constant over a wide axion mass range were obtained by the CERN Axion
Solar Telescope (CAST), $g_{\gamma a}\leq 6.6\times 10^{-11}$ GeV-1 (95% C.L.,
Anastassopoulos et al. 2017). Planned missions such as IAXO [Armengaud et al.,
2014] promise a tremendous improvement over CAST, and have the potential of
reaching down to $g_{\gamma a}\sim 10^{-12}$ GeV-1.

Various astrophysical searches typically attempt to find manifestations of
axions in the radiation from sources observed through regions permeated by
large-scale magnetic fields, such as in the case of galaxy clusters [Ajello et
al., 2016, Berg et al., 2017, Chen & Conlon, 2018, Libanov & Troitsky, 2020],
blazars and quasars [Fairbairn et al., 2011, Galanti et al., 2019], supernovae
[Payez et al., 2015]. Although in most cases the original spectrum of the
source is not precisely known, the presence of mixing should result in the
appearance of distinctive features. Their exact shape is determined by the
choice of axion parameters and configuration of the external magnetic field.
The possible constraints from such studies expand those of CAST, covering the
region $g_{\gamma a}\sim 10^{-12}$ GeV-1.

Another promising approach is to search for manifestations of axions in
strongly magnetized compact objects, such as magnetic white dwarfs (mWDs) and
neutron stars (NSs). First, axions may be produced thermally inside their
cores with energies from several keV to several hundred keV through the
nucleon or electron bremsstrahlung processes. Subsequently, they can escape
the interior due to their feeble interaction with matter and turn into photons
in the magnetic field surrounding the star, creating signatures in its hard
X-ray spectrum. The absence of such features can be used to impose constraints
on $\left(m_{a},g_{\gamma a}\right)$, the neutron-axion coupling constant
$g_{na}$ [Raffelt, 1996, 1990, Fortin & Sinha, 2018, Sedrakian, 2019, Fortin &
Sinha, 2019, Lloyd et al., 2021], and the electron-axion coupling constant
$g_{ea}$ [Dessert et al., 2021]

In the case of white dwarfs, surface photons may undergo conversion as they
propagate through the encompassing gaseous atmosphere [Gill & Heyl, 2011].
This could lead to an enhancement of the observed degree of linear
polarization, which is quite low in most white dwarfs. In particular, for mWD
PG 1015+014 the induced polarization degree exceeds the measured value of 5%
at $g_{\gamma a}\geq 2\times 10^{-11}$ GeV-1. Constraints from other mWDs with
stronger magnetic fields can further expand this region, provided that the
atmospheric plasma density is sufficiently low.

In the case of neutron stars, axions can create a number of distinctive
features. First, their magnetospheres can be a site for the conversion of dark
matter axions [Pshirkov & Popov, 2009, Foster et al., 2020]. Second, spectra
of neutron stars typically include a thermal component, and surface photons
may oscillate into axions as they propagate through the encompassing
magnetosphere [Lai & Heyl, 2006, Perna et al., 2012]. Assuming the spectral
energy distribution of an isotropic blackbody with uniform surface temperature
and a dipolar magnetic field, Zhuravlev et al. 2021 (hereafter Paper I) have
shown that the conversion can reduce the optical O-mode flux by up to 30% at
the currently unconstrained values of $\left(m_{a},g_{\gamma a}\right)$. The
high-energy part of the spectrum remained unchanged in all cases. If the
surface radiation of neutron stars includes a considerable fraction of O-mode
photons, such manifestation of axions can create potentially observable
effects for $g_{\gamma a}\sim 10^{-12}$ GeV-1.

The polarization properties of thermal emission are primarily determined by
the structure of stellar outermost layers. It is commonly believed that the
cooling surface is covered with a gaseous atmosphere that reprocesses the
radiation coming from the underlying crustal layers [Harding & Lai, 2006, van
Adelsberg & Lai, 2006, Potekhin, 2014]. Its opacity for the X-mode photons is
substantially reduced by the strong magnetic field, while that for the O-mode
ones remains almost unaffected. Therefore, in this case, the mixing is
unlikely to create detectable signatures, since the observed photons will be
mostly polarized in the X-mode. On the other hand, a group of sources – the
so-called X-ray dim isolated neutron stars or XDINSs, also known as the
“Magnificent Seven” [see e.g. Turolla, 2009, for a review] – may have a liquid
or solid condensed surface due to a phase transition in the outermost layers
induced by the strong magnetic field and relatively low temperature [Lai &
Salpeter, 1997, Lai, 2001, Medin & Lai, 2007, Suleimanov et al., 2010]. As a
result, no atmosphere will be present, and the polarization properties of the
radiation will be mainly determined by the emissivity of the underlying
metallic layer. Since the latter is of the same order in both X- and O-mode
over a broad energy range, the emission of a condensed surface represents a
promising place to search for manifestations of axions.

Among the “Magnificent Seven”, the best source for such study is RX
J1856.5-3754 (RX J1856). It is the brightest and nearest of XDINSs, with
$V=25.58$, a nearly $\lambda^{-4}$ optical-UV SED [van Kerkwijk & Kulkarni,
2001, Kaplan et al., 2011] and an X-ray spectrum well modeled by two blackbody
components [Sartore et al., 2012]. The period $P\sim 7\text{ s}$ and period
derivative $\dot{P}\sim 3\times 10^{-14}\text{ s s}^{-1}$ translate into a
spindown magnetic field $B\sim 1.5\times 10^{13}$ G [van Kerkwijk & Kaplan,
2008].

More importantly, RX J1856 is the only XDINS for which polarimetric
observations in the optical band were performed with the Very Large Telescope
(VLT). Mignani et al. [2017] measured a phase-averaged polarization degree
$\hat{\Pi}_{\mathrm{L}}=16.43\%\pm 5.26\%$ and a phase-averaged polarization
position angle $\hat{\chi}_{p}=145.39^{\circ}\pm 9.44^{\circ}$, computed east
of the North Celestial Meridian. As shown by González Caniulef et al. 2016
(hereafter Paper II), polarimetric measurements can be effectively used to
determine the structure of stellar outermost layers. The two models produce
different polarization patterns at infinity – in the case of an atmosphere,
the radiation appears much more polarized than in the case of a condensed
surface. As a result, by combining the optical measurements with the data from
upcoming soft X-ray polarimetry missions, such as XPP [Jahoda et al. 2019, the
follow-up mission of IXPE, Weisskopf et al., 2016], it is potentially possible
to unambiguously determine the origin of thermal emission from RX J1856.

The goal of our paper is to describe how the presence of photon-axion
conversion affects the total flux and the polarization observables of a
gaseous atmosphere and a condensed surface, focusing on the case of RX J1856.
In Section 2, we provide an overview of the photon-axion interaction in strong
magnetic fields and summarize the main properties of both models of surface
emission. Section 3 describes our approach to calculating the polarization
observables and the simplified method we adopt to account for the effect of
photon-axion conversion. The results for RX J1856 are presented in Section 4.
The observability of axion effects and methods for deriving constraints on
their parameters using optical and X-ray polarimetry data are discussed in
Section 5.

##  2 Theoretical background

Here we briefly summarize the basic theory behind our calculations; for more
detailed explanations we refer to Papers I and II. We use Lorentz–Heaviside
units with $\hbar=c=1$ and $\alpha=e^{2}/4\pi\sim 1/137$.

###  2.1 Photon-axion mixing

The oscillations of an O-mode photon with energy $\omega$ and an axion in the
external magnetic field are described by the following system:

$\displaystyle\left(\omega+\begin{pmatrix}\Delta_{\parallel}+\Delta_{p}&\Delta_{M}\\\
\Delta_{M}&\Delta_{a}\end{pmatrix}-i\partial_{z}\right)\begin{pmatrix}A_{\parallel}\\\
a\end{pmatrix}=0,$ (1)

$\displaystyle\Delta_{a}$ $\displaystyle=-\frac{m_{a}^{2}}{2\omega},$
$\displaystyle\Delta_{M}$ $\displaystyle=\frac{1}{2}g_{\gamma a}B\sin\Theta,$
(2) $\displaystyle\Delta_{p}$ $\displaystyle=-\frac{\omega_{p}^{2}}{2\omega},$
$\displaystyle\Delta_{\parallel}$
$\displaystyle=\frac{1}{2}\omega\left(n_{\parallel}-1\right),$

where $A_{\parallel}$ and $a$ denote the amplitudes of the ordinary photon
state and the axion, respectively; $\omega_{p}^{2}=4\pi\alpha n_{e}/m_{e}$ is
the plasma frequency squared, the $z-$ axis is along the direction of photon
propagation, and $\Theta$ is the angle between the latter and the direction of
the local magnetic field [Raffelt & Stodolsky, 1988]. We assume that the
electron density is of the order of the Goldreich-Julian value,
$n_{e}=-\boldsymbol{\Omega}\cdot\boldsymbol{B}/e$, where $\boldsymbol{\Omega}$
is the stellar angular velocity [Goldreich & Julian, 1969]. The photon
refractive index $n_{\parallel}$ deviates from unity due to the presence of a
strong magnetic field [see e.g. Harding & Lai, 2006].

As shown in Paper I, for the currently allowed values of the coupling constant
$g_{\gamma a}$ only the first half-period of oscillations can be present in
the magnetosphere of a neutron star (weak-oscillation mode), i.e. photons
convert into axions and do not convert back. This is due to the fact that
either the mixing term is dominated by the QED term
$\left(\Delta_{M}\ll\Delta_{\parallel}\right)$, or the oscillation length
becomes much larger than the scale length of the stellar magnetic field. In
addition, only the low-energy radiation can undergo conversion, while the
X-ray flux always remains unchanged. For $B_{0}=10^{13}$ G, $g_{\gamma
a}=2\times 10^{-11}$ GeV-1, and $m_{a}=10^{-8}$ eV, 30% of the optical O-mode
photons (redshifted energy $\omega_{\infty}=1\text{ eV}$) do not reach the
observer and instead turn into axions. In the following, we express the
resulting change in the spectrum in terms of the differential modification
factor $\Lambda$, which is defined as the relative fraction of the initial
photon flux with a given energy, reaching the observer without conversion. In
particular, a 30% decrease in the optical radiation corresponds to
$\Lambda=0.7$.

###  2.2 Intrinsic polarization degree

The cooling surface of a neutron star is commonly modeled by a gaseous
atmosphere in radiative and hydrostatic equilibrium [see e.g. Potekhin, 2014,
for a review]. A quite general property is that the opacity for X-mode photons
below the electron cyclotron frequency is substantially reduced by the strong
magnetic field, while that for the O-mode ones remains almost unchanged [see
Harding & Lai, 2006, for a comprehensive overview of the physics in strong
magnetic fields]. As a result, the intensity of the X-mode photons becomes
much larger than that of the O-mode ones, and the emergent radiation appears
highly polarized. Here we restrict to the case of a fully ionized pure H
atmosphere and solve the radiative transfer equation using the numerical
method developed by Lloyd [2003, see also ]. The code has four input
parameters: the local magnetic field strength, the effective temperature, the
angle between the local magnetic field and the surface normal, and the surface
gravity. To compute the emergent intensity for a range of photon energies and
emission directions, a complete linearization technique for the two normal
polarization modes is employed in a plane-parallel slab, see Paper II.

The polarization pattern produced by a particular surface emission model can
be conveniently expressed in terms of the intrinsic polarization degree:

$\Pi_{\text{L}}^{\text{EM}}=\frac{F_{\text{X}}-F_{\text{O}}}{F_{\text{X}}+F_{\text{O}}},$
(3)

where $F_{\text{X,O}}$ is the monochromatic, phase-averaged flux in each mode,
obtained by integrating the intensity over the visible part of the star
surface. By definition, $\Pi_{\text{L}}^{\text{EM}}=1$ for the radiation 100%
initially polarized in the X-mode, $\Pi_{\text{L}}^{\text{EM}}=-1$ for that in
the O-mode, and $\Pi_{\text{L}}^{\text{EM}}=0$ for the unpolarized emission.
As shown in Paper II, the intrinsic polarization degree of an atmosphere is
$\Pi_{\text{L}}^{\text{EM}}\approx 0.87$ in the optical band and
$\Pi_{\text{L}}^{\text{EM}}\approx 0.99$ in the X–ray band for all viewing
geometries, so the radiation is dominated by the X-mode photons.

On the other hand, it has been suggested that at the low surface temperature
$\left(\lesssim 100\text{ eV}\right)$ and strong magnetic field $\left(\gtrsim
10^{13}\text{ G}\right)$, typical of XDINSs, a phase transition between a
gaseous and a condensed state may occur. As a result, the underlying metallic
layer of the stellar crust will be exposed, and the emergent spectrum will be
determined by its radiative properties [Lai & Salpeter, 1997, Lai, 2001,
Turolla et al., 2004, Medin & Lai, 2007, Turolla, 2009, Potekhin, 2014]. Two
limiting cases are typically considered: free-ions, in which the effects of
the lattice on the interaction of the electromagnetic waves with ions are
neglected, and fixed-ions, in which the lattice interaction with
electromagnetic waves dominates the ion response. The true spectral properties
are expected to lie between these two limits.

Following Paper II, we adopt the analytic descriptions developed by Potekhin
et al. [2012] to calculate the intensity in two modes. The reflectivity is
first obtained by applying Snell’s law at the interface between the vacuum and
the condensed phase. The emissivity $j_{\omega}$ is related to the latter via
Kirchhoff’s law, and the intensity is

$\begin{gathered}I_{\omega,\mathrm{O}}=j_{\omega,\mathrm{O}}B_{\omega}(T),\\\
I_{\omega,\mathrm{X}}=j_{\omega,\mathrm{X}}B_{\omega}(T),\end{gathered}$ (4)

where
$B_{\omega}(T)=\omega^{3}/\left[4\pi^{3}\left(\exp\left(\omega/T\right)-1\right)\right]$
is the blackbody radiance.

Contrary to the atmospheric model, the emissivity of a condensed surface in
the X- and O-mode is of the same order over a broad energy range [Suleimanov
et al., 2010]. In the optical band, $\Pi_{\text{L}}^{\text{EM}}$ can reach
$-0.3$ in the free-ion and $-0.5$ in the fixed-ion limits for the most
favorable viewing geometries, with the emission being polarized mostly in the
O-mode. In the X-ray band, the radiation is almost unpolarized,
$|\Pi_{\text{L}}^{\text{EM}}|\lesssim 0.07$ in both limits.

##  3 The model

To calculate the effects of photon-axion conversion using realistic models of
surface emission, we combine the methods presented in Papers I and II. The
long computation time implied by such a direct approach can be significantly
reduced by making a number of simplifications which are discussed below.

We restrict to the case in which the stellar magnetic field is dipolar,

$\boldsymbol{B}^{\text{dip}}=\frac{B_{p}}{2}\left(\frac{R_{\rm
NS}}{r}\right)^{3}\begin{pmatrix}2f_{\text{dip}}\cos\theta\\\
g_{\text{dip}}\sin\theta\\\ 0\end{pmatrix},$ (5)

where $B_{p}$ is the polar field strength, $R_{\rm NS}$ is the stellar radius,
and $r,\theta$ are the radial coordinate and the magnetic colatitude,
respectively. The functions $f_{\text{dip }}$ and $g_{\mathrm{dip}}$ account
for the relativistic corrections [Muslimov & Tsygan, 1986, Page & Sarmiento,
1996]. The polar magnetic field is taken to be $B_{p}=10^{13}$ G, which is
compatible with the spin-down estimate for RX J1856 [van Kerkwijk & Kaplan,
2008]. The spin period is the measured one, $P=7$ s [Tiengo & Mereghetti,
2007], while the stellar radius and mass are taken as $R_{\rm NS}=12$ km,
$M_{\rm NS}=1.5M_{\odot}$, respectively. The surface temperature distribution
is that induced by a core-centered dipole, see Paper II for details.

###  3.1 Ray tracing

Thermal emission of a cooling neutron star is calculated using the ray tracing
method presented in Zane & Turolla 2006 (see also Taverna et al. 2015 and
Paper II). First, the source geometry is determined by choosing the angle
$\chi$ between the line of sight $\boldsymbol{\ell}$ and the spin axis
$\mathbf{p}$, as well as the angle $\xi$ between the magnetic (dipole) axis
$\mathbf{b}_{\mathrm{dip}}$ and the spin axis. Then, the monochromatic flux
detected by a distant observer is obtained by summing the contribution of each
surface element which is in view at a given rotational phase, taking into
account the effects of ray bending,

$F_{\nu}(\gamma)=\left(1-\frac{R_{\mathrm{s}}}{R_{\mathrm{NS}}}\right)\frac{R_{\mathrm{NS}}^{2}}{D^{2}}\int_{0}^{2\pi}d\Phi_{\mathrm{S}}\int_{0}^{1}I_{\nu}(\mathbf{k},\theta,\phi)du^{2}.$
(6)

Here $I_{\nu}$ is the specific intensity, which in general depends on the
photon frequency $\nu$, emission direction $\mathbf{k}$, and position on the
stellar surface, $R_{\rm s}$ is the Schwarzschild radius, $D$ the source
distance and $x=R_{\rm s}/R_{\rm NS}$. The intensity is naturally written in
terms of the polar angles $(\theta,\phi)$ of a co-rotating coordinate system
$(x,y,z)$, with the $z$-axis parallel to $\mathbf{b}_{\text{dip}}$ and the
$x$-axis orthogonal to both $\mathbf{p}$ and $\mathbf{b}_{\text{dip}}$.
However, the integration is performed in a fixed reference frame $(X,Y,Z)$,
with the $Z$-axis along the line of sight and the $X$-axis in the
$(\boldsymbol{\ell},\mathbf{p})$ plane, with the associated polar angles
$\left(\Theta_{\mathrm{s}},\Phi_{\mathrm{S}}\right)$. The transformations
linking $\left(\Theta_{\mathrm{s}},\Phi_{\mathrm{S}}\right)$ and
$(\theta,\phi)$ are given in Paper II. Finally, $u=\sin\bar{\Theta}$, and the
angles $\Theta_{\mathrm{S}}$ and $\bar{\Theta}$ are related by the “ray
tracing” integral

$\bar{\Theta}=\int_{0}^{1/2}\frac{dv\sin\Theta_{s}}{\left[(1-x)/4-(1-2vx)v^{2}\sin^{2}\Theta_{s}\right]^{1/2}}.$
(7)

###  3.2 Polarization radius

According to QED, the presence of a strong magnetic field makes the vacuum
around a neutron star anisotropic. The dielectric and magnetic permeability
tensors are modified by virtual electron-positron pairs, which impacts on the
polarization properties of radiation. As a linearly polarized electromagnetic
wave propagates in the magnetized vacuum near the star, its electric field can
instantly follow the direction of the external magnetic field. Up to the
adiabatic (polarization-limiting) radius,

$r_{p}\simeq 7.9\left(\frac{B_{p}}{10^{13}\text{
G}}\right)^{2/5}\left(\frac{\omega_{\infty}}{1\text{
eV}}\right)^{1/5}\left(\frac{R_{\rm NS}}{12\,{\rm
km}}\right)^{1/5}R_{\mathrm{NS}},$ (8)

the photons emitted at the surface maintain their initial polarization state
[Heyl & Shaviv, 2002, Heyl et al., 2003, Taverna et al., 2015]. Around
$r_{p}$, the coupling weakens and the photon electric field is no longer able
to properly follow the variation of the local magnetic field; at $r\gg r_{p}$,
its direction freezes. Using the adiabatic radius approximation, it is
possible to take into account the effects of vacuum polarization without
numerically integrating the wave equations along each photon trajectory. As
shown in Taverna et al. [2015], one can assume mode locking up to $r_{p}$, and
keep the direction of the wave electric field constant at $r>r_{p}$. As a
result, the polarization properties will be determined by the direction of the
local magnetic field at $r_{p}$, as well as by the intrinsic polarization
degree.

To obtain the polarization observables as detected by a distant instrument,
the Stokes parameters of individual photons at $r_{p}$ need to be rotated so
that they refer to the same coordinate system [see Taverna et al., 2015, for
details]. The collective Stokes parameters $I,Q,U$ are then simply the sum of
the individual parameters, and the observed polarization fraction and
polarization angle are given by

$\displaystyle\Pi_{\mathrm{L}}$ $\displaystyle=\frac{\sqrt{Q^{2}+U^{2}}}{I}$
(9) $\displaystyle\chi_{p}$
$\displaystyle=\frac{1}{2}\arctan\left(\frac{U}{Q}\right)$ (10)

###  3.3 Conversion radius

The main effects of the photon-axion interaction can be taken into account
using the conversion-radius approximation. This method is based on several
simplifications which are discussed below.

First, we utilize the fact that the relative converted fraction of the O-mode
flux varies by no more than 2% between different photon trajectories due to
the weak-oscillation mode, which is present at all allowed values of
$g_{\gamma a}$. Having fixed the magnetic field geometry, we use the code of
Paper I to numerically integrate Equations (1) over a large number of ray
paths, taking into account the ordinary photons only, and obtain the average
converted fraction $\bar{\Lambda}$. Furthermore, we assume that the mixing
occurs instantly at a specific point, rather than gradually over a
magnetospheric region. We introduce the conversion radius $r_{c}$, defined as
the distance from the stellar surface within which 80% of the effect takes
place; for the optical radiation of RX J1856, $r_{c}\sim 200$–$300$ km. As a
result, to account for the presence of axions we can simply multiply the
O-mode flux of each photon ray by $\bar{\Lambda}$ at the conversion radius,
instead of performing the time-consuming integration.

The adiabatic-radius method for taking into account the vacuum polarization is
very similar to the conversion-radius approximation. In fact, these two steps
can be combined. The dipolar magnetic field does not change its direction
significantly between $r_{p}$ and $r_{c}\sim\left(2-3\right)r_{p}$; due to the
weak-oscillation mode, the slight depolarization in this region may be safely
neglected when calculating the effect of conversion. As a result, we can
simply assume that $r_{c}=r_{p}$, and simultaneously decrease the O-mode flux
and compute the Stokes parameters of the remaining photons at the adiabatic
radius. By using these simplifications, we can fully utilize the robustness of
the ray-tracing code and avoid solving the photon-axion propagation equations,
which significantly reduces the computational time.

##  4 Results

In this section we present our simulations of the photon-axion interaction in
the emission from a condensed surface in free/fixed ion limits, as well as a
magnetized hydrogen atmosphere. We typically use the axion parameters
$g_{\gamma a}=2\times 10^{-11}$ GeV-1 and $m_{a}=10^{-8}$ eV, which are below
the present experimental constraints obtained by CAST [Anastassopoulos et al.,
2017]. The corresponding average modification factor is shown in Figure 1.

We make specific reference to RX J1856, for which Mignani et al. [2017]
measured a phase-averaged polarization degree
$\hat{\Pi}_{\mathrm{L}}=16.43\%\pm 5.26\%$ and a phase-averaged polarization
position angle $\hat{\chi}_{p}=145.39^{\circ}\pm 9.44^{\circ}$ in the VLT-
FORS2 $\nu_{\rm HIGH}$ filter ($2.0-2.5$ eV). They were able to derive
constraints on the source geometry by comparing the observed values of
$\Pi_{\mathrm{L}}$ and of the X-ray pulsed fraction with the ones predicted by
different surface emission models, as functions of the two angles $\xi$ and
$\chi$ (see Section 2.2). In particular, it is either
$\left(\xi,\chi\right)\sim\left(5^{\circ},21^{\circ}\right)$ or
$\left(3^{\circ},52^{\circ}\right)$ for a fixed-ion condensed surface and
$\left(\xi,\chi\right)\sim\left(14^{\circ},3^{\circ}\right)$ for an
atmosphere.

Figure 1: Average modification factor $\bar{\Lambda}$ for $g_{\gamma
a}=2\times 10^{-11}$ GeV-1 and $m_{a}=10^{-8}$ eV, as a function of the photon
energy. Here the magnetic field geometry is an aligned dipole with
$B_{p}=10^{13}$ G and $\theta=60^{\circ}$.

We first consider the case of a condensed surface in the fixed-ion limit.
Figure 2 illustrates the optical spectra with and without conversion for
$\left(\xi,\chi\right)=\left(5^{\circ},21^{\circ}\right)$; the total flux in
the former case turns out to be reduced by around 15% with respect to that
computed in the absence of mixing. This is expected since the intrinsic
polarization degree is negative, and the radiation is dominated by the O-mode
photons for which conversion may occur (see Paper II). However, the effect may
be difficult to distinguish observationally. The optical counterparts of
XDINSs are faint, and there are considerable uncertainties in their spectral
measurements. In particular, for RX J1856 the change in total intensity by 15%
lies well within the errors van Kerkwijk & Kulkarni [2001]. As a result,
polarimetric observables should be our primary tool in detecting the presence
of axions, while the spectra are best used to validate our calculations.

Figure 2: Number flux (in arbitrary units) for the condensed surface in the
fixed-ion limit as a function of the photon energy: without (solid line) and
with mixing (dashed line). The source geometry is
$\left(\xi,\chi\right)=\left(5^{\circ},21^{\circ}\right)$ and the remaining
parameters are those of RX J1856.

Figure 3 illustrates the phase-averaged polarization fraction with and without
conversion, as well as the difference between the two cases. As we can see,
the predicted $\Pi_{\mathrm{L}}$ can indeed match the value measured in RX
J1856 for some geometrical configurations in the absence of mixing. However, a
20% decrease in the O-mode flux significantly reduces the overall
polarization, so that $\hat{\Pi}_{\mathrm{L}}$ is barely reached if mixing
occurs, and a much smaller portion of the $\xi-\chi$ plane lies within the
$1\sigma$ error region. This result can be interpreted as follows. Without
conversion, the amount of O-mode photons exceeds that of the X-mode ones,
leading to a moderate degree of polarization. When the conversion is taken
into account, the fluxes in the two modes become closer to each other, making
radiation less polarized. According to our calculations, for lower values of
$\bar{\Lambda}$ the radiation becomes dominated by X-mode photons, and the
observed $\Pi_{\mathrm{L}}$ starts increasing. Therefore, our main conclusion
is that for allowed values of $\left(m_{a},g_{\gamma a}\right)$, the
corresponding slight decrease in the O-mode flux can significantly reduce the
polarization fraction.

Figure 3: Contour plots of the phase-averaged polarization fraction in the
VLT-FORS2 $\nu_{\rm HIGH}$ filter ($2.0-2.5$ eV) for the condensed surface in
fixed-ion limit. The solid black line indicates the measured value of
$\bar{\Pi}_{\mathrm{L}}=16.43\%$, while the dashed black lines mark the
$1\sigma$ errors. The white line shows the combinations of
$\left(\xi,\chi\right)$ which reproduce the observed X-ray pulsed fraction of
1.3% [Tiengo & Mereghetti, 2007].

Since the measured $\hat{\Pi}_{\mathrm{L}}=16.43\%$ can still be reached at
some combinations of $\left(\xi,\chi\right)$ despite the presence of
conversion, we can additionally calculate the pulsed fraction induced by a
given magnetic field geometry and compare it with the observational data. For
RX J1856, its value is very low (1.3% in the X-ray band, Tiengo & Mereghetti
2007), which implies that either $\xi$ or $\chi$ is small $\left(\lesssim
15^{\circ}\right)$. As pointed out in Mignani et al. [2017], this method
allows to impose a constraint on the magnetic field geometry and thus exclude
parts of its parameter space. Also, note that the region with $\xi\geq\chi$ is
characterized by a vanishingly low degree of polarization, regardless of the
intrinsic one at the source. The reason for that is the frame rotation of the
Stokes parameters (see Taverna et al. 2015 and Paper II), which prevents us
from distinguishing the presence of conversion in this region. While not
relevant for RX J1856 since the measured optical polarization is non-
negligible, it may become an issue for other sources.

Our calculations of the phase-averaged polarization angle are illustrated in
Figure 4. Here the polarimeter frame was rotated around the line-of-sight by
an angle $\beta=\hat{\chi}_{p}$ with respect to the fixed frame $(X,Y,Z)$. As
shown in Taverna et al. [2015], for the radiation mostly polarized in the
O-mode $\chi_{p}$ becomes equal to the angle $\Psi$ between the polarimeter
reference axis and the $X$-axis, and to $\Psi+90^{\circ}$ for that in the
X-mode. Since the ordinary photons dominate in the region of $\xi-\chi$ plane
constrained by the measured $\hat{\Pi}_{\rm L}$ and X-ray pulsed fraction, we
need to rotate the polarimeter frame by $\hat{\chi}_{p}=145.39^{\circ}$ in
order to ensure consistency with the observations.

As it can be seen, the value of $\chi_{p}$ slightly changes only in the region
with $\xi\geq\chi$. This is due to the fact that the polarization angle
reflects the global direction of the photon electric field, which, in turn,
depends on the direction of the magnetic field at the adiabatic radius (see
Paper II). As a result, the observed $\chi_{p}$ should reflect the “phase-
averaged” direction of the magnetic field at $r_{a}$, and the jumps by
$90^{\circ}$ arise from an even slight predominance of ordinary photons over
extraordinary ones or conversely. Using the intrinsic polarization degree and
our calculations of the intensity, we can infer that the O-mode flux exceeds
that of the X-mode for most magnetic field geometries, both with and without
conversion. As a result, the mixing changes the value of $\chi_{p}$ only in
the region with a vanishingly small degree of polarization, where even a
single additional photon can modify the direction of the phase-averaged wave
electric field. Since $\chi_{p}$ remains constant for the source geometries
which reproduce the measured $\Pi_{\rm L}$, we cannot utilize the measurement
of polarization angle in the same way as that of the polarization fraction.

Figure 4: Contour plots of the phase-averaged polarization angle in the VLT
FORS2 $\nu_{\rm HIGH}$ filter for the condensed surface in fixed-ion limit.
The black line indicates the magnetic field geometries for which the measured
X-ray pulsed fraction is reproduced, while the blue lines mark those for
$\hat{\Pi}_{\rm L}$. The polarimeter frame was rotated by an angle
$\hat{\chi}_{p}=145.39^{\circ}$ with respect to the fixed frame for
consistency with the observations [Mignani et al., 2017].

In the case of a free-ion condensed surface, the emissivity in the O-mode is
closer to that in the X-mode than in the fixed-ion limit. The optical spectra
for $\left(\xi,\chi\right)=\left(35^{\circ},5^{\circ}\right)$ are presented in
Figure 5; the intrinsic polarization degree is still negative, but closer to
zero than in the fixed-ion case. The conversion decreases the total flux by
around 12%, which, again, can hardly be distinguished with observations.

The corresponding phase-averaged polarization fraction is illustrated in
Figure 6. As we can see, without conversion the free-ion model cannot
simultaneously reproduce the measured value of 16.43% and the pulsed fraction
of 1.3%. With 20% of the O-mode flux being removed, $\Pi_{\mathrm{L}}$ can
barely reach the $1\sigma$ error region. Also, note that at some geometries
the difference in polarization fraction is positive, which implies that the
modified flux consists primarily of the X-mode photons. This is better seen
from our calculations of the polarization angle, shown in Figure 7; the
polarimeter frame is still rotated by $\beta=\hat{\chi}_{p}$ with respect to
the fixed frame. Without conversion, $\chi_{p}$ is equal to $145.39^{\circ}$
in the $1\sigma$ region of $\hat{\Pi}_{\rm L}$, and changes to $55.39^{\circ}$
(i.e. $\hat{\chi}_{p}+90^{\circ}$) when the conversion is taken into account.
This confirms confirms our conclusion that the flux becomes dominated by the
X-mode photons. To match the observed value of $145.39^{\circ}$ in the top-
left corner of the $\xi-\chi$ plane for the case with conversion, the
polarimeter frame should be rotated by $235.39^{\circ}$ (or simply by
$55.39^{\circ}$).

Figure 5: Same as Figure 2 but for the free-ion model with
$\left(\xi,\chi\right)=\left(35^{\circ},5^{\circ}\right)$.

Figure 6: Same as Figure 3 but for the free-ion model.

Figure 7: Same as Figure 4 but for the free-ion model. The polarimeter frame
was rotated by an angle $\hat{\chi}_{p}=145.39^{\circ}$ with respect to the
fixed frame.

Finally, in the case of an atmosphere, the effect of conversion was negligible
– the total intensity decreased by less than 5%. This is simply due to the
fact that the ordinary photons make up about 7% of the total flux for
$\Pi_{\text{L}}^{\text{EM}}=0.87$, and removing its fraction cannot produce
observable results. For this reason, we do not show the plots of polarization
fraction and angle, since the former changes by no more than 1%.

##  5 Discussion

In this work we have calculated the effect of photon-axion mixing on thermal
emission of RX J1856.5-3754 in the optical band, using realistic models of
surface layers. The condensed surface in both fixed and free-ion limits turned
out to be a promising target for detecting the change in phase-averaged
polarization fraction, which can be significantly modified by a slight
decrease in the O-mode flux. In the case of a gaseous atmosphere, the
conversion had no measurable effects.

Optical polarimetry alone is clearly insufficient to determine the origin of
surface emission from RX J1856, since both an atmosphere and a condensed
surface can reproduce the observed $\Pi_{\text{L}}$ and $\chi_{p}$ even if
mixing is at work. This degeneracy can be effectively removed using
polarimetric observations at higher energies. In fact, as pointed out in Paper
II, the two models of stellar surface emission produce different polarization
patterns in the soft X-ray band. For a condensed surface in both limits,
$\Pi_{\mathrm{L}}$ is small for all possible source geometries, but non-
negligible for some combinations of $\left(\xi,\chi\right)$, see Figure 8. On
the other hand, for an atmosphere $\Pi_{\mathrm{L}}$ is high for $\xi<\chi$,
and low for $\xi\geq\chi$ due to the geometrical effects, see Section 4 and
Paper II. However, for RX J1856, the latter case is ruled out because such
geometries do not reproduce the measured $\hat{\Pi}_{\mathrm{L}}$ in the
optical band; as a result, in the X-rays the polarization degree of an
atmosphere should always be high.

Figure 8: Contour plots of the phase-averaged polarization fraction in the
soft X-ray band (0.12 - 0.39 keV) for the fixed- and free-ion condensed
surface, respectively. The photon-axion interaction has no effect at these
energies.

Since thermal high-energy photons cannot convert into axions (see Paper I),
measuring the soft X-ray $\Pi_{\mathrm{L}}$ for RX J1856 can provide a simple
method of constraining a range of axion parameters. If the polarization is
significant, the emission is most likely atmospheric and the mixing will have
no effect. On the other hand, if $\Pi_{\mathrm{L}}$ is low, the emission
should originate from a condensed surface. In this case, the constraints can
be obtained using the approach described in Sections 3 and 4. For a given set
of $\left(m_{a},g_{\gamma a}\right)$, we determine the converted fraction of
optical O-mode photons and calculate the polarization fraction in free and
fixed-ion limits. Then, we attempt to find a combination of
$\left(\xi,\chi\right)$ which simultaneously reproduces the high and low-
energy $\Pi_{\mathrm{L}}$, as well as the pulsed fraction of 1.3%, taking into
account the observational errors. If there will be no such geometry in both
limits, we can infer that the considered values of $\left(m_{a},g_{\gamma
a}\right)$ lead to an unrealistic degree of conversion and can be excluded.

The effect of conversion on the optical polarization fraction at lower values
of $g_{\gamma a}$ can be estimated as follows. As discussed in Paper II, the
$\xi-\chi$ parameter space can be divided into two areas: $\xi<\chi$, where
the observed $\Pi_{\mathrm{L}}$ is mainly determined by the intrinsic
polarization at the source, and $\xi\geq\chi$, where the latter is distorted
by the geometrical effects. Since for RX J1856 only the region $\xi<\chi$ is
of interest, we can neglect the geometrical effects and calculate
$\Pi_{\mathrm{L}}$ analytically using the intrinsic polarization degree, which
is simply the relative difference between the fluxes in two modes (Equation
3). As a result, the ratio $F_{O}/F_{X}=1.39$ is required to match the
measured $\hat{\Pi}_{\mathrm{L}}=16.43\%$. The polarization is reduced to
$\hat{\Pi}_{\mathrm{L}}-1\sigma=11.17\%$ if the fluxes relate as
$F_{O}/F_{X}=1.11$, which can be achieved if 10% of the O-mode photons convert
into axions. For a surface magnetic field $B_{p}=10^{13}$ G, this effect is
reached at $g_{\gamma a}\geq 10^{-11}$ GeV-1 and $m_{a}\leq 10^{-6}$ eV, see
Paper I. Obviously, the actual exclusions will depend on the value of X-ray
polarization fraction of RX J1856, and may be greater for other XDINSs with
stronger magnetic fields [Turolla, 2009].

Polarimetric observations are not the only method of determining the origin of
thermal emission from XDINSs. Numerous studies have focused on the spectral
properties of an atmosphere and a condensed surface to explain the well-known
feature of the “Magnificent Seven” – their optical and ultraviolet
counterparts exceed the extrapolation of X–ray blackbodies at low energies by
a factor of 5 – 50 [Kaplan et al., 2011]. Some works restrict to a single
model of surface layers: a strongly magnetized atmosphere was suggested for RX
J1605.3+3249 [Mori & Ho, 2007], while the spectrum of RX J0720.4-3125 was
attributed to the emission of a bare condensed surface [Pérez-Azorín et al.,
2006]. However, recent studies prefer a combination of both, in which the
atmosphere is relatively thin and the spectrum of outgoing radiation is also
affected by the properties of metallic surface underneath it. Although the
mechanisms of creation of such atmosphere are still debated, the model
produced satisfactory results for RX J1308.6+2127 [Suleimanov et al., 2010,
Hambaryan et al., 2011], RX J1605.3-3249 [Pires et al., 2019, Malacaria et
al., 2019], RX J0720.4-3125 [Hambaryan et al., 2017].

In the case of RX J1856, X-ray data reveal different blackbody temperatures,
$T^{\infty}_{c}\sim 40$ eV and $T^{\infty}_{h}\sim 60$ eV [Sartore et al.,
2012], as well as the absence of absorption lines or other spectral features.
Ho et al. [2007] were able to model the broadband spectrum with a thin
hydrogen atmosphere above a condensed iron surface, although their magnetic
field estimate $B_{p}\approx 4\times 10^{12}\text{ G}$ somewhat contradicts
the spin-down value $B_{p}\approx 1.5\times 10^{13}\text{ G}$ obtained by van
Kerkwijk & Kaplan [2008]. Also, the model of Ho et al. [2007] imposes a
constraint on the magnetic field geometry, in addition to that of the X-ray
pulsed fraction. The light curve can be explained if either $\xi$ or $\chi$ is
small $\left(<6^{\circ}\right)$, while the other angle lies between
$20^{\circ}$ and $45^{\circ}$.

In principle, our analysis can be generalized to the thin atmosphere as well,
but implementing this model in the code of Paper II will require a state-of-
the-art approach. Still, the effect of conversion can be estimated on the
basis of the general properties of emission from a condensed surface covered
by a thin atmosphere, see e.g. Potekhin et al. [2012]. At variance with a
conventional atmosphere, optically thick from optical to X-ray energies, a
thin atmosphere is designed in such a way to be thick at low energies and then
become thin above $\approx 0.1$ keV because of the reduced free-free opacity
with increasing frequency. As a result, the emergent radiation in the soft
X-rays (and above) is that produced by the underlying condensate, while in the
optical band it is that of a standard atmosphere. Since the radiation of the
latter is predominantly polarized in the X-mode, we can infer that the mixing
is unlikely to create observable signatures in the case of a thin atmosphere.

In our approach, the only missing component is the high-energy polarization
data for RX J1856. Currently, a number of X-ray polarimetry missions are at an
advanced stage of development. The IXPE [Weisskopf et al., 2016], a NASA Small
Explorer (SMEX) mission expected to fly at the end of 2021, will measure
polarization in the $2-8$ keV energy range. Even though the IXPE bandpass is
at higher energies than those where the emission from XDINSs peaks, the latter
can be reached with its follow-up mission, the X-ray Polarization Probe
[Jahoda et al., 2019], which promises a sensitivity improvement by a factor of
$3-10$ over IXPE and an energy bandpass broadened from $2-8$ keV to $0.1-60$
keV. Other promising missions include eXTP [Zhang et al., 2019] and PiSoX
[Marshall et al., 2020].

Apart from the phase-averaged optical and X-ray polarimetry, no other
diagnostics seems to be as useful for determining the physical state of
stellar surface layers. Phase-resolved polarimetry, on the one hand, can
provide information on the phase dependence of $\Pi_{\mathrm{L}}$ and
$\chi_{p}$, which is particularly sensitive to the source geometry and surface
emission model [Taverna et al., 2015, 2020]. This makes it more descriptive
than phase-averaged measurements, which are largely degenerate with respect to
the angles $\left(\xi,\chi\right)$ since $\chi_{p}$ remains constant in large
regions of the $\xi-\chi$ plane (see Figures 4 and 7). However, in the optical
band, such observations are hardly possible due to the faintness of the
optical counterpart of RX J1856 ($V\sim 25.5$), even with next-generation
telescopes such as the European Extremely Large Telescope (ESO-ELT111See
https://elt.eso.org/.). In the X-rays, it is unclear whether the high-energy
tail of RX J1856 will be bright enough for phase-resolved measurements to be
performed by upcoming soft X-ray polarimetry missions with affordable exposure
times. Spectral fitting, on the other hand, is even more degenerate than
phase-averaged polarimetry, since it generally yields too many sets of
atmospheric and condensed surface models which can satisfactorily reproduce
the observed broadband spectra [see Paper II and Potekhin, 2014].

Our study can be extended to other XDINSs as well, since their low surface
temperatures and strong magnetic fields are sufficient both for the surface
layers to enter a condensation phase and for the mixing to have a significant
effect on the optical emission. However, the low-energy counterparts of the
remaining XDINSs are much fainter than that of RX J1856 [Kaplan et al., 2011],
and polarization measurements are beyond the capabilities of present
generation telescopes, like the VLT. Only the sensitivity improvement promised
by next-generation instruments such as ESO-ELT will allow to extend our study
to more XDINSs. Other classes of neutron stars (e.g. magnetars) do not appear
as suitable, since their optical emission, the origin of which is still
uncertain, is likely due to non-thermal processes in the magnetosphere [see
Turolla et al., 2015, for a review], and the contribution from thermal surface
emission is too faint for mixing to produce an observable effect.

## Acknowledgements

We thank Sergei Popov and Sergey Troitsky for numerous helpful comments. The
work of RT and RT is partially supported by the Italian MUR through grant
UNIAM (PRIN 2017LJ39LM).

## References

  * Ajello et al. [2016] Ajello, M., Albert, A., Anderson, B., et al. 2016, Phys. Rev. Lett., 116, 161101, doi: 10.1103/PhysRevLett.116.161101
  * Anastassopoulos et al. [2017] Anastassopoulos, V., Aune, S., Barth, K., et al. 2017, Nature Physics, 13, 584, doi: 10.1038/nphys4109
  * Armengaud et al. [2014] Armengaud, E., Avignone, F. T., Betz, M., et al. 2014, Journal of Instrumentation, 9, T05002, doi: 10.1088/1748-0221/9/05/T05002
  * Berg et al. [2017] Berg, M., Conlon, J. P., Day, F., et al. 2017, ApJ, 847, 101, doi: 10.3847/1538-4357/aa8b16
  * Chen & Conlon [2018] Chen, L., & Conlon, J. P. 2018, MNRAS, 479, 2243, doi: 10.1093/mnras/sty1591
  * Choi et al. [2020] Choi, K., Im, S. H., & Shin, C. S. 2020, arXiv e-prints, arXiv:2012.05029.  https://arxiv.org/abs/2012.05029
  * Dessert et al. [2021] Dessert, C., Long, A. J., & Safdi, B. R. 2021, arXiv e-prints, arXiv:2104.12772.  https://arxiv.org/abs/2104.12772
  * Di Luzio et al. [2020] Di Luzio, L., Giannotti, M., Nardi, E., & Visinelli, L. 2020, Phys. Rep., 870, 1, doi: 10.1016/j.physrep.2020.06.002
  * Fairbairn et al. [2011] Fairbairn, M., Rashba, T., & Troitsky, S. 2011, Phys. Rev. D, 84, 125019, doi: 10.1103/PhysRevD.84.125019
  * Fortin & Sinha [2018] Fortin, J.-F., & Sinha, K. 2018, Journal of High Energy Physics, 2018, 48, doi: 10.1007/JHEP06(2018)048
  * Fortin & Sinha [2019] —. 2019, Journal of High Energy Physics, 2019, 163, doi: 10.1007/JHEP01(2019)163
  * Foster et al. [2020] Foster, J. W., Kahn, Y., Macias, O., et al. 2020, Phys. Rev. Lett., 125, 171301, doi: 10.1103/PhysRevLett.125.171301
  * Galanti et al. [2019] Galanti, G., Tavecchio, F., Roncadelli, M., & Evoli, C. 2019, MNRAS, 487, 123, doi: 10.1093/mnras/stz1144
  * Gill & Heyl [2011] Gill, R., & Heyl, J. S. 2011, Phys. Rev. D, 84, 085001, doi: 10.1103/PhysRevD.84.085001
  * Goldreich & Julian [1969] Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869, doi: 10.1086/150119
  * González Caniulef et al. [2016] González Caniulef, D., Zane, S., Taverna, R., Turolla, R., & Wu, K. 2016, MNRAS, 459, 3585, doi: 10.1093/mnras/stw804
  * Hambaryan et al. [2017] Hambaryan, V., Suleimanov, V., Haberl, F., et al. 2017, A&A, 601, A108, doi: 10.1051/0004-6361/201630368
  * Hambaryan et al. [2011] Hambaryan, V., Suleimanov, V., Schwope, A. D., et al. 2011, A&A, 534, A74, doi: 10.1051/0004-6361/201117548
  * Harding & Lai [2006] Harding, A. K., & Lai, D. 2006, Reports on Progress in Physics, 69, 2631, doi: 10.1088/0034-4885/69/9/R03
  * Heyl & Shaviv [2002] Heyl, J. S., & Shaviv, N. J. 2002, Phys. Rev. D, 66, 023002, doi: 10.1103/PhysRevD.66.023002
  * Heyl et al. [2003] Heyl, J. S., Shaviv, N. J., & Lloyd, D. 2003, MNRAS, 342, 134, doi: 10.1046/j.1365-8711.2003.06521.x
  * Ho et al. [2007] Ho, W. C. G., Kaplan, D. L., Chang, P., van Adelsberg, M., & Potekhin, A. Y. 2007, MNRAS, 375, 821, doi: 10.1111/j.1365-2966.2006.11376.x
  * Jahoda et al. [2019] Jahoda, K., Krawczynski, H., Kislat, F., et al. 2019, arXiv e-prints, arXiv:1907.10190.  https://arxiv.org/abs/1907.10190
  * Kaplan et al. [2011] Kaplan, D. L., Kamble, A., van Kerkwijk, M. H., & Ho, W. C. G. 2011, ApJ, 736, 117, doi: 10.1088/0004-637X/736/2/117
  * Lai [2001] Lai, D. 2001, Reviews of Modern Physics, 73, 629, doi: 10.1103/RevModPhys.73.629
  * Lai & Heyl [2006] Lai, D., & Heyl, J. 2006, Phys. Rev. D, 74, 123003, doi: 10.1103/PhysRevD.74.123003
  * Lai & Salpeter [1997] Lai, D., & Salpeter, E. E. 1997, ApJ, 491, 270, doi: 10.1086/304937
  * Libanov & Troitsky [2020] Libanov, M., & Troitsky, S. 2020, Physics Letters B, 802, 135252, doi: 10.1016/j.physletb.2020.135252
  * Lloyd [2003] Lloyd, D. A. 2003, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0303561
  * Lloyd et al. [2003] Lloyd, D. A., Hernquist, L., & Heyl, J. S. 2003, The Astrophysical Journal, 593, 1024, doi: 10.1086/376589
  * Lloyd et al. [2021] Lloyd, S. J., Chadwick, P. M., Brown, A. M., Guo, H.-K., & Sinha, K. 2021, Phys. Rev. D, 103, 023010, doi: 10.1103/PhysRevD.103.023010
  * Malacaria et al. [2019] Malacaria, C., Bogdanov, S., Ho, W. C. G., et al. 2019, The Astrophysical Journal, 880, 74, doi: 10.3847/1538-4357/ab2875
  * Marsh [2016] Marsh, D. J. E. 2016, Phys. Rep., 643, 1, doi: 10.1016/j.physrep.2016.06.005
  * Marshall et al. [2020] Marshall, H. L., Heine, S. N. T., Garner, A., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11444, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114442Y, doi: 10.1117/12.2562811
  * Medin & Lai [2007] Medin, Z., & Lai, D. 2007, MNRAS, 382, 1833, doi: 10.1111/j.1365-2966.2007.12492.x
  * Mignani et al. [2017] Mignani, R. P., Testa, V., González Caniulef, D., et al. 2017, MNRAS, 465, 492, doi: 10.1093/mnras/stw2798
  * Mori & Ho [2007] Mori, K., & Ho, W. C. G. 2007, MNRAS, 377, 905, doi: 10.1111/j.1365-2966.2007.11663.x
  * Muslimov & Tsygan [1986] Muslimov, A. G., & Tsygan, A. I. 1986, Soviet Ast., 30, 567
  * Page & Sarmiento [1996] Page, D., & Sarmiento, A. 1996, ApJ, 473, 1067, doi: 10.1086/178216
  * Payez et al. [2015] Payez, A., Evoli, C., Fischer, T., et al. 2015, JCAP, 02, 006, doi: 10.1088/1475-7516/2015/02/006
  * Peccei & Quinn [1977] Peccei, R. D., & Quinn, H. R. 1977, Phys. Rev. Lett., 38, 1440, doi: 10.1103/PhysRevLett.38.1440
  * Pérez-Azorín et al. [2006] Pérez-Azorín, J. F., Pons, J. A., Miralles, J. A., & Miniutti, G. 2006, A&A, 459, 175, doi: 10.1051/0004-6361:20065827
  * Perna et al. [2012] Perna, R., Ho, W. C. G., Verde, L., van Adelsberg, M., & Jimenez, R. 2012, ApJ, 748, 116, doi: 10.1088/0004-637X/748/2/116
  * Pires et al. [2019] Pires, A. M., Schwope, A. D., Haberl, F., et al. 2019, A&A, 623, A73, doi: 10.1051/0004-6361/201834801
  * Potekhin [2014] Potekhin, A. Y. 2014, Physics Uspekhi, 57, 735, doi: 10.3367/UFNe.0184.201408a.0793
  * Potekhin et al. [2012] Potekhin, A. Y., Suleimanov, V. F., van Adelsberg, M., & Werner, K. 2012, A&A, 546, A121, doi: 10.1051/0004-6361/201219747
  * Pshirkov & Popov [2009] Pshirkov, M. S., & Popov, S. B. 2009, Soviet Journal of Experimental and Theoretical Physics, 108, 384, doi: 10.1134/S1063776109030030
  * Raffelt & Stodolsky [1988] Raffelt, G., & Stodolsky, L. 1988, Phys. Rev. D, 37, 1237, doi: 10.1103/PhysRevD.37.1237
  * Raffelt [1990] Raffelt, G. G. 1990, Physics Reports, 198, 1, doi: https://doi.org/10.1016/0370-1573(90)90054-6
  * Raffelt [1996] —. 1996, Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles (University of Chicago press)
  * Sartore et al. [2012] Sartore, N., Tiengo, A., Mereghetti, S., et al. 2012, A&A, 541, A66, doi: 10.1051/0004-6361/201118489
  * Sedrakian [2019] Sedrakian, A. 2019, Phys. Rev. D, 99, 043011, doi: 10.1103/PhysRevD.99.043011
  * Suleimanov et al. [2010] Suleimanov, V., Hambaryan, V., Potekhin, A. Y., et al. 2010, A&A, 522, A111, doi: 10.1051/0004-6361/200913641
  * Taverna et al. [2015] Taverna, R., Turolla, R., Gonzalez Caniulef, D., et al. 2015, MNRAS, 454, 3254, doi: 10.1093/mnras/stv2168
  * Taverna et al. [2020] Taverna, R., Turolla, R., Suleimanov, V., Potekhin, A. Y., & Zane, S. 2020, MNRAS, 492, 5057, doi: 10.1093/mnras/staa204
  * Tiengo & Mereghetti [2007] Tiengo, A., & Mereghetti, S. 2007, ApJ, 657, L101, doi: 10.1086/513143
  * Turolla [2009] Turolla, R. 2009, Isolated Neutron Stars: The Challenge of Simplicity (Berlin, Heidelberg: Springer Berlin Heidelberg), 141–163, doi: 10.1007/978-3-540-76965-17
  * Turolla et al. [2004] Turolla, R., Zane, S., & Drake, J. J. 2004, ApJ, 603, 265, doi: 10.1086/379113
  * Turolla et al. [2015] Turolla, R., Zane, S., & Watts, A. L. 2015, Reports on Progress in Physics, 78, 116901, doi: 10.1088/0034-4885/78/11/116901
  * van Adelsberg & Lai [2006] van Adelsberg, M., & Lai, D. 2006, MNRAS, 373, 1495, doi: 10.1111/j.1365-2966.2006.11098.x
  * van Kerkwijk & Kaplan [2008] van Kerkwijk, M. H., & Kaplan, D. L. 2008, ApJ, 673, L163, doi: 10.1086/528796
  * van Kerkwijk & Kulkarni [2001] van Kerkwijk, M. H., & Kulkarni, S. R. 2001, A&A, 378, 986, doi: 10.1051/0004-6361:20011272
  * Weisskopf et al. [2016] Weisskopf, M. C., Ramsey, B., O’Dell, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, 990517, doi: 10.1117/12.2235240
  * Wilczek [1978] Wilczek, F. 1978, Phys. Rev. Lett., 40, 279, doi: 10.1103/PhysRevLett.40.279
  * Zane & Turolla [2006] Zane, S., & Turolla, R. 2006, MNRAS, 366, 727, doi: 10.1111/j.1365-2966.2005.09784.x
  * Zhang et al. [2019] Zhang, S., Santangelo, A., Feroci, M., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 29502, doi: 10.1007/s11433-018-9309-2
  * Zhuravlev et al. [2021] Zhuravlev, A., Popov, S., & Pshirkov, M. 2021, Physics Letters B, 136615, doi: 10.1016/j.physletb.2021.136615


# Don’t Overlook the Grammatical Gender: Bias Evaluation for Hindi-English
Machine Translation

Pushpdeep Singh
National Institute of Technology, Hamirpur
Anu, Hamirpur, India
<EMAIL_ADDRESS>

###### Abstract

Neural Machine Translation (NMT) models, though state-of-the-art for
translation, often reflect social biases, particularly gender bias. Existing
evaluation benchmarks primarily focus on English as the source language of
translation. For source languages other than English, studies often employ
gender-neutral sentences for bias evaluation, whereas real-world sentences
frequently contain gender information in different forms. Therefore, it makes
more sense to evaluate for bias using such source sentences to determine if
NMT models can discern gender from the grammatical gender cues rather than
relying on biased associations. To illustrate this, we create two gender-
specific sentence sets in Hindi to automatically evaluate gender bias in
various Hindi-English (HI-EN) NMT systems. We emphasise the significance of
tailoring bias evaluation testsets to account for grammatical gender markers
in the source language.

Don’t Overlook the Grammatical Gender: Bias Evaluation for Hindi-English
Machine Translation

Pushpdeep Singh National Institute of Technology, Hamirpur Anu, Hamirpur,
India<EMAIL_ADDRESS>

##  1 Introduction

State-of-the-art NMT models develop gender bias (Prates et al., 2019) during
their training by picking up spurious correlations in their training data. Due
to their poor coreference resolution ability, they often rely on biased
associations with, say, occupation terms to disambiguate the gender of
pronouns (Stanovsky et al., 2019; Saunders et al., 2020). Also, in many cases,
these models translate gender-neutral sentences into gendered ones (Prates et
al., 2019; Cho et al., 2019) especially giving a ‘masculine default’
translation or one which is based on a stereotype. This problem also exists
for HI-EN Machine Translation (Ramesh et al., 2021). When put to use, such
systems can cause various harms (Savoldi et al., 2021). Thus, evaluating and
mitigating such biases from NMT models is critical to ensure fairness.

Prior research evaluating gender bias in machine translation has predominantly
centered around English as the source language (Stanovsky et al., 2019;
Bentivogli et al., 2020; Costa-jussà et al., 2020; Currey et al., 2022;
Rarrick et al., 2023). However, these evaluation methods or testsets don’t
seamlessly extend to other source languages, especially the ones with
grammatical gender. For instance, in Hindi, elements like pronouns,
adjectives, and verbs are inflected with gender. Nonetheless, prior studies in
other source languages often utilise gender-neutral sentences for bias
evaluation using TGBI score (Cho et al., 2019; Ramesh et al., 2021). In our
experiments with TGBI Evaluation, we found that it doesn’t expose gender bias
present in various NMT models.

Therefore, in this work, we propose to evaluate NMT models for bias using
sentences with grammatical gender cues of the source language. This allows us
to ascertain whether NMT models can discern the accurate gender from context
or if they depend on biased correlations. Using Hindi as the source language
in NMT, we develop two bias evaluation benchmarks: _OTSC-Hindi_ and _WinoMT-
Hindi_ and evaluate various HI-EN NMT models for gender bias. We highlight the
importance of creating such benchmarks for source languages with expressive
gender markers.

| IT | GT | MS | AWS
---|---|---|---|---
Sentence Set | $p_{m}$ | $p_{f}$ | $p_{m}$ | $p_{f}$ | $p_{m}$ | $p_{f}$ | $p_{m}$ | $p_{f}$
Female Speaker, Female Friend | 98.41 | 1.59∗ | 1.68 | 98.32∗ | 98.97 | 1.03∗ | 95.61 | 4.39∗
Female Speaker, Male Friend | 99.25∗ | 0.75 | 90.66∗ | 9.34 | 99.72∗ | 0.28 | 95.70∗ | 4.30
Male Speaker, Female Friend | 99.35 | 0.65∗ | 2.43 | 97.57∗ | 66.01 | 33.99∗ | 99.29 | 2.71∗
Male Speaker, Male Friend | 99.91∗ | 0.09 | 96.45∗ | 3.55 | 98.60∗ | 1.40 | 97.48∗ | 2.52
Table 1:  Evaluation using the OTSC-Hindi test set. Here, $p_{m}$ and $p_{f}$ are the percentages of sentences translated as male and female for the speaker’s friend. $*$ corresponds to the percentage of sentences translated into the true label for each sentence set. Bold values indicate the maximum percentage of sentences translated into a single gender class. | $Acc$ | $\Delta_{G}$ | $\Delta_{S}$ | $N$
---|---|---|---|---
IndicTrans | 48.9 | 48.5 | -0.1∙ | 6.2
Google Translate | 69.0⋆ | 10.6⋄ | -3.8 | 5.3
Microsoft Translator | 57.7 | 32.9 | 0.2∙ | 4.1
AWS Translate | 49.9 | 51.9 | -0.2∙ | 2.8
Table 2: Comparison of performance of various NMT Models on WinoMT-Hindi.
${\star}$ indicates significantly highest value, ${\diamond}$ indicates
significantly lowest value, ${\bullet}$ indicates near about values for _A_
cc, $\Delta_{G}$ and $\Delta_{S}$, respectively. _N_ represents percentage of
sentences translated into gender-neutral form.

##  2 Experimental Setup

We create two bias evaluation testsets in Hindi to evaluate gender bias:
Occupation Testset with Simple Context (OTSC-Hindi) and WinoMT-Hindi, which
account for diverse grammatical gender cues in Hindi. Using these, we evaluate
four HI-EN NMT Models: IndicTrans (IT) (Ramesh et al., 2022), Google Translate
(GT), Microsoft Translator (MS) and AWS Translate (AWS). Previous work by
Ramesh et al. (2021) focused only on the gender-neutral side of Hindi;
however, we demonstrate how we can account for gender-specified context for
better evaluation of bias. Code and data are publicly
available111https://github.com/iampushpdeep/Gender-Bias-Hi-En-Eval.

##  3 OTSC-Hindi

Similar to Escudé Font and Costa-jussà (2019), we create a Hindi version with
grammatical gender cues: “ m us kAPF smy s {jAntA,jAntF} h,{mrA,mrF} do-t
[occupation] kA kAm {krtA,krtF} h. ”(I have known [him/her] for a long time,
my friend works as a [occupation].) Unlike the English version, this template
specifies the gender of the speaker using a gender-inflected verb, i.e.
“jAntA(jaanta)”(m) and “jAntF(jaanti)”(f). The possessive pronoun is also
gender inflected based on the gender of the speaker’s friend: “mrA(mera)”(m)
and “mrF(meri)”(f) which decides the use of verb “krtA(karta)”(m) and
“krtF(karti)” (f). However, the English translation only specifies the gender
of speaker’s friend. Based on four possibilities in the template, we construct
these four sets of sentences.

##  4 WinoMT-Hindi

Stanovsky et al. (2019) composed a challenge set called WinoMT for evaluating
gender bias in NMT models. We contextualize this test set for the evaluation
of bias in HI-EN Translation by manually creating “WinoMT-Hindi”, which
consists of 704 WinoBias-like sentences Zhao et al. (2018) in Hindi but
modified to include gender cues of the language, mainly: gender-inflected
adjectives, postpositions, and verbs. We can mark the gender of the target
(English) by simply checking for the presence of male pronouns (he, him or
his) or female pronouns (she or her) in the translation. Interestingly, we
also observe that a few sentences are translated into a gender-neutral form.
For gender bias evaluation, we use the metrics _Acc_ , $\Delta_{G}$ and
$\Delta_{S}$ given by Stanovsky et al. (2019).

##  5 Evaluation Results

The results are shown in Table 1 and Table 2. Based on these results, most
models show heavy bias against the female gender and prefer masculine default
in translations. On WinoMT-Hindi, these models perform only as good as random
guessing, while Google Translate performs better at discerning gender using
the grammatical gender cues in the source. On the other hand, evaluation based
on TGBI (results not shown here) gives almost similar performance on all NMT
systems which shows that it is not good at exposing gender bias.

##  6 Conclusion and Future Work

We highlighted the need for contextualising NMT bias evaluation for non-
English source languages. In future, we plan to extend our evaluation beyond
following a particular template and be accommodating to all gender identities
(Devinney et al., 2022).

## References

  * Bentivogli et al. (2020) Luisa Bentivogli, Beatrice Savoldi, Matteo Negri, Mattia A. Di Gangi, Roldano Cattoni, and Marco Turchi. 2020.  Gender in danger? evaluating speech translation technology on the MuST-SHE corpus.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 6923–6933, Online. Association for Computational Linguistics.
  * Cho et al. (2019) Won Ik Cho, Ji Won Kim, Seok Min Kim, and Nam Soo Kim. 2019.  On measuring gender bias in translation of gender-neutral pronouns.  In _Proceedings of the First Workshop on Gender Bias in Natural Language Processing_ , pages 173–181, Florence, Italy. Association for Computational Linguistics.
  * Costa-jussà et al. (2020) Marta R. Costa-jussà, Pau Li Lin, and Cristina España-Bonet. 2020.  GeBioToolkit: Automatic extraction of gender-balanced multilingual corpus of Wikipedia biographies.  In _Proceedings of the Twelfth Language Resources and Evaluation Conference_ , pages 4081–4088, Marseille, France. European Language Resources Association.
  * Currey et al. (2022) Anna Currey, Maria Nadejde, Raghavendra Reddy Pappagari, Mia Mayer, Stanislas Lauly, Xing Niu, Benjamin Hsu, and Georgiana Dinu. 2022.  MT-GenEval: A counterfactual and contextual dataset for evaluating gender accuracy in machine translation.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 4287–4299, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.
  * Devinney et al. (2022) Hannah Devinney, Jenny Björklund, and Henrik Björklund. 2022.  Theories of “gender” in nlp bias research.  In _Proceedings of the 2022 ACM Conference on Fairness, Accountability, and Transparency_ , FAccT ’22, page 2083–2102, New York, NY, USA. Association for Computing Machinery.
  * Escudé Font and Costa-jussà (2019) Joel Escudé Font and Marta R. Costa-jussà. 2019.  Equalizing gender bias in neural machine translation with word embeddings techniques.  In _Proceedings of the First Workshop on Gender Bias in Natural Language Processing_ , pages 147–154, Florence, Italy. Association for Computational Linguistics.
  * Prates et al. (2019) Marcelo O. R. Prates, Pedro H. C. Avelar, and Luis Lamb. 2019.  Assessing gender bias in machine translation – a case study with google translate.
  * Ramesh et al. (2022) Gowtham Ramesh, Sumanth Doddapaneni, Aravinth Bheemaraj, Mayank Jobanputra, Raghavan AK, Ajitesh Sharma, Sujit Sahoo, Harshita Diddee, Mahalakshmi J, Divyanshu Kakwani, Navneet Kumar, Aswin Pradeep, Srihari Nagaraj, Kumar Deepak, Vivek Raghavan, Anoop Kunchukuttan, Pratyush Kumar, and Mitesh Shantadevi Khapra. 2022.  Samanantar: The largest publicly available parallel corpora collection for 11 Indic languages.  _Transactions of the Association for Computational Linguistics_ , 10:145–162.
  * Ramesh et al. (2021) Krithika Ramesh, Gauri Gupta, and Sanjay Singh. 2021.  Evaluating gender bias in Hindi-English machine translation.  In _Proceedings of the 3rd Workshop on Gender Bias in Natural Language Processing_ , pages 16–23, Online. Association for Computational Linguistics.
  * Rarrick et al. (2023) Spencer Rarrick, Ranjita Naik, Varun Mathur, Sundar Poudel, and Vishal Chowdhary. 2023.  Gate: A challenge set for gender-ambiguous translation examples.
  * Saunders et al. (2020) Danielle Saunders, Rosie Sallis, and Bill Byrne. 2020.  Neural machine translation doesn’t translate gender coreference right unless you make it.  In _Proceedings of the Second Workshop on Gender Bias in Natural Language Processing_ , pages 35–43, Barcelona, Spain (Online). Association for Computational Linguistics.
  * Savoldi et al. (2021) Beatrice Savoldi, Marco Gaido, Luisa Bentivogli, Matteo Negri, and Marco Turchi. 2021.  Gender bias in machine translation.  _Transactions of the Association for Computational Linguistics_ , 9:845–874.
  * Stanovsky et al. (2019) Gabriel Stanovsky, Noah A. Smith, and Luke Zettlemoyer. 2019.  Evaluating gender bias in machine translation.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 1679–1684, Florence, Italy. Association for Computational Linguistics.
  * Zhao et al. (2018) Jieyu Zhao, Tianlu Wang, Mark Yatskar, Vicente Ordonez, and Kai-Wei Chang. 2018\.  Gender bias in coreference resolution: Evaluation and debiasing methods.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers)_ , pages 15–20, New Orleans, Louisiana. Association for Computational Linguistics.