MNRAS 538, 2396–2407 (2025) https://doi.org/10.1093/mnras/staf435 
Advance Access publication 2025 March 17 
Neutrino beaming in ultraluminous X-ray pulsars as a result of 
gravitational lensing by neutron stars 
A. A. Mushtukov , 1 ‹ A. Y. Potekhin , 2 , 3 I. D. Markozov, 2 , 3 S. Nallan, 4 K. Kornacka, 5 I. S. Ognev, 6 
V. Kravtsov , 7 A. A. Dobrynina 6 and A. D. Kaminker 2 
1 Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK 
2 Ioffe Institute, Politekhnicheskaya 26, St Petersburg 194021, Russia 
3 Space Research Institute (IKI) of the Russian Academy of Sciences, Profsoyuznaya 84/32, Moscow 117997, Russia 
4 Carpe Diem Academy, 6712 Tannahill Drive, San Jose, 95120 CA, USA 
5 Christa McAuliffe Academy School of Arts and Sciences, 5200 SW Meadows Rd Ste 150, Lake Oswego, OR 97035, USA 
6 P.G. Demidov Yaroslavl State University, Sovietskaya 14, Yaroslavl 150003, Russia 
7 Department of Physics and Astronomy, University of Turku, FI-20014 Turku, Finland 
Accepted 2025 March 13. Received 2025 February 12; in original form 2025 January 1 
A B S T R A C T 
X-ray pulsars e xperiencing e xtreme mass accretion rates can produce neutrino emission in the MeV energy band. Neutrinos in 
these systems are emitted in close proximity to the stellar surface and subsequently undergo gravitational bending in the space 
curved by a neutron star. This process results in the formation of a distinct beam pattern of neutrino emission and gives rise to the 
phenomenon of neutrino pulsars. The energy flux of neutrinos, when averaged over the neutron star’s pulsation period, can differ 
from the isotropic neutrino energy flux, which impacts the detectability of bright pulsars in neutrinos. We investigate the process 
of neutrino beam pattern formation, accounting for neutron star transparency to neutrinos and gravitational bending. Based on 
simulated neutrino beam patterns, we estimate the potential difference between the actual and apparent neutrino luminosity. We 
show that the apparent luminosity can greatly exceed the actual luminosity, albeit only in a small fraction of cases, depending on 
the specific equation of state and the mass of the star. F or e xample, the amplification can exceed a factor of 10 for ≈0 . 05 per cent 
of typical neutron stars with the mass of 1 . 4 M . Strong amplification is less probable for neutron stars of higher mass. In the 
case of strange stars, a fraction of high-energy neutrinos can be absorbed, and the beam pattern, as well as the amplification of 
apparent neutrino luminosity, depends on neutrino energy. 
Key words: accretion, accretion discs – stars: neutron – X-rays: binaries. 
1
X  
b  
T  
∼  
g  
c  
p  
T  
e  
X  
l  
l  
X  
F
 
t  

A  
fl  
o  
o  
m  
(  
a  
c  
a  
1  
a  
a  
d  
U  
a  
e  
f  
1  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025 I N T RO D U C T I O N  
-ray pulsars (XRPs) are accreting neutron stars (NSs) in close
inary systems (see Mushtukov & Tsygankov 2022 for re vie w).
ypical field strength at the NS surface here is expected to be
10 12 G or even stronger. Such a strong magnetic field modifies the
eometry of accretion flow directing it towards small regions located
lose to the poles of an NS and dramatically influences elementary
rocesses of radiation/matter interaction (see Harding & Lai 2006 ).
he luminosity of XRPs is powered by the accretion process, the most
fficient mechanism of energy release. The apparent luminosity of
RPs co v ers about nine orders of magnitude. The lowest detected
uminosity is known to be ∼10 32 erg s −1 . The brightest XRPs show
uminosity  10 40 erg s −1 and belong to the class of ultraluminous
-ray sources (ULXs; see Bachetti et al. 2014 ; Israel et al. 2017 ;
abrika et al. 2021 for re vie w). 
The geometry of the emission regions at the NS surface is expected
o be dependent on the mass accretion rate (Basko & Sunyaev 1975 ). E-mail: ale xander.mushtuko v@physics.ox.ac.uk 
i  
a  
n  
Published by Oxford University Press on behalf of Royal Astronomical Socie
Commons Attribution License ( https:// creativecommons.org/ licenses/ by/ 4.0/ ), whit a relati vely lo w mass accretion rate (  10 17 g s −1 ), the accretion
ow reaches the stellar surface and is decelerated in the atmosphere
f an NS due to the Coulomb collisions, which leads to the formation
f hotspots located close to the magnetic poles of a star. At higher
ass accretion rates, the luminosity of an NS is sufficiently high
  10 37 erg s −1 ) to cause radiative force that stops accretion flow
bo v e the stellar surface. It leads to the formation of accretion
olumns – extended structures confined by a strong magnetic field
nd supported by the radiation pressure gradient (Wang & Frank
981 ; Mushtukov et al. 2015 ; Zhang, Blaes & Jiang 2022 ). At mass
ccretion rates exceeding ∼10 19 g s −1 , accretion columns can be
dv ectiv e, i.e. X-ray photons are confined inside a sinking region
ue to large optical thickness of the flow (Mushtukov et al. 2018b ).
nder this condition, the temperature of plasma can be as high as
 few hundred keV, which is sufficient to cause intense creation of
lectron–positron pairs (Mushtukov, Ognev & Nagirner 2019 ) and
urther neutrino emission due to their annihilation (Kaminker et al.
992 ). Thus, bright XRPs can manifest themselves as sources of
ntense neutrino emission, where the total energy flux released due to
ccretion is channelled into luminosity in photons and luminosity in
eutrinos. The intrinsic neutrino luminosity of an NS is maximal right© 2025 The Author(s). 
ty. This is an Open Access article distributed under the terms of the Creative 
ch permits unrestricted reuse, distribution, and reproduction in any medium, 
provided the original work is properly cited. 
Neutrino beaming in ULX pulsars 2397 
a  
i  
(  
b
l  
fi  
c
a  
a
m
i
o  
(  
a
e  
(
e  
i
b  
f  
a
i
p
a
o  
g
g
c  
c
p
a
2
2
W
f
fi  
a
r
fi
S
b
d
w  
c
A
M  
N
a  
o  
M
i  
1
w  
e  
s  
p
p
e
w  
t  
c  
f  
  
f
 
R  
a  
o  
a  
l  
w  
t  
i
 
P  
a  
s  
t  
T  
t
v
W  
A  
i  
p
o  
m  
S  
a
m  
f  
S  
T  
b  
B  
a  
B  
i  
 
o  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025fter the supernova explosion and then decreases rapidly with time: it
s expected to be < 10 37 erg s −1 ( < 10 36 erg s −1 ) after a few × 10 2 yr
10 3 yr) after the e xplosion (Yako vlev et al. 2005 ), which is well
elow the expected neutrino luminosity in bright ULXs. 
The higher the mass accretion rate and the total luminosity, the 
arger the fraction of energy released in the form of neutrinos (see
g. 1 in Asthana et al. 2023 ). The highest accretion rates (and
onsequently neutrino luminosities) are expected in ULX pulsars 
nd bright Be X-ray transients (Reig 2011 ) because of their high
pparent photon luminosities. Ho we ver, a high apparent luminosity 
ay not necessarily correspond to a high accretion rate, if radiation 
s strongly collimated by accretion disc winds, as suggested for these 
bjects by King, Lasota & Klu ´zniak ( 2017 ); see also Lasota & King
 2023 ). On the other hand, this hypothesis is still under debate; the
rguments that disfa v our strong collimation have been presented 
.g. by Mushtukov et al. ( 2021 ) and Mushtukov & Portegies Zwart
 2023 ). 
Due to the extragalactic nature of confirmed ULX pulsars, the 
xpected neutrino flux is very low and even expected to be signif-
cantly below the neutrino isotropic background in a MeV energy 
and (Asthana et al. 2023 ). Ho we ver, estimations of neutrino flux
rom ULX pulsars (Asthana et al. 2023 ) were performed under the
ssumption of isotropic neutrino emission, which can be violated by 
nitial non-isotropic emission as well as by gravitational bending of 
article trajectories. 
In this paper, we investigate neutrino beam pattern formation 
ccounting for neutrino propagation in a space curved by the gravity 
f a star and transparency of a star for neutrino emission. The
ravitational bending results in a difference between actual (initially 
enerated) L ν and apparent L ν, app neutrino luminosities. The latter 
an be different for different distant observers. On the basis of
alculated beam patterns we obtain expected distributions of neutrino 
ulsars o v er the amplification factor: 
 ν ≡ L ν, app 
L ν
. (1) 
 M O D E L  SET-UP  
.1 Equation of state and structure of a neutron star 
e assume the NS structure to be spherical. Appreciable deviations 
rom the spherical symmetry can be caused by ultrastrong magnetic 
elds ( B  10 17 G) or by rotation with ultrashort periods (less than
 few milliseconds; see e.g. Haensel, Potekhin & Yakovlev 2007 and 
eferences therein), but we will not consider such cases. 
The general static isotropic metric satisfying the Einstein 
eld equations can be written in the ‘standard form’ using the 
chwarzschild coordinates as follows [see e.g. section 8.1 in Wein- 
erg ( 1972 ) and section 23 in Misner, Thorne & Wheeler ( 1973 )]: 
 s 2 = −B( r ) d t 2 + A ( r ) d r 2 + r 2 (d θ2 + sin 2 θ d φ2 ) , (2) 
here t is the time coordinate, r , θ , and φ are the spherical polar
oordinates, 
 ( r) = 
(
1 − 2 GM r 
c 2 r 
)−1 
, B( r) = e 2 /c 2 , (3) 
 r is the gravitational mass inside a sphere of radius r , G is the
ewtonian constant of gravitation,  is the gravitational potential, 
nd c is the speed of light in a vacuum. Then the mechanical structure
f an NS is go v erned by four first-order differential equations for r ,
 r ,  , and the local pressure P as functions of baryon number βnside a given spherical shell (e.g. Richardson, van Horn & Savedoff
979 ): 
d r 
d β
= 1 
4 πr 2 n¯ 
√ 
A ( r) , (4) 
d M r 
d β
= ρ
n¯ 
√ 
A ( r) , (5) 
d  
d β
= G M r + 4 πr 
3 P /c 2 
4 πr 4 n¯ 
√ 
A ( r) , (6) 
d P 
d β
= −
(
ρ + P 
c 2 
)
d  
d β
, (7) 
here n¯ is the mean number density of baryons. We integrate these
quations from r = 0 and M r = 0 at the centre of the star outwards,
tarting from a predefined baryon density n¯ at the centre, until a
redefined mass density at the outer boundary ρb is reached. 
The boundary condition for the gravitational potential  is 
rovided by the Schwarzschild metric outside the star, 
 
2  b /c 2 = 1 − 2 GM 
Rc 2 
, (8) 
here R and M = M R are the stellar radius and mass, and  b is
he value of  at the stellar surface. Since the value of  at the
entre of the NS is not known in advance, we integrate equation ( 6 )
or a shifted potential ˜  ( β) =  ( β) −  (0), with the initial value
˜ 
 (0) = 0 at the centre of the star, and the value of the shift  (0) is
ound from equation ( 8 ) after the integration has been completed. 
We solve the set of equations ( 4 )–( 7 ) numerically by the classic
unge–Kutta method on a non-uniform grid in β, with variable steps
dapted to provide sufficient accuracy at each grid node for each
f the computed functions. We decrease each step until the desired
ccuracy is reached. In order to prev ent accurac y loss in the outer
ayers of the star, where β is nearly constant as a function of ρ or r ,
e use the difference ( βb − β) as an independent variable, βb being
he value of β at r = R, which is equal to the total number of baryons
n the NS. 
We limit ourselves to three equations of state (EoSs): APR (Akmal,
 andharipande & Rav enhall 1998 ), SLy4 (Douchin & Haensel 2001 ),
nd BSk24 (Pearson et al. 2018 ). These EoSs describe the ground
tate for the nucleon–lepton ( npeμ) composition of matter, which is
he most conserv ati v e assumption, without an y ‘e xotic’ constituents.
he APR EoS for the NS core is based on realistic ef fecti ve
wo- and three-nucleon interactions, which allow one to reproduce 
arious nucleon scattering data and the properties of light nuclei. 
e adopt the version of the APR EoS named A18 + δv + UIX ∗ in
kmal et al. ( 1998 ), which includes a relativistic boost correction,
n the parametrized form of Potekhin & Chabrier ( 2018 ). This
arametrization includes the NS crust described by the BSk24 EoS 
n top of the core described by the APR EoS. The SLy and BSk
odels are based on ef fecti ve nucleon–nucleon interactions of the
kyrme type, adjusted to reproduce the EoS of pure neutron matter
nd the experimental properties of heavy atomic nuclei. These EoS 
odels are unified : they are based on the same microscopic models
or the core and the crust. For the first (SLy) EoS family, we use the
Ly4 EoS in the parametrized form of Haensel & Potekhin ( 2004 ).
he more recent BSk interaction model is more complicated; it is
etter tuned to the recent collection of experimental nuclear data. The
Sk24 and BSk25 versions of this model, which are very similar,
ppear to be preferred, as discussed by Pearson et al. ( 2018 ). The
Sk24 EoS is substantially stiffer than the SLy4 EoS, and we choose
t as a representative example of relatively stiff and soft EoS models.
We also consider the possibility of the EoS of strange matter built
f the u, d , and s quarks (Witten 1984 ; Alcock, Farhi & OlintoMNRAS 538, 2396–2407 (2025) 
2398 A. A. Mushtuk o v et al. 
M
Figure 1. NS mass M versus radius R for four EoS models compared with 
theoretical and observational constraints and estimates. The solid curves 
labelled with numbers 1–3 show mass–radius relations for APR (1), SLy4 (2), 
and BSk24 (3) EoS models, respectively. The dashed–dotted line shows an 
example of mass–radius relation for a strange star. The curves are truncated 
at the hydrostatic stability limits, marked by the heavy dots. The dark and 
light grey shaded triangles are prohibited by general relativity and causality. 
The horizontal light-blue band corresponds to the mass estimate (with 1 σ
uncertainties) for PSR J1810 + 1744. The horizontal dashed lines with arrows 
mark the upper and lower 3 σ limits to the maximum NS mass, derived from 
observations. The error bars show the 1 σ confidence intervals in M and R, 
inferred from observations of five NSs in binary systems, marked by the 
letters from A to E, according to the legend (see text for details). 
1  
E  
a  
6  
o
 
c  
T  
d  
S  
g  
w  
s
 
t  
b  
i  
c  
(  
r  
f  
t  
u  
o  
o  
P  
1  
T  
d  
l  
P
f
 
e  
i  
e  
m  
2  
t  
p  
t  
a  
(  
 
1  
w  
A  
a  
T  
p  
C  
f  
(  
d  
I  
a  
f  
r  
(  
o  
e  
E  
j  
X  
n  
r
 
c  
m  
w
 
N  
r
a  
p  
s  
M
2
2
T  
t  
r  
u  
w  
f  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025986 ; Haensel, Zdunik & Schaefer 1986 ). For the strange matter
oS, we use the approximation proposed by Zdunik ( 2000 ) and
dopt the fiducial parameters in their paper: the bag constant B =
0 MeV fm −3 , the QCD coupling constant αc = 0, and the rest energy
f the strange quark m s c 2 = 100 MeV. 
The solid curves in Fig. 1 show gravitational mass M versus
ircumferential radius R of an NS for the three selected EoSs.
he dot–dashed curve displays M( R) for a strange star model. The
ark grey shaded area ( R < R Sch , where R Sch = 2 GM/c 2 is the
chwarzschild radius) is prohibited by general relativity. The entire
rey-shaded triangle is prohibited by general relativity combined
ith the condition that the speed of sound must be subluminal (e.g.
ection 6.5.7 of Haensel et al. 2007 ). 
The lower horizontal dashed line with the upward arrow marks
he lower 3 σ limit to the maximum NS mass, M max > 2 . 09, obtained
y Romani et al. ( 2022 ) jointly for seven most massi ve kno wn NSs
n binaries with white dwarfs (WDs). The light-blue horizontal band
orresponds to the most accurate of individual estimates for these NSs
 M = 2 . 13 ± 0 . 04 M  for PSR J1810 + 1744). These mass estimates
ely on an analysis of orbital light curves with a specific model
or heating of the WD surface by radiation from a pulsar, hence
hey can be model dependent. Pulsar mass estimates were obtained
sing the effects of general relativity (in particular the Shapiro delay
f the pulsar signal) appear to be less model dependent. Ho we ver,
ne of the largest estimates of this kind ( M = 2 . 01 ± 0 . 04 M  for
SR J0348 + 0432; Antoniadis et al. 2013 ) was recently revised to aNRAS 538, 2396–2407 (2025) 0 per cent lower value ( M = 1 . 806 ± 0 . 037 M ; Saffer et al. 2024 ).
he highest reliable pulsar mass estimates, based on the Shapiro
elay measurements, are currently M = 2 . 073 ± 0 . 069 M  with the
ower bound M > 1 . 95 M  at the 95.4 per cent confidence level for
SR J0740 + 6620 (Fonseca et al. 2021 ) and M = 1 . 908 ± 0 . 016 M 
or PSR J1614 −2230 (Arzoumanian et al. 2018 ). 
Rezzolla, Most & Weih ( 2018 ) used quasi-universal relations
xhibited by equilibrium solutions of rotating relativistic stars to
nfer constraints on the maximum NS mass from an analysis of the
lectromagnetic and gravitational wave signals from the double NS
erger GW170817. Their most conserv ati ve upper limit M max <
 . 59 M  is shown in Fig. 1 by the upper horizontal dashed line with
he do wnward arro w. It relies on the assumption that the merger
roduct in GW170817 has collapsed into a black hole. Ho we ver,
he kilonova produced in this event could also be explained in an
lternative scenario of NS stripping without black hole formation
Blinnikov et al. 2022 ). Thus the indicated limit is model dependent.
The vertical and horizontal error bars in Fig. 1 show the available
 σ confidence intervals in M and R, respectively, for the cases
here these uncertainties are not large (  15 per cent ). The labels
 and B mark such intervals for the bursters SAX J1810.8 −2609
nd 4U 1702 −429, according to the analysis by N ¨attil ¨a et al. ( 2017 ).
he label C corresponds to the nearest and brightest millisecond
ulsar PSR J0437 −4715 in an NS–WD binary system, according to
houdhury et al. ( 2024 ). Labels D and E mark the results obtained
or PSR J0740 + 6620 (Salmi et al. 2024 ) and PSR J0030 + 0451
Miller et al. 2019 ), respectively, using an analysis of the energy-
ependent thermal X-ray wav eform observ ed by the Neutron Star
nterior Composition Explorer ( NICER ). Despite the belief that this
pproach was ‘less subject to systematic errors than other approaches
or estimating neutron star radii’ (Miller et al. 2019 ), a subsequent
eanalysis, performed for PSR J0030 + 0451 by Vinciguerra et al.
 2024 ) with alternative hotspot models and using jointly the data
f NICER and XMM–Newton , resulted in substantially different
stimates, shown in Fig. 1 by the dashed error bars and marked as
 
′ and E ′′ , which demonstrate the strong model dependence. Other
oint mass and radius estimates obtained from spectral analyses of
-ray radiation of NSs exhibit similar model dependence and are
ot plotted here (e.g. Tanashkin et al. 2022 ; see also discussion and
eferences in Potekhin et al. 2020 ). 
Fig. 1 demonstrates that the selected EoS models are reasonably
ompatible with the available observational NS mass and radius esti-
ates, although the softest SLy4 EoS is only marginally compatible
ith the lower limits to M max . 
Solutions of the stellar structure equations for the three selected
S EoSs are illustrated in Fig. 2 . The upper and lower panels show,
espectively, mass density ρ and gravitational mass distributions M r 
s functions of r for the total mass of an NS equal to 1 . 4 M  (left
anels) or 2 M  (right panels). The dot–dashed line in the left panels
hows analogous distributions for the model of a strange star with
 = 1 . 4 M . 
.2 Neutrino emission and propagation 
.2.1 Neutrino opacities 
he distinct feature of neutrino propagation near an NS stems from
he fact that NSs are typically transparent to neutrino emissions of
elati vely lo w energy (Sawyer & Soni 1979 ; Haensel & Jerzak 1987 ),
nless their internal temperature T exceeds ∼10 10 K ( k B T ∼ 1 MeV,
here k B is the Boltzmann constant), which occurs only immediately
ollo wing a supernov a explosion (see e.g. figs 3 and 9 in Potekhin &
Neutrino beaming in ULX pulsars 2399 
(a)
(b)
(c)
(d)
Figure 2. Mass density (panels a and c) and mass (panels b and d) distributions inside an NS/strange star of a given total mass and EoS. Different curves are 
calculated for different EoS: APR (black solid), SLy4 (blue dotted), BSk24 (red dashed), and EoS for a strange star according to Zdunik ( 2000 ). Mass of an NS 
is fixed at M = 1 . 4 M  in panels (a) and (b), and at M = 2 M  in panels (c) and (d). The mass of a strange star was taken to be M = 1 . 4 M  in panel (a). 
C  
e
t  
c
m
s
t
n
K  
n
e
d
a
d
e  
o
o  
M
N
a  
h
1
p
e
 
a
n
a
n
T  
s  
e
λ
w  
s  
P
 
t
o  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025habrier 2018 ). In bright XRPs, it is expected that neutrinos are
mitted near the base of the accretion column. We assume that 
he magnetic field near the NS surface is dominated by the dipole
omponent, leading to neutrinos being initially emitted near the 
agnetic poles of the NS located diametrically opposite on the 
tellar surface. 1 The primary process of neutrino production involves 
he annihilation of electron–positron pairs, although a fraction of 
eutrinos can also be generated via the synchrotron process (see e.g. 
amink er et al. 1992 ; Mushtuk ov et al., in preparation). Emitted
eutrinos propagate both outside and inside the NS, which is cold 
nough to be nearly transparent to them. All emitted neutrinos un- 
ergo gravitational bending in curved space–time. Neutrino emission 
nd their subsequent propagation along curved trajectories form a 
istinct beam pattern. While the vast majority of neutrinos initially 
mitted at the NS surface are of electron fla v our, the composition
f the neutrino flux varies as it propagates, because of neutrino 
scillations. Neutrinos of very high energies E ν (abo v e sev eral
eV; see below) may experience scatterings in the core of the 
S, which influence their trajectories and, consequently, their final 
ngular distrib ution, b ut we will not consider such high energies
ere. 
 Note that in a few XRPs, non-dipole magnetic field structures have been 
roposed to explain observational data (see e.g. Postnov et al. 2013 ; Israel 
t al. 2017 ; Tsygankov et al. 2017 ; M ¨onkk ¨onen et al. 2022 ). 
(  
c
t
∝  
m
(  The primary processes go v erning neutrino opacity in dense matter
re neutrino–neutron scattering (see e.g. Shapiro & Teukolsky 1983 ), 
 + νe ,μ −→ n + νe ,μ, (9) 
nd neutrino absorption, 
 + νe −→ p + e −, n + νμ −→ p + μ−. (10) 
he neutrino mean free path in the elastic limit of neutrino–neutron
cattering ( 9 ) for non-degenerate nucleons can be estimated as (e.g.
quation 11.7.8 in Shapiro & Teukolsky 1983 ) 
sc , 0 ∼ 3 ρnuc 
ρ
(
1 MeV 
E ν
)2 
km , (11) 
here ρnuc ≈ 2 . 5 × 10 14 g cm −3 denotes the mass density at the
aturation number density of baryons n nuc ≈ 0 . 15 fm −3 (Horowitz,
iekarewicz & Reed 2020 ). 
Ho we ver, at temperatures T  10 8 K and densities ρ ∼ (1 –5) ρnuc ,
ypical for mature NSs (unlike, for instance, non-degenerate regions 
f a proto-NS), the neutrino energies most rele v ant for the ULXs
 E ν ∼ 0 . 1 –1 MeV) are large compared with temperature T but small
ompared with neutron Fermi energy E Fn . Under these conditions, 
he mean free path λsc increases compared with λsc , 0 by a factor 
 E Fn /E ν (see e.g. equation 7.2 in Iwamoto 1982 ). Then the neutrino
ean free path λsc in the elastic limit of neutrino–neutron scattering 
equation 9 ) can be estimated from equation (26) in Sawyer & SoniMNRAS 538, 2396–2407 (2025) 
2400 A. A. Mushtuk o v et al. 
M
(
λ
w
 
t  
a  
T  
p  
a  
p
n
w  
t  
1  
S  
n  
f  
fi  
t  
d
i  
E
 
n  
P  
e  
1  
e  
q
 
q  
M  
q  
i  
t
λ
w  
μ  
c  
e  
M  
c
2
T  
t  
(  
p  
s  
l  
I  
i  
t  
t  
c  
e
 
(  
a  
c  
(  
o
w  
p  
e  
F  
A  
d
b
w
 
(  
r  
a  
r  
t
(  
N  
t
2
T
L
w  
b  
t  
p  
t  
b  
s  
T  
t
 
n  
o  
r  
a  
t  
m  
d  
a  
a
T  
o  
h  
c  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025 1979 ) or (similar assessment) equation (7.2) in Iwamoto ( 1982 ), 
sc ≈ 800 
(
n nuc 
n n 
)2 / 3 (1 MeV 
E ν
)3 
km , (12) 
here n n is the number density of neutrons, so that λsc  R. 
The neutrino absorption process (equation 10 ) is forbidden in
he degenerate NS matter at p Fn  p Fp + p Fe (for the same reason
s the direct Urca process; cf. e.g. section 11.2 in Shapiro &
eukolsky 1983 ), where p F(n , p , e) are the Fermi momenta of neutrons,
rotons, and electrons, respectively. Thus in the bulk of a typical NS,
bsorption of neutrinos proceeds (in analogy with the modified Urca
rocesses) via the modified reactions: 
 + n + νe −→ p + n + e −, n + p + νe −→ p + p + e −, 
(13) 
ith the mean free paths λabs still longer than λsc at the low
emperatures (Sawyer & Soni 1979 ; Iwamoto 1982 ; Haensel & Jerzak
987 ). In this case, λabs is given by equation (16) in Sawyer &
oni ( 1979 ) or equation (20) in Haensel & Jerzak ( 1987 ). At the
uclear saturation density, we can estimate λabs using the convenient
ormula (equation 22) in Sawyer & Soni ( 1979 ), which agrees with
g. 1 in Haensel & Jerzak ( 1987 ) at T = 5 × 10 10 K and reproduces
he scaling law (equation 25) in Haensel & Jerzak ( 1987 ). At the
ensities and temperatures typical for the core of a mature NS, λabs 
ncreases with decreasing energy E ν from λabs ∼ (1 –3) × 10 6 km at
 ν = 1 MeV to λabs ∼ (1 –3) × 10 10 km at E ν = 0 . 1 MeV. 
In the quark (degenerate) matter, in contrast to the nucleon matter,
eutrino scattering is less efficient than absorption by quarks (e.g.
al & Dutt-Mazumder 2011 ), because the processes similar to
quation ( 10 ), e.g. d + νe → u + e −, are allowed (see e.g. Iwamoto
982 ). In this case, we have λqabs / λqsc  1, where λqabs and λqsc are
lectron neutrino mean free paths due to absorption and scattering in
uark matter, respectively. 
The absorption coefficient of non-degenerate neutrinos in the
uark matter can be assessed with equation (30) by Pal & Dutt-
azumder ( 2011 ). Let us assume E ν  πk B T , adopt the standard
uark colour factor C F = 4 / 3 and neglect higher order corrections
n this equation, retaining only the leading term, which is equi v alent
o equation (6.8) of Iwamoto ( 1982 ). Then we have 
qabs ≈ 14 km 
αs μu , 500 μd , 500 μe , 11 
(
1 MeV 
E ν
)2 
, (14) 
here αs is the strong coupling constant, μu , 500 = μu / (500 MeV ),
d , 500 = μd / (500 MeV ), μe , 11 = μe / (11 MeV ), μu , d , e being the
hemical potentials of the u and d quarks and of the electrons. For
stimates, following Sch ¨afer & Schwenzer ( 2004 ) and Pal & Dutt-
azumder ( 2011 ), we set αs = μu , 500 = μd , 500 = μe , 11 = 1, which
orresponds to densities ρ ≈ 6 ρnuc . 
.2.2 Neutrino trajectories 
o characterize the geometry of space–time around an NS, we employ
he static spherically symmetric metric in the standard form, equation
 2 ). It is a suitable approximation for NSs in XRPs with typical spin
eriods P spin  0 . 1 s (see e.g. Haensel et al. 2007 ). For a spherically
ymmetric NS in hydrostatic equilibrium, metric (equation 2 ) is
ocally similar to the Schwarzschild metric produced by mass M r .
n this metric, each trajectory of a freely propagating neutrino lies
n one and the same plane. Within this plane, we parametrize the
rajectory using polar coordinates, with r ≥ 0 and ϕ ∈ [0 , 2 π]. Since
he neutrino rest masses are negligibly small compared with theNRAS 538, 2396–2407 (2025) onsidered neutrino energies, we describe a neutrino trajectory by an
quation for a particle with zero rest mass. 
The deri v ation of the general equations of motion in metric
equation 2 ) can be found, e.g. in Weinberg ( 1972 ). In particular,
 trajectory of a particle in the equatorial plane (i.e. with the polar
oordinate θ = π/ 2) is given by equation (8.4.29) of Weinberg
 1972 ). For a massless particle, we should set the right-hand side
f this equation to zero, which leads to 
A ( r) 
r 4 
(
d r 
d ϕ 
)2 
+ 1 
r 2 
= 1 
B( r) b 2 , (15) 
here A ( r) and B( r) are determined by equation ( 3 ), and b has the
hysical meaning of the impact parameter, which is constant along
 very gi ven trajectory due to the angular momentum conservation.
or each trajectory, b is determined by the initial neutrino direction.
t the NS surface, it is related to the angle ζ between the radial
irection and neutrino trajectory as (see Appendix A ) 
 = R sin ζ√ 
1 − R Sch /R 
, (16) 
here R Sch = 2 GM/c 2 is the Schwarzschild radius. 
In the empty space outside an NS (at r > R), A ( r ) B( r ) = 1
see e.g. Weinberg 1972 , section 8.2). In this case, equation ( 15 )
educes to equation (25.55) in Misner et al. ( 1973 ), which describes
 photon trajectory in the Schwarzschild metric. Inside an NS (at
 < R), the functions A ( r) and B( r) are determined by the EoS
hrough the solution of the hydrostatic equilibrium equations ( 4 )–
 7 ). Consequently, trajectories of neutrinos propagating through an
S are contingent upon the mass distribution within the star and are
hus anticipated to vary for different EoSs. 
.3 Neutron star rotation and luminosity distribution 
he apparent luminosity of an NS can be determined as 
 ν, app = 4 πD 
2 
t av 
∫ t av 
0 
F ν( t ) d t , (17) 
here F ν( t) is variable neutrino energy flux density (as registered
y a distant observer), which varies with time t , D is a distance to
he compact object, and t av is a time interval for the averaging. In
ractice, the integration in equation ( 17 ) is performed o v er a long
ime interval ( t av  P spin ) because the mass accretion rate in X-ray
inaries is known to fluctuate o v er a wide range of time-scales, which
hould result in fluctuating pulse profiles in X-rays and neutrinos.
he ratio of the apparent and actual neutrino luminosity determines
he neutrino amplification factor a ν , equation ( 1 ). 
Apparent neutrino luminosity depends on actual neutrino lumi-
osity, neutrino beam pattern, and geometry of NS rotation in the
bserver’s reference frame. Rotation of an NS in the observer’s
eference frame is described by two angles: inclination i (i.e. the
ngle between the rotation axis and observer’s line of sight) and
he magnetic obliquity θB (i.e. the angle between the rotational and
agnetic axis of an NS). The flux is related to the neutrino flux
istribution in the reference frame of an NS, which depends on the
ngle ψ between the observer’s line of sight and NS magnetic axis
t a given phase ϕ p ∈ [0; 2 π] of NS rotation: 
cos ψ = cos i cos θB + sin i sin θB cos ϕ p . (18) 
he angles i and θB are typically unknown for the XRPs. Recent
bservation of X-ray polarization variable o v er NS spin period,
o we ver, sheds light on rotation geometry in some particular ac-
reting strongly magnetized NSs (Doroshenko et al. 2022 , 2023 ;
Neutrino beaming in ULX pulsars 2401 
T  
2  
o  
e
t  
a  
a  
n
i  
t  
M
3
O
t  
N
a
t
f  
d
3
T
e  
m
c
 
m
T
s  
o
t
i
b
a  
n
d  
i  
A
d
t
τ
w
f
λ
w  
q
g  
i  
i  
e  
A
T
b  
t
w
i  
(  
r  
i  
r  
s  
fi
 
S  
n  
c  
c  
b
3
T  
p
i  
t
v  
fi
w  
w  
t  
o  
N  
r
n  
p  
t  
c
p  
A  
o
b  
fl  
(
a
m
 
a  
e
a
i
a
θ
w  
e
I
f  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025sygankov et al. 2022 , 2023 ; Heyl et al. 2024 ; Malacaria et al.
023 ; Mushtukov et al. 2023 ; Forsblom et al. 2024 ), but features
f NS distribution o v er rotation parameters are still uncertain. To
stimate possible deviations of apparent neutrino luminosity from 
he actual one, we assume a random distribution of NSs o v er the
ngles i and θB , simulate neutrino pulse profiles for various i and θB ,
nd calculate theoretical distributions f ( a ν) of NSs o v er the apparent
eutrino luminosity amplification factors a ν . The technique used here 
s similar to the one applied to investigate distributions of XRPs o v er
he apparent luminosity in X-rays (see e.g. Mushtukov et al. 2021 ;
arkozov & Mushtukov 2024 ). 
 N U M E R I C A L  M O D E L  
ur numerical procedure consists of two stages. First, we compute 
he angular distribution of neutrino flux in the reference frame of an
S accounting for neutrino propagation along curved trajectories and 
bsorption inside a star. Then, using the computed angular distribu- 
ion, we simulate neutrino flux variability in the observer’s reference 
rame due to the rotation of an NS and calculate the theoretical
istribution of neutrino pulsars o v er the neutrino amplification factor. 
.1 Neutrino angular distribution 
o obtain the angular distribution of neutrinos, we specify neutrino 
nergy E ν , NS EoS, and mass, which give us NS radius and internal
ass distribution. Then we perform Monte Carlo simulations to 
alculate trajectories of 4 × 10 7 particles in each run. 
As an example, we consider the case when neutrino absorption is
uch more efficient than scattering, so scattering can be neglected. 
hen the simulation includes the following steps: 
(i) We start with the neutrino of energy E ν emitted from the 
urface of an NS near one of its magnetic poles. The initial direction
f particle motion is taken to be random and calculated under 
he assumption that the initial angular distribution of neutrinos is 
sotropic. 
(ii) We choose a random realization of the optical depth travelled 
y the particle before its absorption (the dimensionless free path) 
ccording to the formula τX = − ln X, , where X ∈ (0; 1) is a random
umber having a uniform distribution. 
(iii) We simulate a particle trajectory by numerically solving the 
ifferential equation ( 15 ) for a set of initial parameters, where the
nitial impact factors are calculated according to equation ( 16 ) (see
ppendix B ). If the trajectory crosses the star, we calculate the optical 
epth travelled by the particle, τν , by integrating along the simulated 
rajectory, 
ν( s) = 
∫ s 
0 
d s ′ 
λ( s ′ ) , (19) 
here s is the path-length along the particle trajectory. The mean 
ree path λ is given by 
= 1 
λ
−1 
abs + λ
−1 
sc 
 λabs , (20) 
hich accounts for both absorption ( λabs ) and scattering ( λsc ). These
uantities depend on the neutrino energy E ν , which experiences 
ravitational redshift (equation A11 ) while the particle mo v es along
ts trajectory, and the mass density ρ at each point along the path
nside the star, according to the estimates in Section 2.2 . The length
lement d s ′ in equation ( 19 ) is calculated using equation ( A10 ) in
ppendix A . he optical depth (equation 19 ) is a non-decreasing function, 
ounded from abo v e by some maximum value for every simulated
rajectory. If τX exceeds this maximum, the particle goes to infinity 
ithout absorption and we account for its final momentum direction 
n the simulated angular distribution function. Then we return to step
i) and start the simulation for the next particle. Otherwise, τν( s)
eaches τX at some point in the trajectory. In this case, the particle
s considered to be absorbed, we drop it from further consideration,
eturn to step (i) and start the simulation for the next neutrino. Thus
imulating trajectories of a large number of particles, we arrive at the
nal angular distribution of neutrinos in the NS reference frame. 
The described algorithm is rather general. Ho we ver, as argued in
ection 2.2.1 , the scattering and absorption do not noticeably affect
eutrino flux and can be neglected in a cold NS at the energies
onsidered in this study. On the other side, the neutrino absorption
an play a noticeable role in a strange quark star, as will be seen
elow in Fig. 7 . 
.2 Neutrino amplification factor 
o get the amplification factor (equation 1 ) for given rotation
arameters i and θB , we use pre-calculated neutrino flux distribution 
n the reference frame of a star (Section 3.1 ) and apply equation ( 21 )
o compute the theoretical pulse profile in neutrino emission. The 
ariation of the angle ψ between the line of sight and the magnetic
eld axis is given by 
cos ψ = sin i sin θB cos ϕ p + cos i cos θB , (21) 
here i is the inclination angle of the observer’s line of sight
ith respect to the rotational axis, and θB is the angle between
he magnetic field axis and the rotational axis. The phase angle
f pulsations, ϕ p ∈ [0 , 2 π), represents the rotational phase of the
S, with ϕ p = 0 corresponding to a reference position in the star’s
otation. Calculating the theoretical pulse profile, we assume that 
eutrinos are emitted from two magnetic poles of an NS which are
ositioned on opposite sides of the stellar surface. In our simulations,
he mass accretion rate is assumed to be constant. Under this
ondition, the pulse profile does not experience variations from one 
ulse period to another, and we use t av = P spin in equation ( 17 ).
veraging the neutrino flux variable o v er the NS spin period, we
btain the apparent neutrino luminosity (equation 17 ). Dividing it 
y the actual neutrino luminosity L ν , defined as the initial neutrino
ux integrated over the emission region as seen by a distant observer
i.e. corrected for the gravitational redshift), we obtain the neutrino 
mplification factor (equation 1 ) for any given rotation geometry, 
ass, EoS of the star, and emitted neutrino energy. 
To obtain the theoretical distribution of neutrino pulsars o v er the
mplification factors f ( a ν) we perform Monte Carlo simulations. In
ach simulation, we construct a neutrino pulse profile and calculate 
pparent neutrino luminosity for NS inclination, 
 = arccos (1 − 2 X 1 ) , (22) 
nd magnetic obliquity, 
B = πX 2 , (23) 
here X 1 , X 2 ∈ (0; 1) are random numbers. The constructed differ-
ntial distribution function f ( a ν) is normalized as ∫ ∞ 
0 
f ( a ν) d a ν = 1 . (24) 
n practice, it is useful to consider the cumulative distribution 
unction describing the fraction of objects amplified by a factor largerMNRAS 538, 2396–2407 (2025) 
2402 A. A. Mushtuk o v et al. 
M
Figure 3. Examples of neutrino trajectories calculated for particles emitted 
at the distance 10 6 cm from the Schwarzschild black hole with mass M = 
1 . 4 M . The grey circular area is the region surrounded by the Schwarzschild 
radius R Sch , while the blue dotted circle marks the distance r = 10 km from 
the centre. 
t
F
4
I  
t  
o  
n  
T  
a  
d  
a  
a  
p  
s
4
E  
m  
n  
a  
i  
W  
d  
f  
i  
f  
(
4
U  
d  
(  
s  
(a)
(b)
Figure 4. Examples of neutrino trajectories calculated for particles emitted at 
the NS surface for the case of different NS mass (and radius): (a) M = 1 . 4 M 
and (b) M = 2 M . In these calculations, the SLy4 EoS is assumed. For 
comparison, the same uneven steps for the initial colatitudes of trajectories 
are chosen in both panels. 
t  
o  
fl  
g  
n  
o  
c  
f  
4  
fi  
2  
n  
c  
o  
i  
b  
e  
m  
o  
t  
p
 
b  
a  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025han a ν : 
 ( a ν) ≡
∫ ∞ 
a ν
f ( x ) d x . (25) 
 N U M E R I C A L  RESULTS  
n this section, we demonstrate the results of our numerical simula-
ions of neutrino trajectories (Section 4.1 ), the angular distribution
f neutrino energy flux (Section 4.2 ), and theoretical distributions of
eutrino pulsars o v er the neutrino amplification factor (Section 4.3 ).
he gravitational bending of neutrinos propagating through a star is
ffected by the internal mass distribution. We analyse mass density
istributions calculated for three specific NS EoSs (see Section 2.1 ),
ssuming NS masses of 1 . 4 and 2 M , and for a strange star,
ssuming its mass of 1 . 4 M . In the latter case, to demonstrate the
ossible impact of neutrino absorption in a compact star, we perform
imulations for neutrinos of different energies. 
.1 Neutrino trajectories 
xamples of neutrino trajectories calculated in the Schwarzschild
etric near a 1 . 4 M  black hole are shown in Fig. 3 . Fig. 4 depicts
eutrino trajectories emitted from the surface of 1 . 4 M  (upper panel)
nd 2 M  (lower panel) NS. Unlike photons, neutrinos can penetrate
nto a compact star, where their trajectories follow geodesic paths.
ithin the star, neutrino trajectories are influenced by the mass
istrib ution, gra vitational potential, and pressure (see Appendix A
or details). The greater the mass of a star and the more concentrated
ts matter is toward the centre, the larger the deviation of particles
rom their original propagation direction, i.e. the deflection angle
see Fig. 5 ). 
.2 Angular distribution of neutrino flux 
tilizing the calculated neutrino trajectories, we derive the angular
istribution of neutrino energy flux in the reference frame of an NS
i.e. in a frame where the star does not rotate). Neutrino trajectories
tarted from the magnetic pole at the NS surface are curved and tendNRAS 538, 2396–2407 (2025) o converge in certain directions, leading to a significant amplification
f the neutrino energy flux in those areas. The directions of enhanced
ux depend on the NS mass and internal structure, which are
o v erned by the EoS. We note that the angular distributions of
eutrino flux al w ays show tw o peaks. The first peak is in the direction
pposite to the magnetic pole of a star that produces neutrinos, i.e. at
olatitude ∼π (see Fig. 6 ). Similar peaks have been reported earlier
or photons lensed in the gravitational field of an NS (see figs 3,
, and 9 in Riffert & Meszaros 1988 ; figs 8 and 9 in Kraus 2001 ;
gs 10–12 in Mushtukov et al. 2018a ; and fig. 6 in Mushtukov et al.
024 ). The second peak in the angular distribution corresponds to
eutrinos experiencing the maximal deflection ζ (see Fig. 5 ), at
olatitude ∼[2 π − ( ζ + ζ )] (see Fig. 6 ). The angular distribution
f neutrinos depends on both the EoS (compare different panels
n Fig. 6 ) and the mass of an NS (compare solid red and dotted
lack lines in Fig. 6 ). For NSs with smaller masses, the maximum
nhancement of the neutrino flux is more pronounced. In contrast,
ore massive NSs deflect neutrinos more strongly from their
riginal propagation direction. It results in an angular distribution
hat is closer to the isotropic one, albeit with distinct peaks still
resent. 
In the case of strange stars, unlike the NSs, neutrino absorption can
e noticeable. Nevertheless, angular distribution becomes strongly
nisotropic and neutrino energy flux can exceed the isotropic flux by
Neutrino beaming in ULX pulsars 2403 
(a)
(b)
Figure 5. The deflection angle for neutrinos emitted from the surface of an 
NS at different directions is given by angle ζ (see Section 2.2.2 ). Different 
curves are calculated for different EoSs: APR (solid black), SLy4 (dotted 
blue), BSk24 (dashed red), and for the case of the strange star (dashed–dotted 
green). The upper and lower panels are given for NS mass 1 . 4 and 2 M , 
respectively. 
m  
∼
fl  
4
U
d
f  
f  
i
 
a
e  
a
f  
e  
d
(
f
F
 
a
o
a  
(
(a)
(b)
(c)
Figure 6. Angular distribution of neutrino energy flux from one of the poles 
of an NS is given by a red solid (black dotted) line for 1 . 4 M  (2 M ) NS. 
The horizontal dotted line shows the level of the isotropic neutrino energy 
flux. Different panels correspond to different EoSs: (a) ARP, (b) SLy4, and 
(c) BSk24. In the case of relatively cold NSs under consideration, neutrino 
energy does not affect neutrino transfer within a star and, thus, the angular 
distribution. 
5
W  
e
a  
t
  
a  
n
b
b
N  
b
s
l
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025ore than an order of magnitude (see Fig. 7 ). Only at high energies
1 MeV, some fraction of neutrinos are absorbed, which reduce the 
ux directed towards a star (compare solid and dotted lines in Fig. 7 ).
.3 Luminosity function 
sing the calculated angular distributions of neutrino energy flux, we 
erive theoretical distributions of NSs o v er the neutrino amplification 
actor a ν and calculate the fraction F ( a ν) of NSs with amplification
actors abo v e specific values according to equation ( 25 ), as described
n Section 3.2 . These distributions are shown in Fig. 8 . 
One can see that the distributions of neutrino pulsars o v er the
mplification factor are relatively restricted: the majority of objects 
xhibit amplification factors within the interval a ν ∈ (0 . 5 , 10). The
nticipated population of objects with relatively large amplification 
actors decreases for larger NS masses (see Fig. 8 ). For the consid-
red EoSs, only ∼0 . 1 per cent ( ∼0 . 05 per cent ) of neutrino pulsars
emonstrate an amplification factor a ν > 10 for NS masses of 1 . 4 M 
2 M ). The expected distribution of objects over the amplification 
actor depends on the EoS insignificantly (compare different lines in 
ig. 8 ). 
In the case of strange stars, the distribution of objects o v er the
mplification factor depends on neutrino energy. About 10 per cent 
f ULX hosting strange stars can demonstrate amplification factor 
 ν > 2 and ∼0 . 01 per cent can show amplification factors a ν > 10
see Fig. 9 ).  SUMMARY  
e hav e e xplored the impact of gravitational bending on neutrino
mission in strongly magnetized NSs undergoing extreme mass 
ccretion rates, such as bright X-ray transients or ULX pulsars. NS in-
eriors in the considered class of objects are cold enough (temperature 
 10 keV) to be completely transparent to neutrino emission in keV
nd MeV energy bands (Haensel & Jerzak 1987 ). Thus, a fraction of
eutrino emission is going through an NS experiencing gravitational 
ending. Through Monte Carlo simulations in the metric generated 
y spherically symmetric and quasi-static mass distribution within an 
S, we simulated neutrino beam patterns (Figs 6 and 7 ) influenced
y neutrino gravitational bending. The gravitational bending induces 
trong anisotropy in neutrino emission within the NS reference frame, 
eading to the phenomenon of neutrino pulsars. MNRAS 538, 2396–2407 (2025) 
2404 A. A. Mushtuk o v et al. 
M
Figure 7. The angular distribution of neutrino energy flux from one of 
the poles of a strange star with M = 1 . 4 M  and M r conforming with the 
EoS proposed by Zdunik ( 2000 ). Different lines are calculated for different 
neutrino energy: 100 keV (solid red) and 1 MeV (dotted black). 
(a)
(b)
Figure 8. Fraction of neutrino pulsars of amplification factor exceeding a ν . 
Different lines show results calculated for different EoSs: APR (solid black), 
SLy4 (dotted blue), and BSk24 (dashed red). Different panels correspond to 
different masses of an NS: (a) M = 1 . 4 M  and (b) M = 2 M . One can 
see that a large mass of an NS reduces significantly a fraction of strongly 
amplified sources. 
 
c  
(  
l  
l  
p  
c  
Figure 9. Fraction of neutrino pulsars of amplification factor exceeding a ν
calculated for the case of a strange star with M = 1 . 4 M  and M r conforming 
with the EoS proposed by Zdunik ( 2000 ). Different lines are calculated for 
different neutrino energies: 100 keV (solid red) and 1 MeV (dotted black). 
The line corresponding to higher neutrino energy shows smaller fractions 
because of neutrino absorption in a star. 
o  
N  
k  
r  
fl  
e  
 
q  
a  
F  
a  
N  
a  
e  
A
T  
a  
s  
t  
H  
s  
R
D
T  
p  
t
R
A  
A
A
A
A  
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025Using calculated beam patterns, we have obtained the theoreti-
al distributions of neutrino pulsars o v er the amplification factors
equation 1 ) that show the ratio of apparent (equation 17 ) and actual
uminosity in neutrinos (see Figs 8 and 9 ). These distributions reveal
imited ranges of amplification factors. The majority of neutrino
ulsars are expected to fall within the interval a ν ∈ (0 . 5 , 10). For the
onsidered EoSs, only approximately ∼0 . 1 per cent ( ∼0 . 05 per cent )NRAS 538, 2396–2407 (2025) f neutrino pulsars exhibit an amplification factor a ν > 10 at an
S mass of 1 . 4 M  (2 M ). Thus, the expected neutrino flux from
nown pulsating ULXs and bright Be X-ray transients most likely
emain below the isotropic neutrino background even in the case of
ux amplification due to neutrino gravitational bending (see previous
stimations that neglect gravitational bending in Asthana et al. 2023 ).
In the case of strange stars, where the core of a star is composed of
uark matter, high-energy neutrinos can be subject to absorption. As
 result, the neutrino beam pattern becomes energy dependent (see
ig. 7 ), which affects the expected distribution of objects powered by
ccretion onto strange start o v er the amplification factor (see Fig. 9 ).
ote that presence of quark matter inside a star can cause some
dditional heating of the stellar interiors by the source of neutrino
mission of the surface due to the neutrino absorption in such matter.
C K N OW L E D G E M E N T S  
he authors thank Simon Portegies Zwart for the discussions. We
re grateful to an anonymous referee for their useful comments and
uggestions that helped us fix a mistake in the original version of
he manuscript and impro v e the paper. AAM thanks UKRI Stephen
awking fellowship. The work of AYP and IDM was partially
upported by the Ministry of Science and Higher Education of the
ussian Federation (agreement no. 075-15-2024-647). 
ATA  AVAI LABI LI TY  
he calculations presented in this paper were performed using a
ri v ate code developed and owned by the corresponding author. All
he data appearing in the figures are available upon request. 
EFERENCES  
kmal A. , Pandharipande V. R., Ravenhall D. G., 1998, Phys. Rev. C , 58,
1804 
lcock C. , Farhi E., Olinto A., 1986, ApJ , 310, 261 
ntoniadis J. et al., 2013, Science , 340, 448 
rzoumanian Z. et al., 2018, ApJS , 235, 37 
sthana A. , Mushtukov A. A., Dobrynina A. A., Ognev I. S., 2023, MNRAS ,
522, 3405 
Neutrino beaming in ULX pulsars 2405 
B
B
B
B
C
D
D
D
F  
F
F
H
H
H
H  
H
H
H
I
I
K  
K
K
L
L
M
M
M
M
M  
M  
M
M
M  
M  
M  
M
M  
M
M  
N  
P
P  
P  
P
P  
R
R
R
R
R  
S
S
S
S
S  
T  
T  
T
T
V
W
W
W
Y  
Z
Z
A
E  
i
T
o  
t  
s
w  
(
i  
g
2
u
w  
p  
v  
o  
u
(
L  
a  
c
T
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025achetti M. et al., 2014, Nature , 514, 202 
asko M. M. , Sunyaev R. A., 1975, A&A, 42, 311 
eloborodov A. M. , 2002, ApJ , 566, L85 
linnikov S. , Yudin A., Kramarev N., Potashov M., 2022, Particles , 5, 198 
houdhury D. et al., 2024, ApJ , 971, L20 
oroshenko V. et al., 2022, Nat. Astron. , 6, 1433 
oroshenko V. et al., 2023, A&A , 677, A57 
ouchin F. , Haensel P., 2001, A&A , 380, 151 
abrika S. N. , Atapin K. E., Vinokurov A. S., Sholukhova O. N., 2021,
Astrophys. Bull. , 76, 6 
onseca E. et al., 2021, ApJ , 915, L12 
orsblom S. V. et al., 2024, A&A , 691, A216 
aensel P. , Jerzak A. J., 1987, A&A, 179, 127 
aensel P. , Potekhin A. Y., 2004, A&A , 428, 191 
aensel P. , Zdunik J. L., Schaefer R., 1986, A&A, 160, 121 
aensel P. , Potekhin A. Y., Yakovlev D. G., 2007, Neutron Stars 1: Equation
of State and Structure. Springer, New York 
arding A. K. , Lai D., 2006, Rep. Prog. Phys. , 69, 2631 
eyl J. et al., 2024, Nat. Astron , 8, 1047 
orowitz C. J. , Piekarewicz J., Reed B., 2020, Phys. Rev. C , 102, 044321 
srael G. L. et al., 2017, Science , 355, 817 
wamoto N. , 1982, Ann. Phys. , 141, 1 
aminker A. D. , Levenfish K. P., Yakovlev D. G., Amsterdamski P., Haensel
P., 1992, Phys. Rev. D , 46, 3256 
ing A. , Lasota J.-P., Klu ´zniak W., 2017, MNRAS , 468, L59 
raus U. , 2001, ApJ , 563, 289 
asota J.-P. , King A., 2023, MNRAS , 526, 2506 
indquist R. W. , 1966, Ann. Phys. , 37, 487 
alacaria C. et al., 2023, A&A , 675, A29 
arkozov I. D. , Mushtukov A. A., 2024, MNRAS , 527, 5374 
ihalas D. , Mihalas B. W., 1985, Foundations of Radiation Hydrodynamics. 
Oxford Univ. Press, Oxford 
iller M. C. et al., 2019, ApJ , 887, L24 
isner C. W. , Thorne K. S., Wheeler J. A., 1973, Gravitation. Freeman &
Co., New York 
 ¨onkk ¨onen J. , Tsygankov S. S., Mushtukov A. A., Doroshenko V.,
Suleimanov V. F., Poutanen J., 2022, MNRAS , 515, 571 
ushtukov A. A. , Portegies Zwart S., 2023, MNRAS , 518, 5457 
ushtukov A. A. , Tsygankov S., 2022, preprint ( arXiv:2204.14185 ) 
ushtukov A. A. , Suleimanov V. F., Tsygankov S. S., Poutanen J., 2015,
MNRAS , 454, 2539 
ushtukov A. A. , Verhagen P. A., Tsygankov S. S., van der Klis M.,
Luto vino v A. A., Larchenkova T. I., 2018a, MNRAS , 474, 5425 
ushtuk ov A. A. , Tsygank ov S. S., Suleimanov V. F., Poutanen J., 2018b,
MNRAS , 476, 2867 
ushtukov A. A. , Ognev I. S., Nagirner D. I., 2019, MNRAS , 485, L131 
ushtuko v A. A. , Porte gies Zw art S., Tsygank ov S. S., Nagirner D. I.,
Poutanen J., 2021, MNRAS , 501, 2424 
ushtukov A. A. et al., 2023, MNRAS , 524, 2004 
ushtukov A. A. , Weng A., Tsygankov S. S., Mereminskiy I. A., 2024,
MNRAS , 530, 3051 
 ¨attil ¨a J. , Miller M. C., Steiner A. W., Kajava J. J. E., Suleimanov V. F.,
Poutanen J., 2017, A&A , 608, A31 
al K. , Dutt-Mazumder A. K., 2011, Phys. Rev. D , 84, 034004 
earson J. M. , Chamel N., Potekhin A. Y., Fantina A. F., Ducoin C., Dutta A.
K., Goriely S., 2018, MNRAS , 481, 2994 
ostnov K. , Shakura N., Staubert R., Kochetkova A., Klochkov D., Wilms J.,
2013, MNRAS , 435, 1147 
otekhin A. Y. , Chabrier G., 2018, A&A , 609, A74 
otekhin A. Y. , Zyuzin D. A., Yakovlev D. G., Beznogov M. V., Shibanov Y.
A., 2020, MNRAS , 496, 5052 
eig P. , 2011, Ap&SS , 332, 1 
ezzolla L. , Most E. R., Weih L. R., 2018, ApJ , 852, L25 
ichardson M. B. , van Horn H. M., Savedoff M. P., 1979, ApJS , 39, 29 
iffert H. , Meszaros P., 1988, ApJ , 325, 207 
omani R. W. , Kandel D., Filippenko A. V., Brink T. G., Zheng W., 2022,
ApJ , 934, L17 
affer A. et al., 2024, preprint ( arXiv:2412.02850 ) almi T. et al., 2024, ApJ , 974, 294 
awyer R. F. , Soni A., 1979, ApJ , 230, 859 
ch ¨afer T. , Schwenzer K., 2004, Phys. Rev. D , 70, 114037 
hapiro S. L. , Teukolsky S. A., 1983, Black Holes, White Dwarfs, and Neutron
Stars: The Physics of Compact Objects. Wiley, New York 
anashkin A. S. , Karpova A. V., Potekhin A. Y., Shibanov Y. A., Zyuzin D.
A., 2022, MNRAS , 516, 13 
sygank ov S. S. , Doroshenk o V., Luto vino v A. A., Mushtuko v A. A.,
Poutanen J., 2017, A&A , 605, A39 
sygankov S. S. et al., 2022, ApJ , 941, L14 
sygankov S. S. et al., 2023, A&A , 675, A48 
inciguerra S. et al., 2024, ApJ , 961, 62 
ang Y. M. , Frank J., 1981, A&A, 93, 255 
einberg S. , 1972, Gravitation and Cosmology: Principles and Applications 
of the General Theory of Relativity. Wiley, New York 
itten E. , 1984, Phys. Rev. D , 30, 272 
akovlev D. G. , Gnedin O. Y., Gusakov M. E., Kaminker A. D., Levenfish
K. P., Potekhin A. Y., 2005, Nucl. Phys. A , 752, 590 
dunik J. L. , 2000, A&A , 359, 311 
hang L. , Blaes O., Jiang Y.-F., 2022, MNRAS , 515, 4371 
PPENDI X  A :  GEODESI C  LI NES  
quation( 15 ) can be rewritten in a form appropriate for numerical
ntegration as follows: 
d r 
d ϕ 
= ±r 2 
(
1 
A ( r ) B( r ) b 2 −
1 
A ( r ) r 2 
)1 / 2 
. (A1) 
he numerical modelling of neutrino trajectories based on this first- 
rder differential equation requires a correct choice of the sign on
he right-hand side, as described in Appendix B . One can a v oid this
ign ambiguity by using the second-order equation 
2 A ( r) 
r 4 
d 2 r 
d ϕ 2 
+ 1 
r 4 
(
d A ( r) 
d r 
− 4 A ( r) 
r 
)(
d r 
d ϕ 
)2 
− 2 
r 3 
+ d B( r) / d r 
b 2 B 2 ( r) = 0 , (A2) 
hich is obtained by taking the deri v ati ve of both sides of equation
 15 ) o v er r . 
Based on the appendix in Beloborodov ( 2002 ), we can express b 
n terms of the trajectory variables. The tangent vector for the null
eodesic line associated with the trajectory in the metric (equation 
 ) can be written as 
 
μ = d x 
μ
d λ
, (A3) 
here μ is the index of the coordinate ( t, r, θ, ϕ) and λ is an affine
arameter. We can put θ = π/ 2 without loss of generality. Killing
ectors ∂ / ∂ t and ∂ / ∂ ϕ for equation ( 2 ) correspond to the integrals
f motion u t and u ϕ = b, respectively. If we put λ = 1 /B( r), we get
 t = −1. Then from the condition u μu μ = 0 we obtain 
 u r ) 2 = 1 
A ( r ) B( r ) −
b 2 
A ( r ) r 2 . (A4) 
et us consider the massless particle at the radius r and denote the
ngle between the particle momentum and the radial vector from the
entre of symmetry as ζ . Then 
tan ζ = 
(
u ϕ u ϕ 
u r u r 
)1 / 2 
= b 
r 
(
1 
B( r) −
b 2 
r 2 
)−1 / 2 
. (A5) 
herefore, b can be related to r and ζ as follows: 
sin ζ = b 
r 
√ 
B( r) . (A6) MNRAS 538, 2396–2407 (2025) 
2406 A. A. Mushtuk o v et al. 
M
S  
f  
n  
e
 
n  
B  
(
k
H  
m  
T  
t
i  
N  
o  
c  
a  
t
i  
ν
 
s  
a  
c  
q  
a  
i  
t  
m  
w  
t  
o  
s
 
a  
f  
c  
π  
i  
t
δ
T  
δ  
t  
l  
w
 
b
I
δ
I
d
N  
s  
e
 
a
E
T
A
T
W  
(  
g  
d  
i  
d  
v
S  
n  
 
d  
B
A
w  
t
 
p
a

 
d  
t  
e  
t  
e  
s  
f
 
p  
p
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025ince the vector fields ∂ / ∂ t and ∂ / ∂ ϕ are the Killing fields both
or the Schwarzschild metric and equation ( 2 ), the value of b does
ot change if the particle crosses the NS surface. Thus we arrive at
quation ( 16 ). 
Let us consider opacity transformation in general relativity. The
eutrino transport process can be described by the relativistic
oltzmann equation for massless particles that can be written as
Lindquist 1966 ) 
 
α ∂ I 
∂ x α
−  αβγ k βk γ
∂ I 
∂ k α
= J − κI. (A7) 
ere,  αβγ are the Christoffel symbols, k α is the particle four-
omentum, and I = I ν( ) /ν3 is the invariant specific intensity.
he ordinary specific intensity I ν( ) is usually defined in relation
o the radiative transfer (e.g. Mihalas & Mihalas 1985 ), where ν
s the photon frequency and  is the photon propagation direction.
ote that in general relativity the photon frequency can be defined
nly in the local rest frame associated with an observer. In our
ase, the frequency is ν = k ˆ 0 /h , where h is the Planck constant
nd k ˆ 0 is the neutrino energy measured in the reference frame of
he observer whose ( r, θ, ϕ) coordinates do not change. Since ∂ / ∂ t 
s the Killing vector for the spherical static metric (equation 2 ),√ 
B( r) = constant . 
Furthermore, J = j ν/ν2 in equation ( A7 ) is an invariant emis-
ivity, j ν being an ordinary emissivity; κ = ναν is an invariant
bsorption coefficient, αν ∝ 1 /λν , being an ordinary absorption
oefficient, and λν is a mean free path at the frequency ν. The
uantities j ν , αν , and λν are defined in the same reference frame
s the frequency ν. For an accurate calculation of neutrino transfer
n NSs, it is necessary to take into account the transformation of
he mean free path along the geodesic line due to the change of the
etric coefficient B( r). Since the typical mean free path of a neutrino
ith energy of a few hundred keV is very large in comparison with
he typical NS radius (see Section 2.2.1 ), we can neglect neutrino
pacities in the NSs. In the quark stars, we can neglect neutrino
cattering, but should take into account neutrino absorption. 
In our numerical model, we trace the motion of each neutrino
s it propagates through an NS. Let us consider a neutrino moving
rom ( t, r, π/ 2 , ϕ) to ( t + d t, r + d r, π/ 2 , ϕ + d ϕ) in Schwarzschild
oordinates ( t, r, θ, ϕ) (without loss of generality, we assume θ =
/ 2 = const ). The spatial displacement vector δl = (d r, 0 , d ϕ) lies
n the tangent space at the point ( t, r, π/ 2 , ϕ) and is represented in
he coordinate basis as 
l = d r ∂ 
∂ r 
+ d ϕ ∂ 
∂ ϕ 
. (A8) 
he optical depth d τ associated with this infinitesimal displacement
l is d τ = d s /λ, where λ is the neutrino mean free path and d s is
he length of the spatial motion. Both quantities are e v aluated in the
ocal Minkowski frame corresponding to the element of matter with
hich the neutrino interacts. 
We assume the NS matter is at rest; therefore, the local orthonormal
asis is 
e (0) = 1 √ 
B( r) 
∂ 
∂ t 
, e (1) = 1 √ 
A ( r) 
∂ 
∂ r 
, 
e (2) = 1 √ 
r 
∂ 
∂ θ
, e (3) = 1 √ 
r sin 2 ( θ ) 
∂ 
∂ ϕ 
. (A9) 
n this basis, the displacement vector is 
l = 
√ 
A ( r) d r e (1) + r d ϕ e (3) . NRAS 538, 2396–2407 (2025) ts length is 
 s = 
√ 
A ( r )(d r ) 2 + r 2 (d ϕ) 2 . (A10) 
ote that d s is numerically identical to the length computed in the
patial part of the spherically symmetric metric (equation 2 ), because
 ( 0 ) is parallel to ∂ / ∂ t , as the NS matter is at rest. 
The mean free path depends on the neutrino energy, which changes
long the geodesic due to gravitational redshift: 
 ν ∝ 1 / 
√ 
B( r) . (A11) 
his effect is fully accounted for in our Monte Carlo modelling. 
PPENDI X  B:  SI MULATI ONS  O F  N E U T R I N O  
R A J E C TO R I E S  
e calculate neutrino trajectories, described by differential equation
 15 ), where the mass distribution is spherically symmetric and
iven by M r . The impact factor can be calculated from the initial
irection of particle motion according to equation ( 16 ). A trajectory
s determined by the initial coordinates of a particle r 0 and initial
irection of its motion, which is given by the unit vector of particle
elocity: 
e v0 = v 0 | v 0 | . (B1) 
imulating a trajectory, we choose a spatial separation between the
earest two points of approximate trajectory s and follow the steps:
(i) Using the starting point of particle trajectory r 0 and the
irection of its initial velocity given by the unit vector e v0 (equation
1 ), we calculate the second point of approximate trajectory: 
r 1 = r 0 + e v0 s. (B2) 
t this step i = 1. 
(ii) Then we get the angle between positions r i and r i−1 : 
cos ϕ i ,i −1 = ( r i , r i−1 ) | r i | | r i−1 | , (B3) 
here ( r i , r i−1 ) = 
∑ 3 
j= 1 r 
( j ) 
i r 
( j ) 
i−1 denotes the scalar productions of
wo vectors and r ( j ) i is j th Cartesian coordinate of vector r i . 
(iii) We get direction towards the ( i + 1) th point of approximate
article trajectory: 
e r,i+ 1 = r i + e v,i−1 s | r i + e v,i−1 s| , (B4) 
nd the angle between e r,i+ 1 and r i : 
ϕ ∗i+ 1 ,i = 
( e r,i+ 1 , r i ) 
| r i | . (B5) 
(iv) Using the second-order Runge–Kutta method, applied to the
ifferential equation ( A1 ), we compute the radial distance r i+ 1 at
he next step of the simulation. The sign on the right-hand side of
quation ( A1 ) is determined based on whether the particle is moving
ow ards or aw ay from the centre of the star. If the right-hand side of
quation ( A1 ) becomes zero at any step, the sign changes in the next
tep. This corresponds to the particle reaching its minimum distance
rom the centre for a given impact parameter b. 
(v) We get an estimation of the radial distance towards a new
oint of particle trajectory r i+ 1  0 . 5 R Sch ( M i ) /u i+ 1 and calculate its
osition: 
r ∗i+ 1 = r i+ 1 e r,i+ 1 . (B6) 
Neutrino beaming in ULX pulsars 2407 
o  
o
T  
i
e
s  
o 2  
o  
s
f  
o
 
b  
e  
i
T
©
P
(
D
ow
nloade(vi) Because we want to get trajectory approximated by segments 
f a fixed length s, we recalculate the position of the latest point
f neutrino trajectory as 
r i+ 1 = r i + r 
∗
i+ 1 − r i 
| r ∗i+ 1 − r i | 
s. (B7) 
he unit vector of neutrino velocity at the latest segment of trajectory
s given by 
 v,i = r i+ 1 − r i | r i+ 1 − r i | . (B8) 
(vii) We stop trajectory simulation if the particle experiences 
cattering at a given coordinate or if it is far from the central compact2025 The Author(s). 
ublished by Oxford University Press on behalf of Royal Astronomical Society. This is an Open
 https://cr eativecommons.or g/licenses/by/4.0/), which permits unrestricted reuse, distribution, and repbject: | r i+ 1 | > 5 × 10 R Sch . In this case, we have a final direction
f particle motion given by equation ( B8 ). Otherwise, we return to
tep (iii) and continue the simulation of the trajectory. 
To control the accuracy of trajectory calculations, we perform it 
or smaller spatial step s 1 = 0 . 5 s. In the case of similar results
f the simulation, we stop the impro v ement of accuracy. 
The results of the performed algorithm outside an NS were verified
y comparison of its results with the results of algorithms applied
arlier by Mushtukov et al. ( 2018a , 2024 ) to trace photon trajectories
n XRPs. 
his paper has been typeset from a T E X/L A T E X file prepared by the author. MNRAS 538, 2396–2407 (2025) 
 Access article distributed under the terms of the Creative Commons Attribution License 
roduction in any medium, provided the original work is properly cited. 
d from
 https://academ
ic.oup.com
/m
nras/article/538/4/2396/8082115 by Turun Yliopiston Kirjasto user on 09 June 2025