2023 Vol., X No. XX, 000–000

11institutetext:  1 Purple Mountain Observatory, Chinese Academy of Sciences,
No. 10 Yuan Hua Road, Nanjing 210023, China;
2 Department of Information Engineering, Wuhan Institute of City, Wuhan, Hubei
430083, China;
3 Shanghai Astronomical Observatory, Nandan Road 80, Shanghai 200030, China;
4 National Basic Science Data Center, Building No.2, 4, Zhongguancun South 4th
Street, Haidian District, Beijing 190, China;
5 Zhejiang University-Purple Mountain Observatory Joint Research Center for
Astronomy, Zhejiang University, Hangzhou 327, China;
6 South-Western Institute for Astronomy Research, Yunnan University, Kunming,
650500, China;
7 The Shanghai Key Lab for Astrophysics, Shanghai Normal University, Shanghai,
200234, China;
8 Chinese Academy of Sciences South America Center for Astronomy, National
Astronomical Observatories, CAS, Beijing 100101, China;
9 CAS Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, China;
10 Key Laboratory of Radio Astronomy and Technology, Chinese Academy of
Sciences, A20 Datun Road, Chaoyang District, Beijing, 100101, China;
11 University of Chinese Academy of Sciences, Beijing 100049, China

Received 2023 October XX; accepted XXXX July XX

# How to coadd images: II. Anti-aliasing and PSF deconvolution

Lei Wang E-mail<EMAIL_ADDRESS>Huanyuan Shan 3310101111 Lin Nie 2244
Dezi Liu 66771111 Zhaojun Yan 33 Guoliang Li 114455 Cheng Cheng 8899 Yushan
Xie 33 Han Qu 114455 Wenwen Zheng 114455 Xi Kang 55

###### Abstract

We have developed a novel method for co-adding multiple under-sampled images
that combines the iteratively reweighted least squares and divide-and-conquer
algorithms. Our approach not only allows for the anti-aliasing of the images
but also enables PSF deconvolution, resulting in enhanced restoration of
extended sources, the highest PSNR, and reduced ringing artefacts. To test our
method, we conducted numerical simulations that replicated observation runs of
the CSST/VST telescope and compared our results to those obtained using
previous algorithms. The simulation showed that our method outperforms
previous approaches in several ways, such as restoring the profile of extended
sources and minimizing ringing artefacts. Additionally, because our method
relies on the inherent advantages of least squares fitting, it is more
versatile and does not depend on the local uniformity hypothesis for the PSF.
However, the new method consumes much more computation than the other
approaches.

###### keywords:

Methods: analytical – Techniques: image processing – Gravitational lensing:
weak

##  1 Introduction

As technology has advanced, computer software has become increasingly
important in the field of image processing. To simplify data processing, an
Image Acquisition System (IAS) is often utilized for image digitization.
Essentially, an IAS acts as a digitizer that converts continuous signals into
digital ones by recording data with detectors, such as pixels, which appear in
a mosaic-like pattern known as pixelation. However, due to the diffraction
limit of the optics equipment, the images captured by the IAS are band-
limited. This means that the signal recorded by the image has a maximum
spatial frequency or resolution, as described in Fruchter (2011). A band-
limited signal can be fully reconstructed by high-sampling imaging, according
to the Nyquist sampling theorem. Nevertheless, many astronomical images do not
meet the Nyquist sampling criterion due to technical or economic limitations,
resulting in under-sampled images(Fruchter & Hook 2002).

Under-sampled images suffer from an aliasing effect, which blurs all signals
below a sampling interval. To achieve a resolving power that approaches the
diffraction limit, the sampling rate needs to be increased. This can be
accomplished by combining multiple exposures using various coaddition methods,
such as shift-and-add (Bates & Cady 1980; Farsiu et al. 2004b), Drizzle
(Fruchter & Hook 2002), Super-Drizzle (Takeda et al. 2006), IMCOM (Rowe et al.
2011), iDrizzle (Fruchter 2011), SPRITE (Ngolè Mboula, F. M. et al. 2015), and
fiDrizzle (Wang & Li 2017), as well as iterative back-projection (IBP, Irani &
Peleg 1993; Symons et al. 2021). Drizzle has become the standard for combining
images taken by the Hubble Space Telescope (HST) and the James Webb Space
Telescope (JWST). Some Drizzle-based methods are widely used to restore fine
details of under-sampled multi-exposures and fuse images from different
equipment. SPRITE and IBP are developed to reconstruct the Point Spread
Function (PSF) using stars in even one exposure.

Due to the diffraction of light, the image of a point-like source in the image
plane is not actually a point. Instead, it appears as an extended object,
known as the PSF effect. This degradation of the observed image is caused by
the convolved PSF, but PSF deconvolution technology can improve resolution by
compensating for this numerically. Along with resolution improvement, PSF
deconvolution also enhances contrast and reduces noise(Sage et al. 2017).
There are various algorithms for PSF deconvolution(Starck et al. 2002),
including the Fourier-quotient method, CLEAN method(Högbom & Cornwell 1974),
Bayesian approach[including Landweber method(Landweber 1951), Richardson-Lucy
algorithm(Richardson 1972; Lucy 1974; Shepp & Vardi 1982) ], wavelet-based
deconvolution, and super-resolution techniques(Gerchberg 1974; Hunt 1994;
Lauer 1999; Elad & Feuer 1999; Capel & Zisserman 2003; Park et al. 2003;
Farsiu et al. 2004a; van Ouwerkerk 2006; Tian & Ma 2011; Nasrollahi & Moeslund
2014; Yue et al. 2016; Symons et al. 2021). The PSF deconvolution has numerous
benefits, such as its applicability to even the simplest optical setup,
reduction of financial costs, and streamlining of the acquisition pipeline.
However, when noise contaminates the PSF-convolved image, PSF deconvolution
becomes an ill-posed problem(Starck et al. 2002; Sage et al. 2017).
Regularized methods are often employed to generate an approximation(Takeda et
al. 2007; Ng et al. 2007; Takeda et al. 2009; Yuan et al. 2010; Babacan et al.
2011; Su et al. 2012; Liu & Sun 2014; Zhang et al. 2012), which is effective
and flexible in reducing noise amplification, ringing effect, and flux
divergence.

However, most previous works aim to either the PSF deconvolution for single
exposure (e.g. CLEAN, Landweber, Richardson-Lucy, maximum entropy algorithm,
etc.) or to multiple exposures coaddition but without PSF deconvolution (e.g.
shift-and-add and Drizzle based methods). Only a few works mention the under-
sampled multi-exposures coaddition (MEC) with PSF deconvolution, e.g. Stark-
Pantin: Eq. 53 in Starck et al. (2002) and UPDC: Eq. 10 in Wang et al. (2022).
In this paper, we systematically study how to achieve superresolution in the
iterative MEC by anti-aliasing and PSF deconvolution coaddition (AAPD). We
test several recovery methods with numerical simulations.

##  2 Imaging model and the least squares

The observation image records not only the shape or flux from the objects of
study but also a set of combined observational effects, e.g. blurring effects
(PSF, seeing, vignetting, etc.), sampling effect (pixelation, image field
distortion), noise effect (equipment noise, environment noise), etc.
Therefore, a reasonable imaging model should take those observational effects
into account.

Let $\mathbb{B}$ be an image degrading operator (i.e. blurring) which
represents a combination of a series of image operations e.g. PSF convolution,
pixelation, etc. The formation of an under-sampled observation image
$\mathcal{G}$ with the size of $U\times V$ pixels can be formally expressed as

$\mathcal{G}=\mathbb{N}\left\\{\mathbb{B}\\{\mathcal{O}\\}\right\\},$ (1)

where $\mathcal{O}$ is the intrinsic, continuous surface brightness of the
sky. Assuming that $\mathcal{O}$ contains a fine grid of size $M\times N$, we
can represent the downgraded effects of the image, except for the noise, using
$\mathbb{B}\\{\mathcal{O}\\}$. The noise contamination on the downgraded image
is described by $\mathbb{N}$. To recover the original image $\mathcal{O}$ from
$\mathcal{G}$, we need to know the combined effects of
$\mathbb{N}\\{{\mathbb{B}}\\{\cdot\\}\\}$ in Eq.(1) in advance. In this study,
we focus on the recovery method and assume that the PSFs and the position
shift on all dithered exposures are well measured beforehand. Table 1 provides
a list of important symbols and their representations.

Supposing that we have $L$ blurred (PSF convolved and/or under-sampled)
exposures $\\{\mathcal{G}_{1},\mathcal{G}_{2}...\mathcal{G}_{L}\\}$111Frame
$\mathcal{G}_{k}$ ($k=1...L$) is obtained from the $k$-th observation with
$U\times V$ pixels for the same object and the known blurring effects
$\mathbb{B}\\{\cdot\\}$, can we mock another set of frames
$\\{\mathcal{F}_{1},\mathcal{F}_{2}...\mathcal{F}_{L}\\}$ to approach the
observations, then constrain the original image? To answer this question, we
try to use the least squares to calculate a $\chi^{2}$ for all pixels between
the observation and mock samples. Theoretically, the $\chi^{2}$ can be
represented as

$\displaystyle\chi^{2}=\sum^{L,U,V}_{k,i,j}\frac{(\mathcal{F}_{k,i,j}-\mathcal{G}_{k,i,j})^{2}}{\theta_{k,i,j}^{2}},$
(2)

where $\mathcal{F}={\mathbb{B}}\\{\mathcal{T}\\}$ and $\mathcal{T}$ is a
target image with size of $M\times N$ fine grids, i.e. the estimation of the
original image $\mathcal{O}$. While $\theta_{k,i,j}$ is the noise at the pixel
of the $i$-th row and $j$-th column of the $k$-th observational frame. Note
that the summation of the squares is over $L\times U\times V$ pixels. Since
the blurring effects $\mathbb{B}\\{\cdot\\}$ are known, minimizing the
$\chi^{2}$ results in an estimation of the target image $\mathcal{T}$. By
differentiating $\chi^{2}$ with respect to $\mathcal{T}$ and setting the
derivative to zero, we have

$\displaystyle\frac{\partial\chi^{2}}{\partial\mathcal{T}}=0.$ (3)

It is very difficult to solve the equation Eq. 3 directly, especially for huge
amounts of pixels. Actually, Eq. 3 is a system of linear equations with
$M\times N$ unknowns and $L\times U\times V$ conditions. As a kind of fallback
solution, we employ a divide-and-conquer algorithm to solve Eq. 3 on each
target grid ($m,n$) at a time. The $Drizzled$ grids are used as the initial
values. The target grid ($m,n$) will be renewed after solving Eq. 3. The
program checks if the output meets the preset threshold222In this work, we
suggest using the $\chi^{2}$ as the threshold. when all grids are updated. If
not, the updated grids will take the place of the initial values in Eq. 3.
Then it repeats the above steps until the updated grids meet the threshold.
Therefore, Eq. 3 is solved in an iterative approach like the following:

$\displaystyle\sum^{OS}_{k,i,j}\frac{1}{\theta_{k,i,j}^{2}}\sum^{TS}_{p,q}\mathcal{W}_{k,i,j,p,q}\big{(}\sum^{TS}_{s,t}\mathcal{P}_{k,p,q,s,t}\mathcal{T}^{(r)}_{s,t}+(\mathcal{T}^{(r+1)}_{m,n}-\mathcal{T}^{(r)}_{m,n})\mathcal{P}_{k,p,q,m,n}\big{)}\sum^{TS}_{p,q}\mathcal{W}_{k,i,j,p,q}\mathcal{P}_{k,p,q,m,n}$
(4)
$\displaystyle=\sum^{OS}_{k,i,j}\frac{\mathcal{G}_{k,i,j}}{\theta_{k,i,j}^{2}}\sum^{TS}_{p,q}\mathcal{W}_{k,i,j,p,q}\mathcal{P}_{k,p,q,m,n}$

where $\mathcal{P}_{k}$ is a fine PSF measured from the $k$-th frame. The PSF
must have the same sampling rate as the target image.
$\mathcal{W}_{k,i,j,p,q}$ is a flux allocation weight that stands for how many
areas of an observation pixel ($k,i,j$) are overlapped by the target pixel
($p,q$), shown as Figure 1. In this work, the overlapped polygon of two grids
is clipped using a polygon clipping algorithm
Clipper2333https://github.com/AngusJohnson/Clipper2. The area of the clipped
polygon is calculated using the Green formula. Indexes $r$ and $r+1$ represent
the sequence number of the iteration step. $OS$ is a space that includes all
observation pixels that are overlapped with PSF $\mathcal{P}_{k}$ which is
centred at the target pixel ($m,n$), while $TS$ includes all target pixels
that are overlapped with the $OS$ observation pixels.

Figure 1: Resampling mechanism. The target grid $\mathcal{T}^{(r)}_{m,n}$ is
shown as the fine (in black), the observation pixels $\mathcal{G}_{k,i,j}$ are
shown as coarse (in red). The flux weight $\mathcal{W}_{k,i,j,p,q}$ (in pink)
is determined by the clipped area between the observation pixels and the
target grid. Here the $TS$ space ($p,q$) includes all target cells that
overlap with the observation pixel $\mathcal{G}_{k,i,j}$ (blue region).

Let
$\Psi_{k,i,j,m,n}=\sum^{TS}_{p,q}\mathcal{W}_{k,i,j,p,q}\mathcal{P}_{k,p,q,m,n}$
and
$\Phi_{k,i,j,m,n}^{(r)}=\sum^{TS}_{p,q}\mathcal{W}_{k,i,j,p,q}\sum^{TS}_{s,t}\mathcal{P}_{k,p,q,s,t}\mathcal{T}^{(r)}_{s,t}$,
the Eq. 2 can be rewritten as

$\displaystyle\mathcal{T}^{(r+1)}_{m,n}=\mathcal{T}^{(r)}_{m,n}+\frac{\sum^{OS}_{k,i,j}\frac{\Psi_{k,i,j,m,n}}{\theta_{k,i,j}^{2}}(\mathcal{G}_{k,i,j}-\Phi^{(r)}_{k,i,j,m,n})}{\sum^{OS}_{k,i,j}\frac{\Psi_{k,i,j,m,n}^{2}}{\theta_{k,i,j}^{2}}}$
(5)

which is an iterative solution to the least squares. $\Phi_{k,i,j,m,n}^{(r)}$
is the mimic flux of pixels which are overlapped by the target grid
$\mathcal{T}^{(r)}_{m,n}$. In order to eliminate contaminations from cosmic
rays or abnormal pixels, the weight $\theta_{k,i,j}$ is dynamically adjusted
with the iterations, which results in robust regression of the algorithm. This
is actually an application of the Iteratively Reweighted Least Squares (IRLS)
method in image stacking. We call this divide-and-conquer IRLS method i.e. Eq.
5 as Dirles. The initial value of the iteration $\mathcal{T}^{(0)}$ is set as
the $Drizzled$ real exposures to automatically introduce telescope effects
such as vignetting, instrument noise, and saturation overflow. Therefore,
these effects do not need to be simulated anymore in the Dirles.

In literature, there are two algorithms besides the Dirles which can achieve
the anti-aliasing and PSF deconvolution coaddition simultaneously. These
algorithms are called Starck-Pantin (Starck et al. 2002) and UPDC (Wang et al.
2022). Both Starck-Pantin and UPDC are forward modelling methods that infer
the original image by maximizing the likelihood between observations and
mocks. The algorithms like Richardson-Lucy, Starck-Pantin, and UPDC use FFT to
convolve the PSF. However, the Dirles algorithm convolves the PSF in the real
space, making it more suitable for situations where the PSF changes
dramatically. Although it consumes much more computation than the previous
algorithms, it is a better choice for such situations.

Drizzling the multi-exposures of an ideal point source can form a blurring
effect in the target image (grid), which originates from the pixelation. It is
called pixelation blur to distinguish it from the ordinary PSF(see Wang et al.
2022 for more details). The pixelation blur is not homogeneous, even not
continuously changing in the whole target image. Because different point
sources have different positions therefore different cases of pixels
coaddition. Due to the heterogeneity of pixelation blur, theoretically, it is
not a good choice to deconvolve the PSF by using ordinary methods e.g.
Landweber or Richardson-Lucy algorithm on the $Drizzled$ image.

Figure 2: Cosmic-ray removal. For one of the exposures, there are 4 circles in
the top panels and bottom left. The bottom right panel shows all the
mismatched cases (including 50 exposures). Circles in green stand for the
misidentification cases (They are not real cosmic rays), while blue for the
omissions.

##  3 Image completion: cosmic rays and bad pixels replacement

Instead of using the $Drizzle$ method for cosmic-ray removal(Fruchter & Hook
2002), we apply an improved statistical algorithm that identifies and replaces
abnormal regions with reasonable values. The new algorithm can significantly
reduce the mismatched cosmic rays ($\sim 50\%$) in our CSST-MCI mock test,
e.g. Figure 2. The process of image completion involves the following steps.

  1. 1.

To start, let’s create a target grid, denoted by $\mathcal{T}_{m,n}$, with a
1:1 sampling rate. This grid will be utilized to resample all $L$ exposures
$\mathcal{G}_{k,i,j}$ ($k=1...L$) based on the $WCS$ parameters specified in
their headers.

  2. 2.

For each target pixel located at position $(m,n)$, we collect a statistical
sample of overlapping exposure pixels, denoted by $\mathbb{G}_{0}(m,n)$. We
gather a total of $M\times N$ samples. These statistical samples are then
utilized to calculate the median value, providing the first estimation of an
image that is free of cosmic rays or bad pixels, known as $\mathcal{T}_{\rm
mid}(m,n)$.

  3. 3.

The mean, $\mu_{0}(m,n)$, and standard deviation, $\sigma_{0}(m,n)$, are also
calculated for each sample $\mathbb{G}_{0}(m,n)$. We use these values to
identify cosmic-ray candidates by checking if the exposure pixels with a flux
of $\mathcal{G}_{k,i,j}>\mathcal{T}_{\rm
mid}(m,n)+\nu_{0}\times\sigma_{0}(m,n)$ meet the preliminary voting condition.
Here, $\nu_{0}$ is an artificial threshold set to 5 in our work. If a region
occasionally has multiple cosmic rays, the mean $\mu_{0}(m,n)$ can be
overestimated. Therefore, we replace the mean $\mu_{0}(m,n)$ with the median
$\mathcal{T}_{\rm mid}(m,n)$ in the preliminary voting condition.

  4. 4.

Typically, an exposure pixel, $\mathcal{G}_{k,i,j}$, is covered by 1 to 5
neighbouring target pixels that act as ”voters” in determining whether a
cosmic-ray candidate is present. If a cosmic-ray candidate receives all the
votes of the voters in the district, it is flagged as a cosmic ray. We remove
the cosmic rays from the sample $\mathbb{G}_{0}(m,n)$ to obtain a relatively
pure sample, denoted by $\mathbb{G}_{1}(m,n)$.

  5. 5.

Once again, we calculate the mean, $\mu_{1}(m,n)$, and standard deviation,
$\sigma_{1}(m,n)$, for each sample $\mathbb{G}_{1}(m,n)$. We then perform a
series of similar operations as described above, but with a refreshed
condition $\mathcal{G}_{k,i,j}>\mu_{1}(m,n)+\nu_{1}\times\sigma_{1}(m,n)$,
where $\nu_{1}$ is set to 5. This allows us to remove cosmic rays from sample
$\mathbb{G}_{1}(m,n)$, resulting in a reduced sample called
$\mathbb{G}_{2}(m,n)$. Please note that some exposure pixels with flux
$\mathcal{G}_{k,i,j}>\mu_{1}(m,n)+\nu_{0}\times\sigma_{1}(m,n)$ are removed
for their cosmic-ray identification, but those with flux
$\mathcal{G}_{k,i,j}<\mu_{1}(m,n)-\nu_{0}\times\sigma_{1}(m,n)$ are not. This
means that using the mean $\mu_{2}(m,n)$ and standard deviation
$\sigma_{2}(m,n)$ taken from sample $\mathbb{G}_{2}(m,n)$ to generate a
Gaussian distribution is biased. To avoid that issue, we fit the sample
$\mathbb{G}_{2}(m,n)$ to a Gaussian profile using the non-linear least squares
fitting code $\it MPFIT$444https://pages.physics.wisc.edu/
craigm/idl/cmpfit.html.

  6. 6.

Each bad pixel or cosmic ray exposure pixel will be replaced with a random
value generated by its corresponding Gaussian profile.

The above method is not valid for pixels that are overlapped by only one or
two exposures due to the lack of a statistical sample. These cases usually
occur at the marginal area of the target image $\mathcal{T}_{m,n}$, and can be
ignored. Figure 2 shows a case for cosmic-ray removal. The circles and
ellipticals denote the mismatched cosmic rays. The bottom right panel shows
the total number of mismatched cases for 50 exposures. The blue circles and
ellipticals represent cosmic rays that were omitted due to their low flux
blending with the background, while the green circles indicate cases of
mistaken removal caused by the threshold set-up, i.e. a small probability
event below the threshold.

Figure 3: CSST-MCI mock and multi-exposure coaddtion. The bottom left panel
shows the ground-truth image as the mock input, which is extracted from the
HST observation with at least 15 substructures on the central galaxy. One of
the $g$-band exposures with 2 obvious substructures is illustrated on the top
left. There are about 7 substructures on the top right panel (the Drizzled
image for 50 exposures). All the 15 substructures on the central galaxy can be
found on the last panel, which is coadded by this work.

##  4 Simulation Tests based on the CSST-MCI mock pipeline

By employing the CSST-MCI mock pipeline, we take one HST Wide Field Camera 3
(WFC3) stamp from a real galaxy database of $Galsim$(Rowe et al. 2015) as our
mock input (or original image, ground truth), to generate mock exposures with
random shifts and rotations for $g$-band. In the CSST-MCI mock, four blurring
effects are taken into account: PSF convolution, cosmic ray, down-sampling,
and Gaussian/Poisson noise. Following Wang et al. (2022), a non-strict
positivity constraint is adopted:

$\mathscr{P}_{\mathbb{R}^{+}}\\{\mathcal{T}^{(i+1)}\\}=\begin{cases}\mathcal{T}^{(i+1)}&\mathcal{T}^{(i+1)}\geq
0\\\ \mathcal{T}^{(i)}&\mathcal{T}^{(i+1)}<0\end{cases}$ (6)

where $\mathscr{P}_{\mathbb{R}^{+}}\\{\mathcal{T}^{(i+1)}\\}$ is a component-
wise projection of $\mathcal{T}^{(i+1)}$ onto the set $\mathbb{R}^{+}$ (not
strictly). For the convenience of comparison, we assumed that a PSF is uniform
in a local region. The region size is determined by the spatial growth rate of
the PSF. It is a local uniformity hypothesis for the PSFs.

###  4.1 Mock-I: Faint sources detection

This case centres on the application of mock and multi-exposure coaddition to
CSST-MCI, i.e. Figure 3. The image in the lower left panel serves as the
ground-truth input, extracted from an HST observation and showcasing no less
than 15 substructures555The substructures of galaxies are detected by the
software SourceXtractor++
(https://sourcextractorplusplus.readthedocs.io/en/latest/Introduction.html).
It may be point-like sources or spiral arms in the figure. on the central
galaxy. The upper left panel features one of the $g$-band exposures
($512\times 512$ pixels, with a down-sampling factor $\beta=2$, convolved PSF
size $128\times 128$ pixels), which displays two easily identifiable
substructures. The $Drizzled$ image for 50 exposures is displayed in the upper
right panel, revealing roughly seven substructures. Lastly, the final panel
showcases the coadded image generated by this work, where all 15 substructures
on the central galaxy are clearly visible. There are 7 stars in the field of
view, corresponding to 7 pixels with high flux in the ground-truth input. The
stars are visible in the rest panels but with different full half-maximum
widths (FHWM).

###  4.2 Mock-II: Strong lensing mock and images coaddition

Figure 4 shows a strong gravitational lensing mock and recoveries from four
kinds of coaddition methods. The panel in the bottom left displays the ground-
truth image generated by a lensing model and used as the mock input. Notably,
the strong lensed giant arc within the dashed annulus is known as the Einstein
ring. Moving to the top middle panel, we see an example of one of the CSST-MCI
mocked exposures in the $r$-band with a down-sampling factor $\beta=4$.
Finally, the remaining panels showcase four distinct coaddition methods that
were used to recover images from the exposures ($128\times 128$ pixels,
convolved PSF size $128\times 128$ pixels). These methods are referred to as
Drizzle (top right), UPDC (bottom left), Starck-Pantin (bottom middle), and
Dirles (bottom right). Visually, UPDC and Dirles manifest better results than
the Drizzle or the Starck-Pantin in the Einstein ring recovery from 100 mocked
exposures.

Figure 4: A strong gravitational lensing mock. The bottom left panel shows the
ground-truth image as the mock input, which is generated by a lensing model.
One of the CSST-MCI mocked exposures in $r$-band is illustrated on the top
middle panel (with a down-sampling factor $\beta=4$). Four kinds of coaddition
methods recovered images from exposures are shown on the rest panels: Drizzle
(top right), UPDC (bottom left), Starck-Pantin (bottom middle) and Dirles
(bottom right). Figure 5: A comparison of the PSNR from different coaddition
methods: Drizzle ($PSNR=22.16$), UPDC (in Red), Richardson-Lucy (in Green),
Starck-Pantin (in Blue) and Dirles (in Black). The curves end at the optimal
iteration with the highest PSNR.

Quantitatively, we measure the PSNR within the region of the dashed annulus in
Figure 4 for each iteration of different coaddition methods. The result is
shown in Figure 5. The methods tested are Drizzle (with a PSNR of 22.16), UPDC
(in Red), Richardson-Lucy (in Green), Starck-Pantin (in Blue) and Dirles (in
Black). The curves end at the optimal iteration with the highest PSNR. The
Dirles achieves the highest PSNR among all the recoveries with less than 40
iterations. Richardson-Lucy method consumes the most iterations but the least
computation(see Table 2). Note that the same PSFs are used in the PSF
deconvolutions as those for the mock input.

###  4.3 Mock-III: VST mock and image restoration

To test the recovery methods, we have mocked the 100 frames $r$-band $VOICE-
CDFS-1$ multi-exposures data from the OmegaCAM of the Very Large Telescope
Survey Telescope in the European Southern Observatory. To reduce computation
time, we extracted a square region that is centred at
[$R.A.=53.1511^{\circ},DEC.=-27.7175^{\circ}$] and has a size of $101\times
101$ observation pixels. There is a high-resolution image
(goodss_3dhst_F606W_sci.fits) from the HST in this region. The high-resolution
HST image has been binned $1.77858\times 1.77858$ times along the two
dimensions.

In our simulation process, we use a group of Moffat profiles to replicate
extended sources such as galaxies. To mimic the central point source (the
star) on the target grid, we adopt a single pixel. By adjusting the Moffat
profiles’ parameters and the central point’s brightness, and convolving a PSF
from the binned HST image on the same sky field, we can fit a target image to
the binned HST image. The target image is shown on the top left panel of Fig.
6, and we consider it to be true.

To produce the mock samples, we first combine the fine PSFs and true PSFs
measured in the 100 VOICE frames. The fine PSFs, which have a dimension of
$59\times 59$ pixels, are generated by a Principal Component Analysis (PCA)
based method developed by our team’s cooperators Nie et al. (2021a, b). Next,
we shift the image using the same position shifts as the VOICE frames and
down-sample the true image to the same coarse grids as in the VOICE frames.
Finally, we add the Poisson and Gaussian noise from each VOICE frame to the
corresponding mock sample, as shown in the middle panel of the top array of
Fig. 6.

To make the results more realistic, we create sets of realizations by
modifying the random seed in the Poisson and Gaussian processes. For each set
of realizations, we reconstruct the mock samples using the methods mentioned
in the text, namely Drizzle, UPDC, Starck-Pantin, and Dirles, up to the
optimal number of iterations666The iteration has the highest PSNR. in the
$2\times 2$ finer target grid.

According to the results depicted in Figure 6, the Dirles approach appears to
yield superior reconstructions in comparison to earlier methods. This can be
attributed to its ability to effectively recover object shapes, minimize the
ringing effect, and suppress background noise. It should be noted, however,
that Dirles may disregard certain low SNR regions, such as the tails of
galaxies. To circumvent this potential concern, decreasing the number of
iterations and gradually balancing noise reduction with signal enhancement is
advisable.

Figure 6: Recoveries from the VST mock. The top left panel displays the ground
truth, which is a model fitting to the binned HST observation. One of the VST
mock exposures is shown on the top middle panel. The top right panel displays
the image recovered by using Drizzle method. The three panels at the bottom,
from left to right, represent the recoveries by UPDC, Starck-Pantin and Dirles
methods respectively.

Weak gravitational lensing is primarily concerned with measuring the shear.
This involves measuring the deformation of the background galaxies compared to
the randomly aligned ones. By doing so, it is possible to constrain the
properties of the foreground lensing objects. Therefore, the ellipticity of
galaxies is the signal that researchers are interested in extracting.
Following Hirata & Seljak (2003), the ellipticity of an object is defined as

$\displaystyle e_{+}$ $\displaystyle=(M_{xx}-M_{yy})/(M_{xx}+M_{yy})$
$\displaystyle e_{\times}$ $\displaystyle=2M_{xy}/(M_{xx}+M_{yy})$ (7)

where $M_{ij}$ represents the moments. the spin-2 tensor
$\textbf{e}=(e_{+},e_{\times})$ is the so-called ellipticity tensor. To
address the problem of divergence, we apply a circular Gaussian weighting
function to the largest galaxy, with a weight radius of $r_{w}$. The variance
of shape parameters is then plotted in Fig. 7 against the weight radius
$r_{w}$. By analyzing the visual illustration in Fig. 6 and the quantitative
results in Fig. 7, we can conclude that the Dirles and UPDC methods outperform
other methods in reconstructing shape parameters for extended sources. It is
worth noting that the Drizzle method produced the largest deviation from the
ground truth.

In the field of photometry, the reconstructed source profile plays a pivotal
role in estimating the performance of recovery algorithms. To analyze the
profiles for both the central star and the largest nearby galaxy, we have
presented graphical representations in Fig. 8 and Fig. 9, respectively. Our
analysis based on Fig. 6, 8, and Fig. 9 demonstrates that Dirles outperforms
other methods by effectively recovering high peak intensities in point
sources, providing superior profiles for extended sources, and minimizing
undesirable ringing artefacts. It is noteworthy that, for a given recovery
method, the optimal number of iterations exhibits minimal variation across
different realizations.

Figure 7: Shape parameters comparison between the true and the recoveries: left for the $e_{+}$ and right for the $e_{\times}$. Figure 8: The left panel shows the normalized profiles of the central star drawn in different colours to differentiate between them: Drizzle (in Magenta), UPDC (in Red), Richardson-Lucy (in Green), Starck-Pantin (in Blue) and Dirles (in Black). The right panel is a zoomed-in version that focuses on the area around the zero flux. This view highlights the ringing effect and other details around the point source. Figure 9:  Profiles of the largest galaxy from ground truth (in Cyan) and five recoveries: Drizzle (in Magenta), UPDC (in Red), Richardson-Lucy (in Green), Starck-Pantin (in Blue) and Dirles (in Black) . Table 1: Important Variables Used in the Text Variable |  Description
---|---
$\mathcal{O}$ |  intrinsic, continuous surface brightness of the sky
$\mathcal{G}_{k}$ |  the $k$-th exposure
$\mathcal{T}$ |  fully-sampled high resolution (or definition) image
$\mathcal{P}_{k}$ |  PSF measured from the $k$-th frame
$\mathcal{W}$ |  weight for flux allocation
$\mathbb{G}$ |  sample for pixels on the same WCS position
$\beta$ |  down-sampling factor
$\theta$ |  noise at the pixel
$\mathbb{N}\\{\cdot\\}$ |  noise operator
$\mathbb{B}\\{\cdot\\}$ |  blurring operator
$\mathscr{P}_{\mathbb{R}^{+}}\\{\cdot\\}$ |  component-wise projection operator for constraints set
Table 2: Computational resources consumed by various coaddition approaches during each iteration Simu No. | Drizzle | Richardson-Lucy | Starck-Pantin | UPDC | Dirles
---|---|---|---|---|---
Mock-I ($512\times 512\times 50,\beta=2$) | 3.7s | 220.4ms | 7.9s | 7.8s | 272.2s
Mock-II ($128\times 128\times 100,\beta=4$) | 1.9s | 90.1ms | 3.8s | 3.9s | 144.3s
Mock-III ($101\times 101\times 100,\beta=2$) | 0.3s | 6.8ms | 0.6s | 0.6s | 5.4s

We have conducted a comprehensive study of the computation time required by
different coaddition methods during each iteration. The results are presented
in Table 2. All calculations were performed on a two-socket AMD server
equipped with 2 CPUs of EPYC<EMAIL_ADDRESS>and 1TB DDR4 memory.

The Richardson-Lucy method has the lowest computation cost as it mainly
executes the FFT in the calculation, which increases with an $O(N{\rm log}N)$
trend. On the other hand, Drizzle consumes the most time on the exposure up-
sampling process. Both Starck-Pantin and UPDC involve similar operations such
as PSF convolving (via FFT), down-sampling, comparing, and up-sampling. Hence,
they have almost the same computation cost, which is mainly contributed by the
sampling.

Dirles has the highest computation cost compared to the other methods. Instead
of the FFT-based convolution, its PSF convolution is performed in real space,
which significantly increases computational time. Apart from the size of input
exposures, the computation also depends on the dimension of the PSF provided
in advance. However, the Dirles method stands out from other approaches as it
does not assume local uniformity for the PSF. Consequently, it is a better fit
for scenarios where the PSF undergoes significant changes within an exposure
space. With Dirles, we can convolve a PSF map that is different everywhere in
the field of view. Liu Dezi et al. (2023, in preparation) are constructing
such a PSF map by using a set of basic functions for fitting.

##  5 Conclusions

In this article, we propose a novel AAPD method, named Dirles, which is based
on the least squares fitting technique. Following a series of rigorous
simulation tests conducted on these algorithms, our findings can be summarized
as follows.

  1. 1.

In detecting faint sources, Dirles recovered all 15 substructures while
Drizzle missed 8.

  2. 2.

The Dirles approach outperforms previous methods by achieving the highest
PSNR, recovering object shapes, minimizing the ringing effect and suppressing
noise. However, it may ignore low SNR regions. To address this, gradually
reduce iterations and balance noise reduction with signal enhancement.

  3. 3.

The Dirles method outperforms other methods in reconstructing shape parameters
for extended sources. In radiometry, the Dirles method performs the best in
terms of recovering the highest peak in point sources, providing the best
profile for extended sources, and reducing the ringing effect.

  4. 4.

Compared to the other methods, Dirles has the highest computation cost.
However, it stands out from other approaches as it does not assume local
uniformity for the PSF, making it a better fit for scenarios where the PSF
undergoes significant changes within an exposure space.

After extensive simulations, it can be concluded that Dirles outperforms
previous works in restoring point/extended sources, maintaining extended
source shapes, and reducing ringing despite high computation consumption.

In the future, numerous upcoming telescopes will commence astronomical
observations, such as NASA’s Wide Field Infrared Survey Telescope (WFIRST),
the European Space Agency’s Euclid mission, the National Science Foundation-
funded Large Synoptic Survey Telescope (LSST), and China’s Space Station
Optical Telescope (CSST). These advanced instruments will generate a vast
volume of imaging data. Effectively processing these images while maintaining
high fidelity is an imminent necessity. As the implementation of Dirles
algorithm uses a divide-and-conquer approach, it may be possible to accelerate
it to a considerable extent using GPU-based parallel computing. We are
confident that there is a significant potential for improving the speed and
effectiveness of the Dirles method even further.

## Acknowledgements

We are grateful to the anonymous reviewers for their helpful comments, which
greatly improved the presentation of this work. This work was supported by the
GHfund A(202302017475). This work is also supported by the Foundation for
Distinguished Young Scholars of Jiangsu Province (No. BK20140050), the
National Natural Science Foundation of China (Nos. 11973070, 11333008,
11273061, 11825303, and 11673065), the China Manned Space Project with No.
CMS-CSST-2021-A01, CMS-CSST-2021-A03, CMS-CSST-2021-B01 and the Joint Funds of
the National Natural Science Foundation of China (No. U1931210). HYS
acknowledges the support from Key Research Program of Frontier Sciences, CAS,
Grant No. ZDBS-LY-7013 and Program of Shanghai Academic/Technology Research
Leader. We acknowledge the support from the science research grants from the
China Manned Space Project with CMS-CSST-2021-A04, CMS-CSST-2021-A07.

In this study, a cluster is used with the SIMT accelerator made in China. The
cluster includes many nodes each containing 2 CPUs and 4 accelerators. The
accelerator adopts a GPU-like architecture consisting of a 16GB HBM2 device
memory and many compute units. Accelerators connected to CPUs with PCI-E, the
peak bandwidth of the data transcription between main memory and device memory
is 16GB/s.

## References

  * Babacan et al. (2011) Babacan, S. D., Molina, R., & Katsaggelos, A. K. 2011, IEEE Trans Image Process, 20, 984
  * Bates & Cady (1980) Bates, R., & Cady, F. 1980, Optics Communications, 32, 365
  * Capel & Zisserman (2003) Capel, D., & Zisserman, A. 2003, IEEE Signal Processing Magazine, 20, 75
  * Elad & Feuer (1999) Elad, M., & Feuer, A. 1999, IEEE Trans On Pattern Analysis and Machine Intelligence (PAMI), 21, 817
  * Farsiu et al. (2004a) Farsiu, S., Robinson, D., Elad, M., & Milanfar, P. 2004a, International Journal of Imaging Systems and Technology, 14, 47
  * Farsiu et al. (2004b) Farsiu, S., Robinson, M., Elad, M., & Milanfar, P. 2004b, IEEE Transactions on Image Processing, 13, 1327
  * Fruchter (2011) Fruchter, A. S. 2011, Publications of the Astronomical Society of the Pacific, 123, 497
  * Fruchter & Hook (2002) Fruchter, A. S., & Hook, R. N. 2002, Publications of the Astronomical Society of the Pacific, 114, 144
  * Gerchberg (1974) Gerchberg, R. 1974, Optica Acta: International Journal of Optics, 21, 709
  * Hirata & Seljak (2003) Hirata, C., & Seljak, U. 2003, Monthly Notices of the Royal Astronomical Society, 343, 459
  * Högbom & Cornwell (1974) Högbom, J. A., & Cornwell, T. J. 1974, Astronomy and Astrophysics, 500, 55
  * Hunt (1994) Hunt, B. 1994, International Journal of Modern Physics C, 05, 151
  * Irani & Peleg (1993) Irani, M., & Peleg, S. 1993, J Vis. Commun. Image Represent, 4, 324
  * Landweber (1951) Landweber, L. 1951, Am. J. Math., 73, 615
  * Lauer (1999) Lauer, T. R. 1999, Publications of the Astronomical Society of the Pacific, 111, 227
  * Liu & Sun (2014) Liu, C., & Sun, D. 2014, IEEE Trans. Pattern Anal. Mach. Intell., 36, 346
  * Lucy (1974) Lucy, L. B. 1974, AJ, 79, 745
  * Nasrollahi & Moeslund (2014) Nasrollahi, K., & Moeslund, T. 2014, Machine Vision and Applications, 25, 1423
  * Ng et al. (2007) Ng, M. K., Shen, H., Lam, E. Y., & Zhang, L. 2007, EURASIP Journal on Advances in Signal Processing, 2007, 074585
  * Ngolè Mboula, F. M. et al. (2015) Ngolè Mboula, F. M., Starck, J.-L., Ronayette, S., Okumura, K., & Amiaux, J. 2015, A&A, 575
  * Nie et al. (2021a) Nie, L., Li, G., Peterson, J. R., & Wei, C. 2021a, Monthly Notices of the Royal Astronomical Society, 503, 4436
  * Nie et al. (2021b) Nie, L., Li, G., Zhang, J., Fan, Z., & Peterson, J. R. 2021b, Monthly Notices of the Royal Astronomical Society, 508, 3785
  * Park et al. (2003) Park, S. C., Park, M. K., & Kang, M. G. 2003, IEEE Signal Processing Magazine, 20, 21
  * Richardson (1972) Richardson, W. H. 1972, J. Opt. Soc. Am. A, 62, 55
  * Rowe et al. (2011) Rowe, B., Hirata, C., & Rhodes, J. 2011, The Astrophysical Journal, 741, 46
  * Rowe et al. (2015) Rowe, B., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121
  * Sage et al. (2017) Sage, D., Donati, L., Soulez, F., et al. 2017, Methods, 115, 28
  * Shepp & Vardi (1982) Shepp, L. A., & Vardi, Y. 1982, IEEE Trans Med Imaging, 1, 113
  * Starck et al. (2002) Starck, J. L., Pantin, E., & Murtagh, F. 2002, Publications of the Astronomical Society of the Pacific, 114, 1051
  * Su et al. (2012) Su, H., Tang, L., & Wu, Y. 2012, IEEE Trans Image Process, 21, 1031
  * Symons et al. (2021) Symons, T., Zemcov, M., Bock, J., et al. 2021, The Astrophysical Journal Supplement Series, 252, 24
  * Takeda et al. (2006) Takeda, H., Farsiu, S., Christou, J., & Milanfar, P. 2006, in The Maui Economic Development Board, ed. S. Ryan, Vol. id.E27, Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference
  * Takeda et al. (2007) Takeda, H., Farsiu, S., & Milanfar, P. 2007, IEEE Trans Image Process, 16, 349
  * Takeda et al. (2009) Takeda, H., Milanfar, P., & Protter, M. 2009, IEEE Trans Image Process, 18, 1958
  * Tian & Ma (2011) Tian, J., & Ma, K.-K. 2011, Signal, Image and Video Processing, 5, 329
  * van Ouwerkerk (2006) van Ouwerkerk, J. 2006, Image and Vision Computing, 24, 1039
  * Wang et al. (2022) Wang, L., Li, G., & Kang, X. 2022, Monthly Notices of the Royal Astronomical Society, 517, 787
  * Wang & Li (2017) Wang, L., & Li, G.-L. 2017, Research in Astronomy and Astrophysics, 17, 100
  * Yuan et al. (2010) Yuan, Q., Zhang, L., & Shen, H. 2010, IEEE Trans Image Process, 19, 3157
  * Yue et al. (2016) Yue, L., Shen, H., Li, J., et al. 2016, Signal Process., 128, 389
  * Zhang et al. (2012) Zhang, H., Zhang, L., & Shen, H. 2012, Signal Process, 92, 2082


# SCATTERING OF NEUTRINOS BY A ROTATING BLACK HOLE ACCOUNTING FOR THE
ELECTROWEAK INTERACTION WITH AN ACCRETION DISK

MAXIM DVORNIKOV  Pushkov Institute of Terrestrial Magnetism, Ionosphere
and Radiowave Propagation (IZMIRAN),
108840 Moscow, Troitsk, Russia
<EMAIL_ADDRESS>
###### Abstract

We study spin effects in the neutrino gravitational scattering by a
supermassive black hole with a magnetized accretion disk having a finite
thickness. We exactly describe the propagation of ultrarelativistic neutrinos
on null geodesics and solve the spin precession equation along each neutrino
trajectory. The interaction of neutrinos with the magnetic field is owing to
the nonzero diagonal magnetic moment. Additionally, neutrinos interact with
plasma of the accretion disk electroweakly within the Fermi approximation.
These interactions are obtained to change the polarization of incoming
neutrinos, which are left particles. The fluxes of scattered neutrinos,
proportional to the survival probability of spin oscillations, are derived for
various parameters of the system. In particular, we are focused on the matter
influence on the outgoing neutrinos flux. The possibility to observe the
predicted effects for astrophysical neutrinos is briefly discussed.

##  1 Introduction

Neutrino interactions with external fields can significantly modify the
dynamics of neutrino oscillations, which were recently confirmed
experimentally (see, e.g., Ref. Ace22). We deal here mainly with neutrino spin
oscillations, in which active left handed neutrinos become sterile right
particles under the influence of some external fields. Gravitational
interaction, in spite of its weakness, can also contribute to neutrino
oscillations (see, e.g., Ref. PirRoyWud96).

The propagation of a spinning particle in curved spacetime was first
considered in Ref. Pap51. Then, the quasiclassical approach for the
description of the spin evolution of a point-like elementary particle in a
gravitational field was formulated in Ref. PomKhr98. This method was applied
for the studies of neutrino spin oscillations in a curved spacetime both in
frames of the General Relativity (GR) in Refs. AlaNod13,BayPen21 and in
various alternative gravity theories in Refs. AlaNod15,Cha15,MasLam21. The
perturbative quantum treatment of neutrino spin oscillations in the
gravitational field of a black hole (BH) was carried out in Refs.
SorZil07,Sor12. Quasiclassical and quantum approaches for the studies of
spinning particles in a gravitational field were reviewed in Ref. Ver22.

We considered neutrino spin oscillations in curved spacetimes in Refs.
Dvo06,Dvo13,Dvo20a,Dvo20,Dvo21,Dvo23,Dvo19. Both stationary [13, 14, 15, 16,
17, 18] and time dependent metrics, like a gravitational wave [19], were
analyzed. We studied spin effects in the neutrino gravitational scattering by
BH in Refs. Dvo20a,Dvo20,Dvo21,Dvo23. In a scattering problem, both in- and
out-states are in the asymptotically flat spacetime. Thus, we can attribute
certain polarization to these states. We studied the cases of both nonrotating
and rotating BHs, as well as analyzed the contribution of various external
fields to spin oscillations of scattered neutrinos. In the present work, we
continue to examine this problem.

The recent observations of the shadows of supermassive BHs (SMBHs) in the
centers of M87 and our Galaxy by the Event Horizon Telescope (EHT) [20, 21]
were the motivation for the present research. Those results are the unique
tests of GR in the strong field limit. There are searches of the high energy
neutrinos emission from the centers of some active galaxies [22]. If such
neutrinos are produced in the vicinity of the SMBH surface, they are strongly
lensed analogously to photons observed by the EHT collaboration. External
fields, which can exist near SMBHs, will change the neutrino polarization
leading to the distortion of the outgoing neutrinos flux. Alternatively,
neutrinos emitted in a core-collapsing supernova (SN) explosion can be
gravitationally lensed by the SMBH in the center of our Galaxy. Such
possibility was discussed in Refs. MenMocQui06,Los21. In this situation, we
can also expect the influence of external fields in the curved spacetime on
the dynamics of neutrino spin oscillations.

This work is organized in the following way. We start in Sec. 2 with a brief
description of the ultrarelativistic neutrinos motion in the Kerr metric. We
also describe the neutrino spin evolution in background matter under the
influence of electromagnetic and gravitational fields. Then, in Sec. 3, we set
up the neutrino interaction with plasma of an accretion disk. The rest of the
parameters of the system is chosen in Sec. 4. We present the results for the
calculation of the outgoing neutrino fluxes, which account for neutrino spin
oscillations, in Sec. 5. The conclusion is given in Sec. 6 We find the
component of the four velocity of plasma in the accretion disk in A. Some
computational details are provided in B.

##  2 Motion of a neutrino in the Kerr metric and its spin evolution

We study the propagation of a test neutrino in the gravitational field of a
rotating BH which is described by the Kerr metric,

$\mathrm{d}s^{2}=g_{\mu\nu}\mathrm{d}x^{\mu}\mathrm{d}x^{\nu}=\left(1-\frac{rr_{g}}{\Sigma}\right)\mathrm{d}t^{2}+2\frac{rr_{g}a\sin^{2}\theta}{\Sigma}\mathrm{d}t\mathrm{d}\phi-\frac{\Sigma}{\Delta}\mathrm{d}r^{2}-\Sigma\mathrm{d}\theta^{2}-\frac{\Xi}{\Sigma}\sin^{2}\theta\mathrm{d}\phi^{2},$
(2.1)

where the Boyer-Lindquist coordinates $x^{\mu}=(t,r,\theta,\phi)$ are
utilized. In Eq. (2.1), we use the following notations:

$\Delta=r^{2}-rr_{g}+a^{2},\quad\Sigma=r^{2}+a^{2}\cos^{2}\theta,\quad\Xi=\left(r^{2}+a^{2}\right)\Sigma+rr_{g}a^{2}\sin^{2}\theta.$
(2.2)

The mass of BH in Eqs.(2.1) and (2.2) is $M=r_{g}/2$ and its angular momentum
is $J=Ma$, where $0<a<M$. The BH spin is directed upward from the equatorial
plane $\theta=\pi/2$.

The arbitrary motion of a test ultrarelativistic particle in the Kerr metric
is characterized by two conserved quantities: the angular momentum $L$ and the
Carter constant $Q$. If we study the scattering problem, $Q>0$. The law of
motion and the form of a trajectory for such a particle can be found in
quadratures. The corresponding quantities are expressed in terms of elliptic
integrals. The detailed description of this problem can be found, e.g., in
Ref. Cha83.

If a test particle is spinning, like a neutrino, the invariant equation for
the evolution of the particle spin $S^{\mu}$ in a curved spacetime under the
influence of the electromagnetic field $F_{\mu\nu}$ and the background matter
has the form [14],

$\displaystyle\frac{\mathrm{D}S^{\mu}}{\mathrm{d}\tau}=$ $\displaystyle
2\mu\left(F^{\mu\nu}S_{\nu}-U^{\mu}U_{\nu}F^{\nu\lambda}S_{\lambda}\right)+\sqrt{2}G_{\mathrm{F}}E^{\mu\nu\lambda\rho}G_{\nu}U_{\lambda}S_{\rho},$
(2.3)

where $U^{\mu}=\tfrac{\mathrm{d}x^{\mu}}{\mathrm{d}\tau}$ is the four velocity
in the world coordinates, $\tau$ is the proper time, $\mathrm{D}S^{\mu}$ is
the covariant differential,
$E^{\mu\nu\lambda\rho}=\tfrac{1}{\sqrt{-g}}\varepsilon^{\mu\nu\lambda\rho}$ is
the covariant antisymmetric tensor in a curved spacetime, $g=\det(g_{\mu\nu})$
is the determinant of the metric tensor, $\mu$ is the neutrino magnetic
moment, $G_{\mathrm{F}}=1.17\times 10^{-5}\,\text{GeV}^{-2}$ is the Fermi
constant, and $G_{\mu}$ is the vector which incorporates the characteristics
of the background matter. We discuss its form in details in Sec. 3 below.
Although Eq. (2.3) is written down for massive particles, it has the
appropriate limit for an ultrarelativistic neutrino. In Eq. (2.3), we suppose
that a neutrino is a Dirac particle. Despite some claims about the Majorana
nature of neutrinos, this issue is still unresolved [26].

It is more convenient to describe the evolution of the neutrino polarization
in a locally Minkowskian frame $x_{a}=e_{a}^{\ \mu}x_{\mu}$, where $e_{a}^{\
\mu}$ are the vierbein vectors which have the form,

$\displaystyle e_{0}^{\ \mu}=$
$\displaystyle\left(\sqrt{\frac{\Xi}{\Sigma\Delta}},0,0,\frac{arr_{g}}{\sqrt{\Delta\Sigma\Xi}}\right),\quad
e_{1}^{\ \mu}=\left(0,\sqrt{\frac{\Delta}{\Sigma}},0,0\right),$ $\displaystyle
e_{2}^{\ \mu}=$ $\displaystyle\left(0,0,\frac{1}{\sqrt{\Sigma}},0\right),\quad
e_{3}^{\
\mu}=\left(0,0,0,\frac{1}{\sin\theta}\sqrt{\frac{\Sigma}{\Xi}}\right).$ (2.4)

These vectors diagonalize the metric in Eq. (2.1), $g_{\mu\nu}e_{a}^{\
\mu}e_{b}^{\ \nu}=\eta_{ab}$, where $\eta_{ab}=(1,-1,-1,-1)$ is the Minkowski
metric tensor. If we define the three vector of the polarization $\bm{\zeta}$
in the particle rest frame in the locally Minkowskian frame, it obeys the
equation,

$\frac{\mathrm{d}\bm{\bm{\zeta}}}{\mathrm{d}t}=2(\bm{\bm{\zeta}}\times\bm{\bm{\Omega}}),\quad\bm{\bm{\Omega}}=\bm{\bm{\Omega}}_{g}+\bm{\bm{\Omega}}_{\mathrm{em}}+\bm{\bm{\Omega}}_{\mathrm{matt}}$
(2.5)

where

$\displaystyle\bm{\bm{\Omega}}_{g}$
$\displaystyle=\frac{1}{2U^{t}}\left[\mathbf{b}_{g}+\frac{1}{1+u^{0}}\left(\mathbf{e}_{g}\times\mathbf{u}\right)\right],$
$\displaystyle\bm{\bm{\Omega}}_{\mathrm{em}}$
$\displaystyle=\frac{\mu}{U^{t}}\left[u^{0}\mathbf{b}-\frac{\mathbf{u}(\mathbf{u}\mathbf{b})}{1+u^{0}}+(\mathbf{e}\times\mathbf{u})\right],$
$\displaystyle\bm{\bm{\Omega}}_{\mathrm{matt}}$
$\displaystyle=\frac{G_{\mathrm{F}}}{\sqrt{2}U^{t}}\left[\mathbf{u}\left(g^{0}-\frac{(\mathbf{gu})}{1+u^{0}}\right)-\mathbf{g}\right].$
(2.6)

Here $u^{a}=(u^{0},\mathbf{u})=e_{\ \mu}^{a}U^{\mu}$,
$G_{ab}=(\mathbf{e}_{g},\mathbf{b}_{g})=\gamma_{abc}u^{c}$,
$\gamma_{abc}=\eta_{ad}e_{\ \mu;\nu}^{d}e_{b}^{\ \mu}e_{c}^{\ \nu}$ are the
Ricci rotation coefficients, the semicolon stays for the covariant derivative,
$f_{ab}=e_{a}^{\ \mu}e_{b}^{\ \nu}F_{\mu\nu}=(\mathbf{e},\mathbf{b})$ is the
electromagnetic field tensor in the locally Minkowskian frame, and
$g^{a}=(g^{0},\mathbf{g})=e_{\ \mu}^{a}G^{\mu}$. The details of the derivation
of Eqs. (2.5) and (2) can be found in Ref. Dvo13.

Instead of dealing with Eq. (2.5), we can rewrite it in the form of the
effective Schrödinger equation $\mathrm{i}\dot{\psi}=H\psi$, where
$H=-(\bm{\sigma}\bm{\Omega})$ and $\bm{\sigma}$ are the Pauli matrices. It is
more convenient to label a point on a neutrino trajectory with $r$ rather than
with $t$. Thus, we transform the Schrödinger equation to the form,

$\mathrm{i}\frac{\mathrm{d}\psi}{\mathrm{d}r}=H_{r}\psi,\quad
H_{r}=-\frac{\mathrm{d}t}{\mathrm{d}r}\mathcal{U}_{2}(\bm{\sigma\Omega})\mathcal{U}_{2}^{\dagger},$
(2.7)

where the matrix $\mathcal{U}_{2}=\exp(\mathrm{i}\sigma_{2}\pi/4)$ accounts
for the fact that incoming and outgoing neutrinos move oppositely and along
the first axis in the locally Minkowskian frame. The derivative
$\mathrm{d}t/\mathrm{d}r$ is taken on the basis of the law of motion of
neutrinos; cf., e.g., Refs. Dvo23,Cha83.

The initial spin wavefunction is $\psi_{-\infty}^{\mathrm{T}}=(1,0)$. It
corresponds to a left polarized neutrino. After reconstructing the neutrino
trajectory and solving Eq. (2.7), we get the wave function
$\psi_{+\infty}^{\mathrm{T}}=(\psi_{+\infty}^{(\mathrm{R})},\psi_{+\infty}^{(\mathrm{L})})$
of an outgoing neutrino. The survival probability for a neutrino to remain
left polarized after the scattering is
$P_{\mathrm{LL}}=\left|\psi_{+\infty}^{(\mathrm{L})}\right|^{2}$. Here we
account for the fact that the neutrino velocity changes its direction in the
locally Minkowskian frame.

##  3 Neutrino interaction with background matter

In this section, we discuss in details the electroweak interaction of
scattering neutrinos with plasma of an accretion disk around BH. The
properties of the accretion disk are also considered.

We suppose that the accretion disk consists of the hydrogen plasma. The
neutrino interaction with electrons and protons is treated within the Fermi
model in the forward scattering approximation. The effective Lagrangian for
the interaction of the neutrino bispinor $\nu$ with matter has the form,

$\mathcal{L}_{\mathrm{matt}}=-\frac{G_{\mathrm{F}}}{\sqrt{2}}\bar{\nu}\gamma^{\mu}(1-\gamma^{5})\nu\cdot
G_{\mu},$ (3.1)

where $\gamma^{\mu}$ and $\gamma^{5}$ are the Dirac matrices. The matter
characterictics are in the four vector potential $G^{\mu}$ which has the form
[14],

$G^{\mu}=\sum_{f=e,p}\left(q_{f}^{(1)}J_{f}^{\mu}+q_{f}^{(2)}\Lambda_{f}^{\mu}\right),$
(3.2)

where $J_{f}^{\mu}$ is the invariant hydrodynamics current of plasma fermions,
$\Lambda_{f}^{\mu}$ is the plasma invariant polarization, and [27]

$q_{f}^{(1)}=I_{\mathrm{L}3}^{(f)}-2Q_{f}\sin^{2}\theta_{\mathrm{W}}+\delta_{ef},\quad
q_{f}^{(2)}=-I_{\mathrm{L}3}^{(f)}-\delta_{ef}.$ (3.3)

Here $I_{\mathrm{L}3}^{(f)}$ is the third component of the weak isospin of
background fermions, $Q_{f}$ is the value of their electric charge,
$\theta_{\mathrm{W}}$ is the Weinberg angle, $\delta_{ef}=1$ for electrons and
vanishes for protons. The coefficients $q_{f}^{(1,2)}$ in Eq. (3.3) correspond
to the scattering of electron neutrinos.

The plasma motion in an accretion disk around a rotating BH is quite complex
[28]. A circular orbit is possible only when an accretion disk is thin and is
situated in the BH equatorial plane. Since we study the general gravitational
scattering of neutrinos, which involves the neutrino motion both above and
below the equatorial plane (see Sec. 4 below), and consider the accretion disk
with a nonzero thickness, we assume that plasma in the accretion disk is
nonrelativistic and unpolarized. In this case, only
$J_{f}^{0}=n_{f}U_{f}^{t}\neq 0$, whereas $J_{f}^{i}=0$ and
$\Lambda_{f}^{\mu}=0$. Here $n_{f}$ is the invariant number density of
background fermions, measured by a comoving observer, and $U_{f}^{t}$ is the
time component of their four velocity. We compute $U_{f}^{t}$ in A; see Eq.
(A.4). Eventually, we get that the only nonzero component of $G^{\mu}$ in Eq.
(3.2) is $G^{t}=n_{e}U_{f}^{t}$, where $n_{e}$ is the electron number density
and we suppose that plasma is electroneutral.

Using Eqs. (2.4) and (A.4), we obtain that the vector $g^{a}=e_{\
\mu}^{a}G^{\mu}$ has the form,

$g^{a}=\frac{n_{e}U_{f}^{t}}{\sqrt{\Xi}}\left(\sqrt{\Sigma\Delta},0,0,-\frac{arr_{g}\sin\theta}{\sqrt{\Sigma}}\right).$

The distribution of matter in an accretion disk is quite model dependent. We
suppose that [29]

$n_{e}(r,\theta)=n_{e}^{(c)}\exp\left(-\frac{r^{2}\cos^{2}\theta}{2H^{2}}\right),$
(3.4)

where $H$ is the disk thickness and $n_{e}^{(c)}$ is the central density, i.e.
the density at the equatorial plane. We take that [30] $n_{e}^{(c)}(r)\propto
n_{0}r^{-3/2}$, where $n_{0}=10^{18}\,\text{cm}^{-3}$ for SMBH with
$M=10^{8}M_{\odot}$ [31].

The matter contribution to the vector $\bm{\Omega}$, which determines the spin
precession of ultrarelativistic neutrinos in Eq. (2), takes the form,

$\bm{\Omega}_{\mathrm{matt}}=\frac{G_{\mathrm{F}}}{\sqrt{2}U^{t}}\mathbf{u}\left(g^{0}-\frac{(\mathbf{gu})}{1+u^{0}}\right),$
(3.5)

where we neglect the term $\propto\mathbf{g}/U^{t}$ for ultrarelativistic
neutrinos. Nevertheless both $\mathbf{u}/U^{t}$ and $\mathbf{u}/(1+u^{0})$ are
finite in Eq. (3.5).

##  4 Parameters of the system

We study the gravitational scattering of ultrarelativistic neutrinos by a
rotating SMBH surrounded with a magnetized accretion disk with a finite
thickness. The incoming flux of neutrinos is parallel to the equatorial plane
of BH. Nevertheless, unlike Refs. Dvo20a,Dvo20,Dvo21, we do not restrict
ourselves by the equatorial neutrino motion. Incoming neutrinos are emitted
from the direction $(\theta,\phi)_{s}=(\pi/2,0)$. A remote neutrino detector
is in the arbitrary position with the angular coordinates
$(\theta_{\mathrm{obs}},\phi_{\mathrm{obs}})$.

The mass of SMBH is taken to be $10^{8}M_{\odot}$. The accretion disk around
this SMBH consists of hydrogen plasma with the electron density distribution
in Eq. (3.4) (see also Refs. SteSpr99,Jia19). The main uncertainty among the
disk parameters is its thickness $H$. We vary it in the range $(1-10)r_{g}$ in
our simulations (see, e.g., Ref. Sad09), considering disks with a nonzero
thickness. We study the simplified model of the accretion disk which does not
account for the disk rotation. The neutrino interaction with plasma is within
the Fermi approximation of the standard model. We consider only the forward
scattering of electron neutrinos on plasma fermions.

We take into account the poloidal component of the magnetic field which is
based on the following vector potential in the world coordinates [33]:

$A_{t}=Ba\left[1-\frac{rr_{g}}{2\Sigma}(1+\cos^{2}\theta)\right],\quad
A_{\phi}=-\frac{B}{2}\left[r^{2}+a^{2}-\frac{a^{2}rr_{g}}{\Sigma}(1+\cos^{2}\theta)\right]\sin^{2}\theta,$
(4.1)

where the amplitude $B$ is supposed to scale with radius as $B\propto
B_{0}r^{-5/4}$ [34]. The strength $B_{0}$, which is the magnetic field at the
inner radius of the disk, is taken to be $B_{0}=3.2\times 10^{2}\,\text{G}$.
This value of $B_{0}$ is below the Eddington limit for $M=10^{8}M_{\odot}$
[35]. The toroidal magnetic field, which is inevitably generated in a thick
disk and can be rather strong, is not considered in our model. We plan to
account for both the realistic plasma motion and the toroidal magnetic field
in one of the forthcoming works.

A neutrino is taken to be a Dirac particle with a nonzero magnetic moment
$\mu$. We consider the values of $\mu$ in the range
$\mu=(10^{-14}-10^{-13})\mu_{\mathrm{B}}$. Such magnetic moments are within
the theoretical and astrophysical constraints on neutrino magnetic moments
established in Refs. Bel05,Via13.

To reconstruct the neutrino motion, first, we find the $\theta(r)$ dependence.
One can find the details of this procedure in Refs. Dvo23,Cha83. The knowledge
of $\phi(r)$ is not necessary since the spin evolution does not depend on
$\phi$. We just need to obtain the final scattering angle
$\phi_{\mathrm{obs}}$. This fact significantly accelerates numerical
simulations allowing us to involve more test particles and perform more
frequent meshing in the radial direction. Now we use $2500$ test neutrinos in
the incoming flux. However, some of these particles are not involved in the
scattering since they fall to the BH shadow region.

The right hand side of Eq. (2.7) is given only in the discrete nodes $r_{i}$.
Moreover, the grid $r_{i}$ is irregular. Thus, we apply the two-step
Adams–Bashforth method for the numerical integration of Eq. (2.7). The
corresponding iterative procedure is provided in B; see Eqs. (B.4) and (B.5).
It allows us to increase the accuracy of numerical simulations compared to
Ref. Dvo23 where the Euler method was used.

##  5 Results

In this section we present the results of the numerical solution of Eq. (2.7)
accounting for the values of the parameters in Sec. 4.

Ultrarelativistic neutrinos are emitted as left polarized particles in frames
of the standard model. A terrestrial neutrino telescope can detect only left
neutrinos as well. Thus, if $F_{0}\propto\mathrm{d}\sigma/\mathrm{d}\varOmega$
is the flux of scalar test particles, the observed flux of neutrinos is
$F_{\nu}=P_{\mathrm{LL}}F_{0}$. Here $\mathrm{d}\sigma/\mathrm{d}\varOmega$ is
the differential cross section of the gravitational scattering. Our main goal
is to find the ratio $F_{\nu}/F_{0}$ for neutrinos gravitationally scattered
by BH and to account for the neutrino spin precession in all external fields.
We normalize the neutrino flux by $F_{0}$ in order to avoid the singularities
in the cross section which are inherent in the gravitational scattering of
both spinning and spinless particles [38]. The detailed study of the
gravitational scattering of scalar test particles off a rotating BH can be
found in Ref. Boz08. We utilize the term ‘scalar particles’ meaning their
motion along null geodesics.

In the present work, we use the modified method for the integration of Eq.
(2.7); cf. B. We should check the behavior of the solution when only the
gravitational interaction is taken into account. We expect that
$F_{\nu}/F_{0}\approx 1$, i.e. no spin oscillations occur in a purely
gravitational scattering. The corresponding general theorem has been proven in
Ref. Dvo23. We show $F_{\nu}/F_{0}=P_{\mathrm{LL}}$ in Fig. 1 for BHs with
different angular momenta.

Figure 1: The ratio of the fluxes of spinning ultrarelativistic neutrinos
$F_{\nu}$ and the scalar particles $F_{0}$ scattered off BH with different
spins when only the gravitational interaction is accounted for. Panel (a):
$a=0.02M$; panel (b): $a=0.98M$.

One can see in Fig. 1 that the deviation of $F_{\nu}/F_{0}$ from one is less
than $1\%$ for the majority of $(\theta_{\mathrm{obs}},\phi_{\mathrm{obs}})$.
The higher accuracy of simulations compared to that in Ref. Dvo23 is owing to
the two-step Adams–Bashforth method represented in B. The result that
$F_{\nu}/F_{0}\approx 1$ means that our simulations are reliable.

One has white areas Fig. 1 (see also Figs. 2 and 3 below), which are present
because of the insufficient number of neutrinos scatter to these regions.
Thus, 2D cubic interpolation, used in our simulations, is unable to create
contour plots there. This shortcoming may be eliminated by a significant
enhancement of the test particles number. We plan to implement this task in a
future work.

Now, we add the neutrino interaction with the magnetic field and the accretion
disk to Eq. (2.7). First, we consider an almost nonrotating SMBH with
$a=0.02M$. We show the ratio $F_{\nu}/F_{0}$ in Fig. 2. We plot
$F_{\nu}/F_{0}$ in case when the magnetic interaction is dominant, i.e.
$V_{m}=0$, in Fig. 2 for $\mu=10^{-14}\mu_{\mathrm{B}}$ and in Fig. 2 for
$\mu=10^{-13}\mu_{\mathrm{B}}$. In both cases, $B_{0}=3.2\times
10^{2}\,\text{G}$. Then, in Figs. 2-2, we account for a nonzero interaction
with background matter by setting the electron number density near the SMBH
surface to $n_{0}=10^{18}\,\text{cm}^{-3}$. Figures 2 and 2 correspond to the
disk thickness $H=r_{g}$, whereas Figs. 2 and 2 to $H=10r_{g}$.

Figure 2: The ratio of the fluxes $F_{\nu}/F_{0}$ for ultrarelativistic
$\nu_{e}$ scattered off the almost nonrotating SMBH with $a=0.02M$ and
$M=10^{8}M_{\odot}$. Panel (a): $V_{\mathrm{B}}=\mu B_{0}r_{g}=2.7\times
10^{-2}$, $V_{m}=G_{\mathrm{F}}n_{0}r_{g}/\sqrt{2}=0$. Panel (b):
$V_{\mathrm{B}}=2.7\times 10^{-1}$, $V_{m}=0$. Panel (c):
$V_{\mathrm{B}}=2.7\times 10^{-2}$, $V_{m}=10^{-1}$, $H=r_{g}$. Panel (d):
$V_{\mathrm{B}}=2.7\times 10^{-1}$, $V_{m}=10^{-1}$, $H=r_{g}$. Panel (e):
$V_{\mathrm{B}}=2.7\times 10^{-2}$, $V_{m}=10^{-1}$, $H=10r_{g}$. Panel (f):
$V_{\mathrm{B}}=2.7\times 10^{-1}$, $V_{m}=10^{-1}$, $H=10r_{g}$.

We can see in Figs. 2 and 2 that the neutrino interaction with plasma of the
accretion disk suppresses neutrino spin oscillations. We remind that spin
oscillations of neutrinos are in the resonance when only the magnetic
interaction is taken into account (see, e.g., Ref. VolVysOku86, where this
problem is discussed in the flat spacetime). The electroweak interaction with
matter shifts spin oscillations out of the resonance. It explains the behavior
of $P_{\mathrm{LL}}$ in Figs. 2 and 2. Moreover, the consideration of a
thicker disk in Figs. 2 and 2 results in the further enhancement of the
survival probability since the regions with the great $P_{\mathrm{LL}}$ become
wider. This result is explained by the fact that more scattered neutrinos
interact with matter if the disk is thicker.

Finally, we plot $F_{\nu}/F_{0}$ in Fig. 3 for an almost maximally rotating
SMBH with $a=0.98M$. Qualitatively, the behavior of the survival probability
resembles that in Fig. 2, which is discussed above. We just mention that the
interaction with matter is more pronounced for a small neutrino magnetic
moment. Indeed, one can see in Figs. 3 and 3 that neutrino spin oscillations
are almost suppressed.

Figure 3: The same as in Fig. 2 for $a=0.98M$, i.e. SMBH is almost maximally
rotating.

##  6 Conclusion

In the present work, we have studied the scattering of ultrarelativistic
neutrinos off a rotating SMBH surrounded by a magnetized accretion disk with a
nonzero thickness. The flux of left polarized incoming neutrinos was taken to
propagate parallel to the equatorial plane of BH. However, we did not restrict
ourselves to a purely equatorial neutrino motion. These neutrinos have been
strongly gravitationally lensed towards a distant observer with the arbitrary
angular coordinates $0<\theta_{\mathrm{obs}}<\pi$ and
$0<\phi_{\mathrm{obs}}<2\pi$.

Unlike scalar particles or photons, the polarization of neutrinos is important
for their detection since we cannot observe sterile right neutrinos. That is
why, besides the reconstruction of the particle trajectories in the curved
spacetime, we had to study neutrino spin precession in the external fields. We
have supposed that a neutrino is a Dirac particle having a nonzero magnetic
moment, which provides the interaction with a poloidal magnetic field in the
accretion disk. The neutrino electroweak interaction with plasma in the disk
is considered in frames of the Fermi approximation.

In Sec. 2, we have outlined the neutrino motion in the Kerr metric and the
description of the spin evolution in external fields in curved spacetime. The
neutrino interaction with background matter and the model of the accretion
disk have been described in Sec. 3. The characteristics of the magnetic field
are the same as in Ref. Dvo23. Therefore we did not discuss them in details
here. Then, in Sec. 4, we have specified the rest of the parameters in the
system and considered some details of the numerical simulations.

We have presented our results in Sec. 5. The fluxes of spinning neutrinos have
been found for various parameters of the system. We have confirmed the general
theorem, proven in Ref. Dvo23, that the polarization of ultrarelativistic
neutrinos is unchanged in a purely gravitational scattering. Then, we have
obtained that the neutrino interaction with plasma of an accretion disk
suppresses neutrino spin oscillations. This feature, known in flat spacetime,
has been demonstrated for the neutrino gravitational scattering. Besides the
consideration of the accretion disk with a finite thickness, one of the
advantages of the present work, compared to Ref. Dvo23, was the more precise
numerical integration of Eq. (2.7). It allowed us to significantly increase
the accuracy of the simulations.

The derivation of the time component of the four velocity of matter in the
accretion disk has been provided in A. In B, we have described the more
precise method of the integration of Eq. (2.7) in a nonuniform grid.

The results obtained allow one to explore the characteristics of inner parts
of the accretion disk around SMBH in our Galaxy, e.g., with neutrinos emitted
in a SN explosion. The huge number of such neutrinos are expected to be
observed with the existing or future neutrino telescopes. If such neutrinos
are gravitationally lensed by the central SMBH, we have to account for the
interaction with strong external fields in curved spacetime near this SMBH.

##  Appendix A Derivation of $U^{t}$ for an arbitrary position of a plasma
particle

The action for a test particle moving in the Kerr metric has the form,
$S=-Et+L\phi+\dotsc,$ where we show only the terms responsible for the
conserved energy, $E$, and the projection of the particle angular momentum on
the BH spin, $L$. The canonical four momentum of a massive particle reads
$p_{\mu}=mU_{\mu}=-\partial_{\mu}S$. Thus, we have that

$m\left(g_{tt}U^{t}+g_{t\phi}U^{\phi}\right)=E,\quad m\left(g_{\phi
t}U^{t}+g_{\phi\phi}U^{\phi}\right)=-L.$ (A.1)

Solving Eq. (A.1), we get that [41]

$U^{t}=-\frac{(Eg_{\phi\phi}+Lg_{\phi t})}{m(g_{\phi
t}^{2}-g_{\phi\phi}g_{tt})}.$ (A.2)

Using the Kerr metric components in Eq. (2.1), one rewrites Eq. (A.2) in the
form,

$U^{t}=\frac{E\left[\left(r^{2}+a^{2}\right)\Sigma+rr_{g}a^{2}\sin^{2}\theta\right]-Lrr_{g}a\sin^{2}\theta}{m\Sigma\Delta\sin^{2}\theta}.$
(A.3)

If a particle moves in the equatorial plane with $\theta=\pi/2$, we reproduce
the expression for $U^{t}$ obtained in Ref. Dvo20.

Finally, if a particle is nonrelativistic, it has $E=m$ and $L=0$. We obtain,
basing on Eq. (A.3), that

$U^{t}=\frac{\Xi}{\Sigma\Delta\sin^{2}\theta},$ (A.4)

which is used in Sec. 3 as $U_{f}^{t}$.

##  Appendix B Two-step Adams–Bashforth method for an irregular grid

We discuss the solution of the single first order differential equation

$y^{\prime}=f(x,y).$ (B.1)

The generalization of the results to a system of equations is straightforward.
If the function $f(x,y)$ in Eq. (B.1) is given only in discrete points
$f_{i}$, we cannot use the precise Runge-Kutta method.

The value of the function at the $(i+2)$th node is

$y_{i+2}=y_{i+1}+\int_{x_{i+1}}^{x_{i+2}}f(x,y(x))\mathrm{d}x.$ (B.2)

We replace the integrand in Eq. (B.2) with the linear polynomial

$f(x,y(x))\to
p(x)=\frac{x-x_{i}}{x_{i+1}-x_{i}}f_{i+1}+\frac{x-x_{i+1}}{x_{i}-x_{i+1}}f_{i}.$
(B.3)

One can see that $p(x_{i})=f_{i}$ and $p(x_{i+1})=f_{i+1}$.

Integrating in Eq. (B.2) and accounting for Eq. (B.3), we get that

$y_{i+2}=y_{i+1}+C_{1}f_{i+1}+C_{0}f_{i},$ (B.4)

where

$C_{1}=\frac{(x_{i+2}-x_{i+1})}{2(x_{i+1}-x_{i})}(x_{i+2}+x_{i+1}-2x_{i}),\quad
C_{0}=\frac{(x_{i+2}-x_{i+1})}{2(x_{i+1}-x_{i})}(x_{i+1}-x_{i+2}).$ (B.5)

If the grid is regular with the step $h$, we obtain in Eq. (B.5) that
$C_{1}=\tfrac{3}{2}h$ and $C_{0}=-\tfrac{1}{2}h$. It reproduces the result
known previously (see, e.g., Ref. KorKor68).

Equations (B.4) and (B.5) are more precise than the Euler method. We use them
in Sec. 5 to solve Eq. (2.7), where the grid is irregular.

## References

  * [1] F. P. An, et al. (Daya Bay collaboration), arXiv:2211.14988.
  * [2] D. Píriz, M. Roy, and J. Wudka, Phys. Rev. D 54 (1996) 1587, hep-ph/9604403.
  * [3] A. Papapetrou, Proc. Roy. Soc. Lond. A 209 (1951) 248.
  * [4] A. A. Pomeranskiĭ and I. B. Khriplovich, J. Exp. Theor. Phys. 86 (1998) 839, gr-qc/9710098.
  * [5] S. A. Alavi and S. Nodeh, Grav. Cosmol. 19 (2013) 129, arXiv:1108.3593.
  * [6] G. Baym and J.-C. Peng, Phys. Rev. D 103 (2021) 123019, arXiv:2103.11209.
  * [7] S. A. Alavi and S. Nodeh, Phys. Scr. 90 (2015) 035301, arXiv:1301.5977.
  * [8] S. Chakraborty, J. Cosmol. Astropart. Phys. 10 (2015) 019, arXiv:1506.02647.
  * [9] L. Mastrototaro and G. Lambiase, Phys. Rev. D 104 (2021) 024021, arXiv:2106.07665.
  * [10] F. Sorge and S. Zilio, Class. Quantum Grav. 24 (2007) 2653.
  * [11] F. Sorge, Class. Quantum Grav. 29 (2012) 045002.
  * [12] S. N. Vergeles, N. N. Nikolaev, Yu. N. Obukhov, A. Ya. Silenko, and O. V. Teryaev, Phys.—Usp. 66 (2) (2023), arXiv:2204.00427.
  * [13] M. Dvornikov, Int. J. Mod. Phys. D 15 (2006) 1017, hep-ph/0601095.
  * [14] M. Dvornikov, J. Cosmol. Astropart. Phys. 06 (2013) 015 arXiv:1306.2659.
  * [15] M. Dvornikov, Phys. Rev. D 101 (2020) 056018, arXiv:1911.08317.
  * [16] M. Dvornikov, Eur. Phys. J. C 80 (2020) 474, arXiv:2006.01636.
  * [17] M. Dvornikov, J. Cosmol. Astropart. Phys. 04 (2021) 005, arXiv:2102.00806.
  * [18] M. Dvornikov, Class. Quantum Grav. 40 (2023) 015002, arXiv:2206.00042.
  * [19] M. Dvornikov, Phys. Rev. D 99 (2019) 116021, arXiv:1902.11285.
  * [20] K. Akiyama et al. (Event Horizon Telescope Collaboration), Astrophys. J. Lett. 875 (2019) L1, arXiv:1906.11238.
  * [21] K. Akiyama et al. (Event Horizon Telescope Collaboration), Astrophys. J. Lett. 930 (2022) L12.
  * [22] R. Abbasi et al. (IceCube Collaboration), Phys. Rev. D 106 (2022) 022005, arXiv:2111.10169.
  * [23] O. Mena, I. Mocioiu, and C. Quigg, Astropart. Phys. 28 (2007) 348, astro-ph/0610918.
  * [24] J. M. LoSecco, Universe 7 (2021) 335, arXiv:2109.01957.
  * [25] S. Chandrasekhar, The mathematical theory of black holes (Clarendon Press, Oxford, 1983).
  * [26] S. M. Bilenky, Universe 6 (2020) 134, arXiv:2008.02110.
  * [27] M. Dvornikov and A. Studenikin, J. High Energy Phys. 09 (2002) 016, hep-ph/0202113.
  * [28] M. A. Abramowicz and P. C. Fragile, Living Rev. Relativ. 16 (2013) 1, arXiv:1104.5499.
  * [29] R. Stehle and H. C. Spruit, Mon. Not. R. Astron. Soc. 304 (1999) 674.
  * [30] R. Narayan and I. Yi, Astrophys. J. Lett. 428 (1994) L13, astro-ph/9403052.
  * [31] J. Jiang, A. C. Fabian,T. Dauser, L. Gallo, J. A. Garcia, E. Kara, M. L. Parker, J. A. Tomsick, D. J. Walton, and C. S. Reynolds, Mon. Not. R. Astron. Soc. 489 (2019) 3436, arXiv:1908.07272.
  * [32] A. Sądowski, Astrophys. J. Suppl. 183 (2009) 171, arXiv:0906.0355.
  * [33] R. M. Wald, Phys. Rev. D 10 (1974) 1680.
  * [34] R. D. Blandford and D. G. Payne, Mon. Not. R. Astron. Soc. 199 (1982) 883.
  * [35] V. S. Beskin, MHD Flows in Compact Astrophysical Objects: Accretion, Winds and Jets (Springer, Heidelberg, 2010).
  * [36] N. F. Bell, V. Cirigliano, M. J. Ramsey-Musolf, P. Vogel, and M. B. Wise, Phys. Rev. Lett. 95 (2005) 151802, hep-ph/0504134.
  * [37] N. Viaux, M. Catelan, P. B. Stetson, G. G. Raffelt, J. Redondo, A. A. R. Valcarce, and A. Weiss, Astron. Astrophys. 558 (2013) A12, arXiv:1308.4627.
  * [38] S. K. Bose and M. Y. Wang, J. Math. Phys. 15 (1974) 957.
  * [39] V. Bozza, Phys. Rev. D 78 (2008) 063014, arXiv:0806.4102.
  * [40] M. B. Voloshin, M. I. Vysotskiĭ, and L. B. Okun’, Sov. Phys. JETP 64 (1986) 446,
  * [41] C. M. Will, Class. Quantum Grav. 29 (2012) 217001, arXiv:1208.3931.
  * [42] G. A. Korn and T. M. Korn, Mathematical Handbook for Scientists and Engineers: Definitions, Theorems and Formulas for Reference and Review (McGraw-Hill, New York, 1968), 2nd ed.