# Improved description of $e^{+}e^{-}\to p\bar{p}\pi^{0}$ cross section line
shape and more stringent constraints on $\psi(3770)$ and $\psi(4230)\to
p\bar{p}\pi^{0}$

Ya-Qian Wang<EMAIL_ADDRESS>Department of Physics, Hebei University,
Baoding 071002, China  Institute of High Energy Physics, Chinese Academy of
Sciences, Beijing 100049, China  Chang-Zheng Yuan<EMAIL_ADDRESS>Institute
of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China
University of Chinese Academy of Sciences, Beijing 100049, China

###### Abstract

By using the cross sections of $e^{+}e^{-}\to p\bar{p}\pi^{0}$ measured at the
center-of-mass energies from 3.65 to 4.60 GeV, we find that the line shape
well agrees with pure continuum production parametrized by a power law
function and the significance of both $\psi(3770)$ and $\psi(4230)\to
p\bar{p}\pi^{0}$ is less than $0.4\sigma$. We set more stringent constraints
on the charmless decays of the charmonium state, $\psi(3770)$, and the
charmoniumlike state, $\psi(4230)$, to $p\bar{p}\pi^{0}$ than in previous
measurements. The data are also used to estimate the cross sections of
$p\bar{p}\to\pi^{0}\psi(3770)$ and $p\bar{p}\to\pi^{0}\psi(4230)$ that are
essential for planning the data taking of the $\overline{\text{P}}\text{ANDA}$
experiment.

###### pacs:

14.20.Dh, 13.25.Gv, 13.66.Bc

The charmonium state $\psi(3770)$, lying about 40 MeV above the open charm
threshold, is expected to decay dominantly into the Okubo-Zweig-Iizuka allowed
final state, $D\bar{D}$, together with a small fraction of hadronic and
radiative transitions into final states with a lower mass charmonium state.
However, the total cross section $\sigma(e^{+}e^{-}\to\psi(3770))$ measured by
counting the total number of inclusive hadronic events is not saturated by
$D\bar{D}$ measured by counting the number of charged and neutral $D\bar{D}$
pairs c_nonDD_bes ; c_nonDD_cleo . The hadronic transitions and the radiative
transitions of the $\psi(3770)$ can be calculated reliably in the quark models
quarkmodels  and have been measured precisely in experiments c_pdg_2020 . This
indicates that noticeable charmless decays of the $\psi(3770)$, i.e.,
$\psi(3770)\to\text{\it light hadrons}$, may exist. On the other hand, people
did elaborate studies to search for exclusive charmless decays and failed to
find any mode with a large branching fraction c_nonDD_ex_BES ; c_nonDD_ex_CLEO
.

The large rate of charmless decays of the $\psi(3770)$ can also be
accommodated in theoretical models. Among many theoretical efforts trying to
solve the “$\rho\pi$ puzzle” in $J/\psi$ and $\psi(3686)$ decays observed by
the Mark-II experiment rhopi_markII , the 2S-1D charmonium mixing scenario
c_JLRosener  relates a partial width in $\psi(3770)$ decays into any final
state to the corresponding partial widths in $J/\psi$ and $\psi(3686)$ decays
c_Rosener_WYM , as both $\psi(3686)$ and $\psi(3770)$ are the mixtures of the
2S and 1D charmonium states. As a result, a large $\psi(3770)\to\text{\it
light hadrons}$ is allowed in this model, and the measurement of the
$\psi(3770)$ decays can test the model predictions. The process $\psi(3770)\to
p\bar{p}\pi^{0}$ is one of the possible exclusive channels contributing to the
charmless decays of the $\psi(3770)$.

The observation of the charmoniumlike states such as the $X(3872)$ X3872 , the
$Y(4260)$ Y4260 , and the $Z_{c}(3900)$ Zc3900  indicates that the hadrons are
more complicated than the expectations in quark models. The exotic properties
of these states such as very close to open charm thresholds and large coupling
to hidden-charm final states may suggest they are hadronic states beyond the
conventional quark model reviews . Search for their decays into light hadrons
will also shed light on their nature. Considering the small coupling to open
charm final states, we may even expect a state like the $\psi(4230)$ Y4230
(the dominant component of the $Y(4260)$ structure) has a larger decay rate to
the charmless final state than a charmonium state does. These discussions can
be extended to other vector charmoniumlike states such as the $Y(4360)$ and
$Y(4660)$ Y4660 .

If indeed $\psi(3770)\to p\bar{p}\pi^{0}$ and/or $\psi(4230)\to
p\bar{p}\pi^{0}$ are observed, one may calculate
$\sigma(p\bar{p}\to\pi^{0}\psi(3770))$ and
$\sigma(p\bar{p}\to\pi^{0}\psi(4230))$ by using the cross symmetry as has been
done in Ref. PhysRevD.73.096003 . These cross sections serve as an essential
input for the $\overline{\text{P}}\text{ANDA}$
(Anti$\underline{\text{P}}$roton $\underline{\text{An}}$nihilations experiment
at $\underline{\text{Da}}$rmstadt), which plans to study charmonium and
charmoniumlike states produced in $p\bar{p}$ annihilation.

To study the resonance contribution in the $e^{+}e^{-}\to p\bar{p}\pi^{0}$
channel, the process $\psi(3770)\to p\bar{p}\pi^{0}$ can hardly be determined
without a good knowledge of the continuum production guoyp . The continuum
cross sections of $e^{+}e^{-}\to p\bar{p}\pi^{0}$ are calculated in the energy
range from threshold up to $\sqrt{s}=4.2$ GeV, applying the conservation of
the hadron electromagnetic currents and the P-invariance of the hadron
electromagnetic interaction c_ppbpiz_th . We focus on the study of
$e^{+}e^{-}\to p\bar{p}\pi^{0}$ in this paper, and aim to extract both the
resonance and the continuum contributions.

The process of $e^{+}e^{-}\to p\bar{p}\pi^{0}$ is studied in the vicinity of
the $\psi(3770)$ by the BESIII experiment in 2014 BESIII-2014  (labeled with
“2014” hereinafter). Two indistinguishable solutions c_zhuk  for the cross
section of $\psi(3770)\to p\bar{p}\pi^{0}$ are extracted, and the maximum
cross section of $p\bar{p}\to\psi(3770)\pi^{0}$ is expected at a center-of-
mass energy (CME) of 5.26 GeV using a constant decay amplitude approximation
PhysRevD.73.096003 . The same process is studied in the energy range
$\sqrt{s}=4.008\sim 4.600$ GeV with the BESIII data taken at 13 CMEs
BESIII-2017  (labeled with “2017” in the context), and no significant
resonance is observed.

The above two analyses study the $e^{+}e^{-}\to p\bar{p}\pi^{0}$ process from
$\sqrt{s}=3.65\sim 4.60$ GeV, and in principle the data can be combined to
extract a better estimation of the physics quantities. By combining the
measurements, we benefit not only from a better description of the continuum
contribution but also from improvement on the fit result in the vicinity of
the $\psi(3770)$ and the $\psi(4230)$. In a word, besides feeding the
$\overline{\text{P}}\text{ANDA}$ experiment by estimating the cross sections
$\sigma(p\bar{p}\to\pi^{0}\psi)$, the results could also provide further
insights into the puzzling question on the mechanisms of non-$D\bar{D}$
transitions for $\psi(3770)$ and shed light on the understanding of the
$\psi(4230)$ and other vector charmoniumlike states.

Both the BESIII analyses fit the cross sections with the formula

$\displaystyle\sigma(m)=$ $\displaystyle\left|\sqrt{\sigma_{\rm
con}}+\sqrt{f_{\sigma_{\psi}}}\frac{M\Gamma}{m^{2}-M^{2}+iM\Gamma}e^{i\phi}\right|^{2},$
(1)

where $\sigma_{\rm con}=C/s^{\lambda}$ represents the continuum amplitude with
the unknown exponent $\lambda$ and a constant $C$; the resonance $\psi$ has a
fixed mass $M$ and total width $\Gamma$ c_pdg_2020 ; the factor
$f_{\sigma_{\psi}}=\sigma(e^{+}e^{-}\to\psi\to p\bar{p}\pi^{0})$ is the peak
cross section of the resonance $\psi$; the parameter $\phi$ describes the
relative phase between the continuum and resonance amplitudes.

The two analyses study the same final state but at different energy ranges.
The cross sections are fitted with the same model but with different
resonances. The parameters ($C$, $\lambda$) describing the continuum
contribution are expected to be the same. However, an obvious discrepancy
exists between the two best fits in the regions dominated by the continuum
amplitude, as shown in Fig. 1 and indicated by the very different common
parameters listed in Table 1.

Figure 1:  Extrapolation of the BESIII fit results, the 2014 BESIII-2014  and 2017 BESIII-2017  analyses give very different estimation of the continuum amplitude.  Table 1:  Fit parameters in the two BESIII analyses BESIII-2014 ; BESIII-2017 . The resonance parameters of the $\psi(3770)$ and $\psi(4230)$ are fixed in the fits. Data |  $\sqrt{s}$ (GeV) |  $C$ ($\text{ GeV}^{2\lambda}\text{pb}$) | $\lambda$ |  $M_{\psi}$ (GeV) |  $\Gamma_{\psi}$ (GeV)
---|---|---|---|---|---
2014 | $3.650-3.804$ | $(0.4\pm 0.6)\times 10^{3}$ | $1.4\pm 0.6$ | $3.77315$ | $0.0272$
2017 | $4.008-4.600$ | $(5.4\pm 5.3)\times 10^{5}$ | $4.2\pm 0.4$ | $4.251$ | $0.12$

We try to do a combined fit to the measurements presented in the two analyses
BESIII-2014 ; BESIII-2017 . Table 2 shows a compilation of the cross sections.
We do a least $\chi^{2}$ fit with

$\chi^{2}=\sum_{i=1}^{N}\frac{(\sigma^{\rm meas}_{i}-\sigma^{\rm
fit}(m_{i}))^{2}}{(\Delta\sigma^{\rm meas}_{i})^{2}},$

where $\sigma^{\rm meas}_{i}\pm\Delta\sigma^{\rm meas}_{i}$ is the dressed
cross section from experimental measurement, and $\sigma^{\rm fit}(m_{i})$ is
the cross section value calculated from the model below with the parameters
obtained from the fit. Here $m_{i}$ is the CME that corresponds to the $i$th
one of all the $N$ energy points. We only use the statistical errors in our
fits since the systematic errors ($\sim 6.5\%$) for all the data points are
correlated.

Table 2: $e^{+}e^{-}\to p\bar{p}\pi^{0}$ Born cross sections measured by the BESIII experiment BESIII-2014 ; BESIII-2017 . The first errors are statistical and second ones systematic. $\frac{1}{|1-\Pi|^{2}}$ is the vacuum polarization factor to obtain the dressed cross sections c_vp_cmd . $\sqrt{s}$ (GeV) |  $\sigma(e^{+}e^{-}\to p\bar{p}\pi^{0})$ (pb) |  $\frac{1}{|1-\Pi|^{2}}$
---|---|---
3.650 | $10.09\pm 0.84\pm 0.16$ | 1.0196
3.746 | $9.60^{+2.45}_{-2.12}\pm 0.16$ | 1.0596
3.753 | $7.28\pm 1.38\pm 0.12$ | 1.0573
3.757 | $10.44\pm 1.77\pm 0.17$ | 1.0564
3.765 | $8.73\pm 1.35\pm 0.14$ | 1.0558
3.773 | $7.71\pm 0.09\pm 0.13$ | 1.0598
3.780 | $7.92^{+2.66}_{-2.31}\pm 0.13$ | 1.0611
3.791 | $9.03\pm 1.35\pm 0.15$ | 1.0592
3.804 | $8.44\pm 1.48\pm 0.14$ | 1.0573
4.008 | $5.09\pm 0.18_{-0.24}^{+0.26}$ | 1.004
4.085 | $4.47\pm 0.46_{-0.21}^{+0.27}$ | 1.052
4.189 | $3.64\pm 0.43_{-0.19}^{+0.18}$ | 1.056
4.208 | $3.52\pm 0.39_{-0.22}^{+0.17}$ | 1.057
4.217 | $3.24\pm 0.37\pm 0.18$ | 1.057
4.226 | $3.15\pm 0.08\pm 0.14$ | 1.056
4.242 | $3.30\pm 0.36_{-0.15}^{+0.19}$ | 1.056
4.258 | $3.08\pm 0.10_{-0.15}^{+0.14}$ | 1.054
4.308 | $2.32\pm 0.33_{-0.10}^{+0.15}$ | 1.053
4.358 | $2.48\pm 0.11_{-0.12}^{+0.13}$ | 1.051
4.387 | $1.92\pm 0.26\pm 0.10$ | 1.051
4.416 | $2.16\pm 0.10_{-0.11}^{+0.10}$ | 1.053
4.600 | $1.63\pm 0.08\pm 0.08$ | 1.055

First, we fit the cross sections from 3.65 to 4.60 GeV with the continuum
amplitude only. The fitted result is shown as the red curve in Fig. 2, and the
fitted parameters are listed in Table 3. The goodness-of-fit is
$\chi^{2}/$NDF=13.0/20 (NDF is the number of degrees of freedom),
corresponding to a confidence level of 88%, a very good fit.

Figure 2: Combined fits to the $e^{+}e^{-}\to p\bar{p}\pi^{0}$ cross sections
BESIII-2014 ; BESIII-2017 . The red curve is the fit with continuum amplitude
only, the green curve is the fit with coherent sum of the continuum and the
$\psi(3770)$ amplitudes, and the dashed blue curve is the fit with coherent
sum of the continuum and the $\psi(4230)$ amplitudes.

We then fit the cross sections considering possible resonance contributions.
The continuum amplitude and the resonance ($\psi(3770)$ or $\psi(4230)$)
amplitudes are added coherently,

$\sigma(m)=\left|\sqrt{\sigma_{\rm con}}+BW(m)e^{i\phi}\right|^{2},$

where $BW(m)=\frac{\sqrt{12\pi\Gamma_{e^{+}e^{-}}\Gamma_{\rm tot}\times{\cal
B}(\psi\to p\bar{p}\pi^{0})}}{m^{2}-M^{2}+iM\Gamma_{\rm tot}}$ is a Breit-
Wigner function to describe the resonance amplitude with a mass $M$, total
width $\Gamma_{\rm tot}$, electronic partial width $\Gamma_{e^{+}e^{-}}$, and
the branching fraction ${\cal B}(\psi\to p\bar{p}\pi^{0})$. The continuum term
$\sigma_{\rm con}$ is defined the same as before. In the fits, both the mass
and the total width of the resonances are fixed according to the PDG
c_pdg_2020 , and the product $\Gamma_{e^{+}e^{-}}\times{\cal B}(\psi\to
p\bar{p}\pi^{0})$ is a free parameter. The continuum term $\sigma_{\rm con}$
and the relative phase between the continuum and the resonance amplitudes,
$\phi$, are also float parameters.

The $\psi(3770)$ or the $\psi(4230)$ resonance amplitude is added to the fit.
The fit results are shown as the solid green and dashed blue curves in Fig. 2
and listed in Table 3. The goodness-of-fit, $\chi^{2}$/NDF, is $12.3/18$ and
$12.8/18$ for $\psi(3770)$ and $\psi(4230)$, corresponding to a confidence
level of 83% and 80%, respectively. The statistical significance of the
$\psi(3770)$ is $0.36\sigma$ and that for the $\psi(4230)$ is $0.11\sigma$, by
comparing the differences in the $\chi^{2}$ and in the NDF. Two solutions are
found in each fit, one of which has a well determined magnitude and interferes
with the continuum amplitude destructively ($\phi\sim 270^{\circ}$), and the
other solution comes with large uncertainty and agrees with zero with
statistical uncertainties. The two solutions of
$\Gamma_{e^{+}e^{-}}\times{\cal B}$ may differ by several orders of magnitude
for both $\psi(3770)$ and $\psi(4230)$.

Table 3: Combined fit results for the three configurations. The subscripts of $\psi(3770)$ and $\psi(4230)$ indicate the two solutions of the fits. The errors are statistical only. There is an additional 6.5% systematic error in $\Gamma_{e^{+}e^{-}}\times{\cal B}$ and $C$ due to the common systematic error in the cross section measurements.  Fit |  $\chi^{2}$/NDF |  $\Gamma_{e^{+}e^{-}}\times{\cal B}$ (eV) |  $\phi$ (∘) |  $C$ ($10^{5}\text{GeV}^{2\lambda}$ pb) | $\lambda$
---|---|---|---|---|---
Continuum | 13.0/20 | $\ldots$ | $\dots$ | $3.07\pm 0.58$ | $3.96\pm 0.07$
$\psi(3770)_{1}$ | 12.3/18 | $(0.28\pm 1.40)\times 10^{-3}$ | $340\pm 54$ | $4.1\pm 1.8$ | $4.07\pm 0.15$
$\psi(3770)_{2}$ | 12.3/18 | $0.875\pm 0.017$ | $270.1\pm 2.4$ | $4.1\pm 1.8$ | $4.07\pm 0.15$
$\psi(4230)_{1}$ | 12.8/18 | $(0.88\pm 4.20)\times 10^{-5}$ | $78\pm 202$ | $3.22\pm 0.83$ | $3.98\pm 0.10$
$\psi(4230)_{2}$ | 12.8/18 | $0.820\pm 0.016$ | $268.7\pm 1.1$ | $3.23\pm 0.84$ | $3.98\pm 0.10$

By combining the two BESIII measurements, the fit parameters, especially those
describing the continuum are further constrained. Our first fit without any
resonance indicates in another way the extremely small significance of the
$\psi(3770)$ or $\psi(4230)$. We notice that the parameters $C$ and $\lambda$
are relatively stable in all the fits within uncertainties. The theoretical
predictions on the continuum amplitude vary significantly in both the shape
and the magnitude with different parametrizations (“new” and “old”)
c_ppbpiz_th . The new one favored by the BESIII data follows a power law with
$C=2.1\times 10^{3}~{}\text{GeV}^{2\lambda}$ and $\lambda=3.58$ in the CME
from 3.6 to 4.2 GeV c_ppbpiz_th . However, the absolute cross sections,
reflected by $C$, is still below the BESIII data by two orders of magnitude.
The $\lambda$, describing the slope of the line shape, is sensitive to the
choice of parametrization and can be further constrained by our combined fit.

Although not significant, the central value of the branching fraction of
$\psi(3770)\to p\bar{p}\pi^{0}$ is also extracted to be either $(1.1\pm 5.4\pm
0.1\pm 0.1^{\Gamma_{e^{+}e^{-}}})\times 10^{-6}$ or $(0.33\pm 0.01\pm 0.03\pm
0.03^{\Gamma_{e^{+}e^{-}}})$% using $\Gamma_{e^{+}e^{-}}=(0.262\pm 0.018)$ keV
from the PDG c_pdg_2020 , where the errors are statistical, systematic, and
from the uncertainty of quoted $\Gamma_{e^{+}e^{-}}$ c_pdg_2020 . Very
recently, authors of Ref. Achasov:2021adv  obtained
$\Gamma_{e^{+}e^{-}}=(0.19\pm 0.04)$ keV for the $\psi(3770)$ by taking into
account the contributions of the mixed $\psi(3770)$ and $\psi(3686)$
resonances. With this $\Gamma_{e^{+}e^{-}}$, the corresponding branching
fraction is calculated as $(1.5\pm 7.4\pm 0.1\pm
0.4^{\Gamma_{e^{+}e^{-}}})\times 10^{-6}$ or $(0.46\pm 0.01\pm 0.03\pm
0.10^{\Gamma_{e^{+}e^{-}}})$%. All the results are within their large
uncertainties either from our fit or from the $\Gamma_{e^{+}e^{-}}$.
Obviously, the branching fraction of $\psi(3770)\to p\bar{p}\pi^{0}$ highly
depends on the choice of the $\Gamma_{e^{+}e^{-}}$ in our calculation, and any
reliable $\Gamma_{e^{+}e^{-}}$ can be taken as input to extract the branching
fraction. The two solutions suggest either noticeable or negligible charmless
decays of $\psi(3770)\to p\bar{p}\pi^{0}$. As one of the charmless decays of
the $\psi(3770)$, the branching fraction can also be estimated considering the
mixing with the $\psi(3686)$ due to tensor forces and coupling to charmed
meson pairs c_JLRosener . This S- and D-wave charmonium mixing model accounts
for the “$\rho\pi$ puzzle” by converting expected $\psi(3686)\to\rho\pi$ to
corresponding $\psi(3770)$ partial widths. Taking the $J/\psi$ and the
$\psi(3686)$ decays into the same final state as input, the branching fraction
of $\psi(3770)\to p\bar{p}\pi^{0}$ is predicted in the range [$(2.9\pm
4.9)\times 10^{-7}$, $(15.3\pm 1.5)\times 10^{-5}$] assuming the mixing angle
$\theta=(12\pm 2)^{\circ}$. One of our fit results falls in the range of this
prediction, while the other solution is beyond the range by one order of
magnitude. If we use another mixing angle $\theta=-(27\pm 2)^{\circ}$, the
model predicts the branching fraction to be in the range [$(2.4\pm 1.9)\times
10^{-7}$, $(2.9\pm 0.3)\times 10^{-5}$]. Compared to the results with
$\theta=12^{\circ}$, the range becomes narrower, and the maximum decay rate of
$\psi(3770)\to p\bar{p}\pi^{0}$ is much smaller.

The $\overline{\text{P}}\text{ANDA}$ experiment produces neutral state with
any conventional quantum numbers ($J^{PC}$) through $p\bar{p}$ reactions.
While exotic states is also available in association with an additional meson,
e.g., $p\bar{p}\to\pi^{0}X$, where $X$ is an exotic hybrid or a charmonium
PhysRevD.73.096003 ; Pire:2013jva . The branching fraction of $\psi(3770)\to
p\bar{p}\pi^{0}$ can be taken as input for the constant decay amplitude model
PhysRevD.73.096003  to calculate the cross section of
$p\bar{p}\to\pi^{0}\psi(3770)$. The maximum cross sections predicted by the
model is $\sigma(p\bar{p}\to\pi^{0}\psi(3770))=(0.04\pm 0.20\pm 0.01\pm
0.01^{\Gamma_{e^{+}e^{-}}})$ nb or $(118\pm 6\pm 12\pm
12^{\Gamma_{e^{+}e^{-}}})$ nb at $\sqrt{s}=5.26$ GeV, both agree with previous
results BESIII-2014 . The cross sections are sensitive to the
$\Gamma_{e^{+}e^{-}}$ in the same way as the branching fraction of
$\psi(3770)\to p\bar{p}\pi^{0}$. If we use $\Gamma_{e^{+}e^{-}}=(0.19\pm
0.04)$ keV Achasov:2021adv , the above cross section becomes $(0.05\pm 0.27\pm
0.01\pm 0.02^{\Gamma_{e^{+}e^{-}}})$ nb or $(164\pm 7\pm 13\pm
37^{\Gamma_{e^{+}e^{-}}})$ nb, which are larger by 25% and 39% compared with
the results with $\Gamma_{e^{+}e^{-}}$ from the PDG c_pdg_2020 . The
production rate of $p\bar{p}\to\pi^{0}\psi(4230)$ reaches the peak at
$\sqrt{s}=5.9$ GeV, which exceeds the energy of
$\overline{\text{P}}\text{ANDA}$ by 0.4 GeV. The cross sections at
$\sqrt{s}=5.5$ GeV is predicted as
$\sigma(p\bar{p}\to\pi^{0}\psi(4230))=\frac{1}{\Gamma_{e^{+}e^{-}}}(29.4\pm
0.59\pm 5.3)\text{~{}nb}\cdot\text{keV}$ or
$\frac{1}{\Gamma_{e^{+}e^{-}}}(0.32\pm 1.6\pm
0.06)\text{~{}nb}\cdot\text{eV}$, where the currently unknown electron partial
width of the $\psi(4230)$ is expected in the future.

In summary, we do a combined fit to the $e^{+}e^{-}\to p\bar{p}\pi^{0}$ cross
sections measured by the BESIII experiment and obtain a better description of
the continuum production of $e^{+}e^{-}\to p\bar{p}\pi^{0}$. It shows that the
interference effect between continuum and the resonance amplitudes plays a
very important role in determining the resonance contributions in this mode.
BESIII has collected data at much more energy points from 4.130 to 4.946 GeV,
and the cross sections are expected to be further measured and better examined
in the future white .

This work is supported in part by National Key Research and Development
Program of China under Contract No. 2020YFA0406301, Advanced Talents
Incubation Program of the Hebei University under Contract No. 521100221013,
and National Natural Science Foundation of China (NSFC) under Contract Nos.
11961141012, 11835012, and 11521505.

## References

  * (1) M. Ablikim et al. [BES], Phys. Rev. Lett. 97, 121801 (2006); Phys. Lett. B 641, 145-155 (2006); Phys. Lett. B 652, 238-244 (2007); Phys. Rev. D 76, 122002 (2007); Phys. Lett. B 659, 74-79 (2008).
  * (2) D. Besson et al. [CLEO], Phys. Rev. Lett. 96, 092002 (2006) [erratum: Phys. Rev. Lett. 104, 159901 (2010)].
  * (3) S. Godfrey and N. Isgur, Phys. Rev. D 32, 189-231 (1985).
  * (4) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020).
  * (5) G. S. Huang et al. [CLEO], Phys. Rev. Lett. 96, 032003 (2006).
  * (6) M. Ablikim et al. [BES], Phys. Lett. B 650, 111-118 (2007).
  * (7) M. E. B. Franklin, G. J. Feldman, G. S. Abrams, M. S. Alam, C. A. Blocker, A. Blondel, A. Boyarski, M. Breidenbach, D. L. Burke and W. C. Carithers, et al. Phys. Rev. Lett. 51, 963-966 (1983).
  * (8) J. L. Rosner, Phys. Rev. D 64, 094002 (2001).
  * (9) P. Wang, C. Z. Yuan and X. H. Mo, Phys. Rev. D 70, 114014 (2004).
  * (10) S. K. Choi et al. [Belle], Phys. Rev. Lett. 91, 262001 (2003).
  * (11) B. Aubert et al. [BaBar], Phys. Rev. Lett. 95, 142001 (2005).
  * (12) M. Ablikim et al. [BESIII], Phys. Rev. Lett. 110, 252001 (2013); Z. Q. Liu et al. [Belle], Phys. Rev. Lett. 110, 252002 (2013) [erratum: Phys. Rev. Lett. 111, 019901 (2013)].
  * (13) For recent reviews, see N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873, 1-154 (2020); F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao and B. S. Zou, Rev. Mod. Phys. 90, no. 1, 015004 (2018); H. X. Chen, W. Chen, X. Liu and S. L. Zhu, Phys. Rept. 639, 1 (2016).
  * (14) M. Ablikim et al. [BESIII], Phys. Rev. Lett. 118, no.9, 092001 (2017).
  * (15) X. L. Wang et al. [Belle], Phys. Rev. Lett. 99, 142002 (2007).
  * (16) A. Lundborg, T. Barnes and U. Wiedner, Phys. Rev. D 73, 096003 (2006).
  * (17) Y. P. Guo and C. Z. Yuan, Phys. Rev. D 105, no.11, 114001 (2022).
  * (18) G. I. Gakh, M. I. Konchatnij, N. P. Merenkov and E. Tomasi-Gustafsson, Phys. Rev. D 105, no.9, 094029 (2022).
  * (19) M. Ablikim et al. [BESIII], Phys. Rev. D 90, no.3, 032007 (2014).
  * (20) K. Zhu, X. H. Mo, C. Z. Yuan and P. Wang, Int. J. Mod. Phys. A 26, 4511 (2011).
  * (21) M. Ablikim et al. [BESIII], Phys. Lett. B 771, 45-51 (2017).
  * (22) F. Ignatov, http://cmd.inp.nsk.su/~ignatov/vpl/.
  * (23) N. N. Achasov and G. N. Shestakov, Phys. Rev. D 103, no.7, 076017 (2021).
  * (24) B. Pire, K. Semenov-Tian-Shansky and L. Szymanowski, Phys. Lett. B 724, 99-107 (2013) [erratum: Phys. Lett. B 764, 335-335 (2017)].
  * (25) M. Ablikim et al. [BESIII], Chin. Phys. C 44, no.4, 040001 (2020).


# ALIFE2024 template

Michael Crosscombe∗, Ilya Horiguchi, Norihiro Maruyama, Shigeto Dobata
Takashi Ikegami

Department of General Systems Studies, Graduate School of Arts and Sciences,
The University of Tokyo, Japan
<EMAIL_ADDRESS>
# A Simulation Environment for the Neuroevolution of
Ant Colony Dynamics

Michael Crosscombe∗, Ilya Horiguchi, Norihiro Maruyama, Shigeto Dobata
Takashi Ikegami

Department of General Systems Studies, Graduate School of Arts and Sciences,
The University of Tokyo, Japan
<EMAIL_ADDRESS>
###### Abstract

We introduce a simulation environment to facilitate research into emergent
collective behaviour, with a focus on replicating the dynamics of ant
colonies. By leveraging real-world data, the environment simulates a target
ant trail that a controllable agent must learn to replicate, using sensory
data observed by the target ant. This work aims to contribute to the
neuroevolution of models for collective behaviour, focusing on evolving neural
architectures that encode domain-specific behaviours in the network topology.
By evolving models that can be modified and studied in a controlled
environment, we can uncover the necessary conditions required for collective
behaviours to emerge. We hope this environment will be useful to those
studying the role of interactions in emergent behaviour within collective
systems.

## Introduction

In the pursuit of novel forms of Collective Intelligence (CI) living systems
provide us with extensive datasets for understanding the fundamental
properties of collective behaviour. That is, living systems are useful for
studying how local interactions between individuals lead to the emergence of
global behaviours in the population through a process of self-organisation
(Couzin and Krause,, 2003). However, a full understanding of the necessary
conditions for collective behaviour cannot come from merely observing living
systems. Instead, we must study said conditions through controlled experiments
that, while inspired by natural systems, allow for the evolution and
manipulation of models and their variables to observe the effects.

Simulations afford the opportunity to model individual behaviours and
interactions within a controlled setting, enabling the study of emergent
properties in a systematic and repeatable manner. Recently, efforts to merge
deep learning and complex systems has greatly improved CI research
capabilities by enabling simulations to scale to hundreds or thousands of
individuals; many more than was previously possible Ha and Tang, (2022). By
leveraging real-world data from living systems, we can simulate environments
with the same conditions in which collective behaviours naturally occur,
allowing for the exploration of factors such as communication, information
aggregation (Kameda et al.,, 2022), coarse-graining (Flack,, 2017), and their
impacts on the emergence of collective behaviours.

Here we introduce a Gymnasium (Towers et al.,, 2023) environment to facilitate
research into the necessary conditions for emergent collective
behaviour111https://github.com/TooHuman/AntDynamics-Gymnasium. Our environment
presents a target ant trial (from real data) that the learning agent attempts
to replicate, mimicking the dynamics based on relative sensory data observed
by the target ant along the trail (e.g., nearby ants).

## Simulating Ant Dynamics

Figure 1:  Left: Top-down video recordings, cropped and scaled to
$1280$x$1280$ resolution. The ants are filmed from above in an evenly-lit
$100$ mm diameter arena. Individual ant positions are extracted from a $4$
hour recording. Right: Our simulation environment using Pygame and Gymnasium
to reproduce ant colony dynamics with an agent that can interface with a
policy network.

We chose to reproduce the dynamics of the ant species Pristomyrmex punctatus
which is characterised by its lack of a queen caste, with mature colonies
often consisting of thousands of female workers. Instead of constructing
nests, the colony exhibits a strong clustering behaviour (Tsuji and Dobata,,
2011). (For further information on the experimental setup and the
detection/tracking pipeline, see Maruyama et al., (2023).) Based on data
extracted from video recordings, we built a simulation environment in which
the agent (ant) begins each trial at position $P_{\alpha}(0)=P_{\tau}(0)$, the
beginning of a (real) target ant trail. Figure 1 shows the comparison between
the video recordings and the simulation environment. Currently, we select a
target ant at random from the data whilst ensuring that the target moved a
modest amount during the $t_{\text{lim}}=30$ second time period. The target
ant and corresponding trail are depicted in red, while the agent and agent
trail are shown in blue. This is seen more clearly in Figure 2 where we also
depict the ant’s vision which we further separate into four segments: $5$ in
the forward-facing direction (indicated by the blue line depicting the agent’s
direction $\theta$), and $3$ facing backwards. Each ant in the environment
begins at some position translated to pixel coordinates $(x,y)$ from the real
data, relative to the arena’s dimensions and scale.

Figure 2:  Close-up view of the controllable agent (blue) and the agent’s
corresponding trail $P_{\alpha}$(light blue) next to the target trail (red).
The blue circle surrounding the agent represents the agent’s vision divided
into segments.

Given two paths, the agent’s path $P_{\alpha}(x)$ and the target path
$P_{\tau}(x)$, the reward function $\mathcal{R}$ is defined as the scaled
negative area between these paths over time, given by:

$R=-\sum_{t=1}^{t_{\text{lim}}}\left(1-\frac{A_{t}}{\sqrt{1+A_{t}^{2}}}\right)$
(1)

where $A_{t}$ represents the area of the trapezoid formed between the points
along both trails at time $t$ and $t-1$. Time $t$ represents a discrete time
step in our simulation and reflects the real-time data underlying the ant
movements. Each trial runs while $t\leq t_{\text{lim}}$. This reward function
is designed such that perfect alignment of the two trails yields a reward of
$0$, with the reward becoming increasingly negative as their deviation
increases. However, alternative reward functions ought to be designed for more
effective learning.

There are currently $13$ inputs for the agent observations:
$\\{x,y,s,\theta,\dot{\theta},V_{fl_{1}},V_{fl_{2}},V_{fc},V_{fr_{2}},V_{fr_{1}},V_{r},V_{b},V_{l}\\}$
denoting the agent’s position, speed, angle, rotation speed and visual
observations. Each visual input provides a proportional measure of the number
of ants detectable within the range and angle of vision. The outputs of the
policy network are simply: forward, backward, turn-left, turn-right. (Due to
our inability to detect pheromone trails from video recordings, these have
been omitted.)

## Encoding Collective Behaviours

We intend to use this environment for the neuroevolution of models which
accurately capture the natural dynamics of ant colonies. The process of
neuroevolution allows us to investigate how collective behaviours (both
existing and novel) might be evolved through a series of gradual changes to
the neural network that controls an individual such that, when replicated in a
collective of individuals, intelligent collective behaviour emerges from their
interactions. More specifically, rather than focusing on adapting the weights
of a network for a given task, we intend to focus on constructing the network
topology itself, with the goal of encoding the task into the structure of the
network to make inferring such a model more efficient, such as in the works of
Gaier and Ha, (2019) and Stanley and Miikkulainen, (2002). An alternative
approach to growing neural networks for specific tasks was recently introduced
by Najarro et al., (2023). Our purpose in evolving such policy networks is to
understand the role that interaction bottlenecks within and between
individuals play in the formation of their collective behaviours (Crosscombe
and Sato,, 2023), so that we might better develop novel collective behaviours
in the future.

Our environment is designed to not only replicate the intricate social
interactions among ants but also to understand the emergent properties of
their collective intelligence. By providing a simulation environment that
captures real-time observational data from the target domain, researchers can
iteratively evolve neural network configurations, promoting the emergence of
sophisticated behavioural patterns akin to those found in nature.

## Conclusions

We have introduced a new Gymnasium-based environment which features the
reproduction of real ant colony behaviours that can be used to study the
facets of emergent collective behaviour. Eventually, we aim to encode a target
collective behaviour into neural networks through a process of neuroevolution,
with emphasis on the evolved topology as introduced by Gaier and Ha, (2019).
The evolved network can then be used as a policy network to control multiple
agents which exhibit realistic collective behaviours inline with those
observed in colonies of Pristomyrmex punctatus. We hope that this work
contributes to research into the necessary conditions for such collective
behaviours to emerge.

## Acknowledgements

Michael Crosscombe is supported by a Grant-in-Aid for JSPS Fellows 23KF0108.
This work was partially supported by JSPS KAKENHI grants 21H04885 and
24H00707.

## References

  * Couzin and Krause, (2003) Couzin, I. and Krause, J. (2003).  Self-organization and collective behavior in vertebrates.  Advances in The Study of Behavior, 32(1):1–75.
  * Crosscombe and Sato, (2023) Crosscombe, M. and Sato, H. (2023).  On the existence of information bottlenecks in living and non-living systems.  arXiv [physics.soc-ph].
  * Flack, (2017) Flack, J. C. (2017).  Coarse-graining as a downward causation mechanism.  Philos. Trans. A Math. Phys. Eng. Sci., 375(2109).
  * Gaier and Ha, (2019) Gaier, A. and Ha, D. (2019).  Weight agnostic neural networks.  arXiv [cs.LG].
  * Ha and Tang, (2022) Ha, D. and Tang, Y. (2022).  Collective intelligence for deep learning: A survey of recent developments.  Collective Intelligence, 1(1):263391372211148.
  * Kameda et al., (2022) Kameda, T., Toyokawa, W., and Tindale, R. S. (2022).  Information aggregation and collective intelligence beyond the wisdom of crowds.  Nature Reviews Psychology, 1(6):345–357.
  * Maruyama et al., (2023) Maruyama, N., Crosscombe, M., Dobata, S., and Ikegami, T. (2023).  Emergence of differentiation of deterministic/stochastic behavior in ants’ collectives.  In The 2023 Conference on Artificial Life, page 116. MIT Press.
  * Najarro et al., (2023) Najarro, E., Sudhakaran, S., and Risi, S. (2023).  Towards self-assembling artificial neural networks through neural developmental programs.  In The 2023 Conference on Artificial Life. MIT Press.
  * Stanley and Miikkulainen, (2002) Stanley, K. O. and Miikkulainen, R. (2002).  Evolving neural networks through augmenting topologies.  Evol. Comput., 10(2):99–127.
  * Towers et al., (2023) Towers, M., Terry, J. K., Kwiatkowski, A., Balis, J. U., Cola, G. d., Deleu, T., Goulão, M., Kallinteris, A., KG, A., Krimmel, M., Perez-Vicente, R., Pierré, A., Schulhoff, S., Tai, J. J., Shen, A. T. J., and Younis, O. G. (2023).  Gymnasium.
  * Tsuji and Dobata, (2011) Tsuji, K. and Dobata, S. (2011).  Social cancer and the biology of the clonal ant pristomyrmex punctatus (hymenoptera: Formicidae).  Myrmecological News, 15:91–99.


11institutetext: Tuorla Observatory, Department of Physics and Astronomy,
FI-20014 University of Turku, Finland 11email<EMAIL_ADDRESS>22institutetext: Nordita, KTH Royal Institute of Technology and Stockholm
University, Hannes Alfvéns väg 12, SE-10691 Stockholm, Sweden 33institutetext:
Physics Department and Columbia Astrophysics Laboratory, Columbia University,
538 West 120th Street New York, NY 10027 44institutetext: Center for
Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY
10010, USA 55institutetext: Institut für Astronomie und Astrophysik,
Universität Tübingen, Sand 1, D-72076 Tübingen, Germany

#  artpol: Analytical ray-tracing method for spectro-polarimetric properties
of accretion disks around Kerr black holes

Vladislav Loktev 11 Alexandra Veledina 1122 Juri Poutanen 11 Joonas Nättilä
3344 Valery F. Suleimanov 55

Spectro-polarimetric signatures of accretion disks in X-ray binaries and
active galactic nuclei contain information about the masses and spins of their
central black holes, as well as the geometry of matter close to the compact
objects. This information can be extracted using the means of X-ray
polarimetry. In this work, we present a fast analytical ray-tracing technique
for polarized light artpol that helps obtain the spinning black hole
parameters from the observed properties. This technique can replace the
otherwise time-consuming numerical ray-tracing calculations. We show that
artpol proves accurate for Kerr black holes with dimensionless spin parameter
$a\leq 0.94$ while being over four orders of magnitude faster than direct ray-
tracing calculations. This approach opens broad prospects for directly fitting
the spectro-polarimetric data from the Imaging X-ray Polarimetry Explorer.

###### Key Words.:

accretion, accretion disks – galaxies: active – gravitational lensing: strong
– methods: analytical – polarization – stars: black holes

##  1 Introduction

The spin of a black hole (BH) is a fundamental parameter that controls the
behavior of the inflowing matter, the accretion disk, and the properties of
the outflowing material, the relativistic jets. The magnitude of spin
determines the curvature of space-time close to the BH, determines the energy
dissipation profile, and affects the spectral energy distribution of the
observed emission. The spin values for both Galactic BHs in X-ray binaries and
their supermassive counterparts have been probed using the distinct imprints
in spectral and timing properties via the continuum-fitting method and the
iron line reflection/reverberation method (Miller et al. 2009; McClintock et
al. 2014; Uttley et al. 2014). The methods are based on obtaining the radius
of the innermost stable circular orbit (ISCO), which, in turn, is related to
the BH spin. For X-ray binaries, an additional constraint on the parameters
comes from the relativistic precession model (Motta et al. 2014, 2022), which
links the mass, spin, and radius with the characteristic frequencies found in
the X-ray light curves. The statistical distribution of spins probed by these
methods, as well as spin values obtained by different methods for the same
source, do not match (Draghis et al. 2023). This calls for an alternative
method to verify spin determination measurements.

Polarization of radiation escaping from the BH vicinity can be used as a fine
tool to determine the curvature of space-time. In this context, X-ray
polarimetric signatures of accretion disks have long been anticipated to carry
important information about the BH spin, and the launch of Imaging X-ray
Polarimetry Explorer (IXPE, Weisskopf et al. 2022) has opened these new
exciting possibilities. At the same time, the observed polarimetric signatures
of BH X-ray binaries have proved many previously proposed models to fail in
describing the data (Krawczynski et al. 2022; Rodriguez Cavero et al. 2023;
Ratheesh et al. 2023), requiring to alter basic assumptions on the geometry
and radiative mechanisms producing local spectra, which are often impossible
to alter in the data-fitting models. A fast tool that relates the local
spectra with the observed spectro-polarimetric signatures is required.

The accretion disk polarization is produced by multiple scatterings in the
upper layers of its atmosphere. The first predictions of the disk polarization
(Rees 1975) have been made using the results of calculations in the case of
pure electron-scattering, semi-infinite plane-parallel atmospheres
(Chandrasekhar & Breen 1947; Sobolev 1949; Chandrasekhar 1960; Sobolev 1963),
and were limited to Newtonian approximation. The polarization degree (PD) may
be altered due to the presence of absorption effects in the atmosphere, which
were considered using Monte-Carlo (Lightman & Shapiro 1975) and analytical
(Loskutov & Sobolev 1979, 1981) means. Further, the effects of general and
special relativity introduce important modifications to both PD and
polarization angle (PA; Connors & Stark 1977; Stark & Connors 1977; Pineault &
Roeder 1977a, b). Aberration and light deflection lead to a rotation of the PA
and alter the viewing angle of different parts of the accretion disk,
affecting the PD. Additional rotation of PA along the photon trajectory is
expected for the case of spinning BHs, described by the Kerr (1963) metric,
thanks to the frame-dragging effects. Therefore, the observed spectral
dependence of polarimetric signatures can act as an independent probe of the
BH spin (Connors et al. 1980; Dovčiak et al. 2008).

Precise computations of the effects of general and special relativity on the
polarization properties can be done using the parallel transport of the
polarization vector along null geodesics, which often involves computations of
the Walker & Penrose (1970) constant of motion (Connors et al. 1980; Dovčiak
et al. 2008; Ingram et al. 2015). The geodesics, in turn, are computed using
the ray-tracing techniques (e.g., Dexter 2016; Zhang et al. 2019; Chan et al.
2013; Pihajoki et al. 2018). It is also essential to keep track of the
convergence of the computed flux, PD and PA down to the characteristic scale
of the observed errors on these quantities–this often means that several
simulations with increased resolution must be performed for each parameter
set. An additional source of computational errors comes from the assumptions
of a small outer radius of the disk, which is enforced by the high
computational costs. Typically, the value $\lesssim 100R_{\rm S}$ (where
$R_{\rm S}$ is the Schwarzschild radius) is considered, which allows for the
(highly polarized) secondary images of the disk to be visible in the region
surrounding the outer radius of the simulated disk – in reality, a much higher
extend of the disk, $\gtrsim 10^{5}R_{\rm S}$, completely covers those from
the line of sight.

The ray-tracing technique is too computationally expensive for a direct data
fitting. Instead, pre-computed geodesics have been used to accelerate
calculations (e.g., Li et al. 2005; Krawczynski 2012; Krawczynski &
Beheshtipour 2022). Alternatively, analytical and semi-analytical approaches
can be applied to solve geodesic equations (Dexter & Agol 2009; Yang & Wang
2014; Cárdenas-Avendaño et al. 2023). Their applicability is, however,
limited, e.g., the semi-analytical expressions for geodesics can only be used
for the equatorial plane of the BH. Finally, in the Schwarzschild metric, the
calculation of geodesics may be omitted, and the photon trajectory can be
treated using an approximation to the light bending relation (Beloborodov
2002; Poutanen 2020a). The latter approach gave reliable results for the low-
energy synchrotron emission observed from the supermassive BH M87*, with an
almost face-on disk (Narayan et al. 2021).

Physical understanding of the modifications to polarization signatures caused
by the curved space-time and fast motions of matter is difficult to achieve if
we use the implicit Walker-Penrose constant. For this, one can use the
explicit analytical expression for the rotation of the polarization plane
along the photon path, which we call the analytical ray-tracing technique for
polarized light (artpol hereafter). This approach has been used to extract
polarimetric properties of spinning spherical (Poutanen 2020b) and oblate
(Loktev et al. 2020) neutron stars, as well as accretion disks around
Schwarzschild BHs (Loktev et al. 2022). The method was first proposed in
Pineault (1977) and applied, in the context of accretion disks, to the PA
rotation caused by, separately, general and special relativity, while the
expression for their combined effects was first derived in Loktev et al.
(2022).

In this work, we apply artpol to the accretion disks around spinning BH. In
Sect. 2, we describe the formalism that can be used in spectro-polarimetric
modeling and imaging of accretion disks around Kerr BHs. In Sect. 3, we
compare the results of the artpol technique to those obtained using explicit
ray-tracing calculations. We show that the PD and PA computed using this
method remain accurate to the level of current observational (IXPE)
uncertainties for BH spin parameters up to $a=0.94$. We summarize our findings
and discuss a broad range of applications of the technique in Sect. 4.

##  2 Methods

###  2.1 Local model of the disk emission

We consider a BH with mass $M$ and dimensionless spin $a=Jc/GM^{2}$, where $J$
is the angular momentum. For the sake of simplicity, in this paper, we
consider a standard equatorial geometrically thin accretion disk (Novikov &
Thorne 1973). The BH is situated at the origin with the spin directed along
the $z$-axis and is orthogonal to the disk plane. We assume that the state of
the fluid only depends on the Boyer-Lindquist dimensionless radial coordinate
$r=R/{R}_{\rm S}$, which expresses the distance from the central object $R$ in
units of the BH Schwarzschild radius ${R}_{\rm S}=2GM/c^{2}$. At a given
radius $r$, the matter moves, relative to the locally non-rotating observer,
with Keplerian velocity, which we express in the units of the speed of light
$c$ as (Kato et al. 2008):

$\beta=\frac{\cal F}{\cal B\sqrt{\cal D}}\sqrt{\frac{1}{2r}},$ (1)

where

$\displaystyle\cal B$ $\displaystyle=1+\frac{a}{\sqrt{8r^{3}}},$ (2)
$\displaystyle\cal D$ $\displaystyle=1-\frac{1}{r}+\frac{a^{2}}{4r^{2}},$ (3)
$\displaystyle{\cal F}$
$\displaystyle=1-\frac{a}{\sqrt{2r^{3}}}+\frac{a^{2}}{4r^{3}}.$ (4)

The disk is considered to be optically thick. The energy flux from a surface
element is described by the effective temperature, which has the radial
profile in the form

$T_{\rm eff}^{4}(r)=\frac{3GM\dot{M}}{8\pi\sigma_{\rm
SB}R^{3}}f(r,a)=T_{*}^{4}\frac{f(r,a)}{r^{3}},$ (5)

where

$T_{*}^{4}=\frac{3GM\dot{M}}{8\pi\sigma_{\rm SB}{R}_{\rm S}^{3}},$ (6)

combines the BH mass and the accretion rate $\dot{M}$, and $f(r,a)$ is the
factor accounting for the relativistic and boundary condition corrections
(Page & Thorne 1974; Eq. (3.191) in Kato et al. 2008). In the Newtonian case,
assuming the inner disk edge at $r=3$, the correction factor is

$f(r)=1-\sqrt{\frac{3}{r}}.$ (7)

The disk extends from the ISCO up to an outer edge at $r_{\rm{out}}=3000$,
which was chosen to reproduce the correct spectra down to the photon energy
$E\sim 0.01kT_{*}$. The radius of the ISCO depends on the spin as

${r}_{\rm{ISCO}}=\frac{1}{2}\left(3+Z_{2}\pm\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})}\right),$
(8)

where the minus sign corresponds to a corotation of the disk and the BH and
the plus sign is for the case of the retrograde rotation, and

$\displaystyle Z_{1}$ $\displaystyle=1+\sqrt[\scriptstyle
3]{1-a^{2}}\left(\sqrt[\scriptstyle 3]{1+a}+\sqrt[\scriptstyle
3]{1-a}\right),$ (9) $\displaystyle Z_{2}$
$\displaystyle=\sqrt{3a^{2}+Z_{1}^{2}}.$ (10)

The local spectrum of the disk is assumed to be a diluted blackbody with the
color correction $f_{\rm col}=1.7$ (Shimura & Takahara 1995) with the flux
being

$F_{E^{\prime}}=\frac{\pi}{f_{\rm col}^{4}}B_{E^{\prime}}(f_{\rm col}T_{\rm
eff}),$ (11)

where ${E^{\prime}}$ is photon energy measured in the frame comoving with the
matter of the disk (where all quantities are denoted with primes) and
$B_{E^{\prime}}$ is the Planck function. The angular distribution of the
specific intensity and the polarization is assumed to correspond to the case
of electron-scattering-dominated semi-infinite atmosphere (Chandrasekhar &
Breen 1947; Chandrasekhar 1960; Sobolev 1949, 1963). The three-component
Stokes vector fully describes a linearly polarized radiation field as a
function of photon energy $E^{\prime}$ and zenith angle $\zeta^{\prime}$
(measured in the comoving frame):

$\bm{I}^{\prime}_{E^{\prime}}(\zeta^{\prime})=\begin{bmatrix}I^{\prime}_{E^{\prime}}\\\
Q^{\prime}_{E^{\prime}}\\\
U^{\prime}_{E^{\prime}}\end{bmatrix}=\frac{1}{{\pi}}F_{E^{\prime}}a_{\rm
es}(\zeta^{\prime})\begin{bmatrix}1\\\ p_{\rm es}(\zeta^{\prime})\cos
2\chi_{0}\\\ p_{\rm es}(\zeta^{\prime})\sin 2\chi_{0}\end{bmatrix},$ (12)

where the angular distribution can be approximated as (Nättilä & Pihajoki
2018)

$a_{\rm
es}(\zeta^{\prime})=\frac{60}{143}(1+2.3\cos{\zeta^{\prime}}-0.3\cos^{2}{\zeta^{\prime}})$
(13)

and the PD (Viironen & Poutanen 2004)

$p_{\rm
es}(\zeta^{\prime})=0.1171\frac{1-\cos{\zeta^{\prime}}}{1+3.582\cos{\zeta^{\prime}}}.$
(14)

The Stokes parameter $U$ is zero in the comoving frame because of the
azimuthal symmetry. Furthermore, in our case, the escaping radiation is
polarized perpendicular to the meridional plane formed by the normal to the
local surface and the photon momentum, resulting in the local PA of
$\chi_{0}=\pi/2$ (Stokes $Q<0$).

Once the model for the disk structure (i.e., radial dependence of the velocity
profile) and the energy and angular dependence of the local Stokes vector are
established, we need to consider how the polarized radiation is modified
toward the observer in the curved space-time. One way to account for that is
to compute a set of geodesics from the vicinity of the BH to the observer (we
focus on this in more detail in Sect. 2.3). The other possibility is to use a
faster approach that exploits analytical approximations, which we describe in
Sect. 2.2.

Figure 1: Geometry of a flat accretion disk ring. An element of the flat
equatorial disk, defined by the unit radius-vector $\hat{r}$ at azimuthal
angle $\varphi$, moves along the direction $\hat{v}$. The coordinate system is
based on the normal of the disk $\hat{n}$ and the observer vector $\hat{o}$. A
photon is emitted from the disk element along $\mbox{\boldmath$\hat{k}$}_{0}$.
The important angles used in this work are also shown.

###  2.2 Analytical ray-tracing in Schwarzschild metric

We previously tested the analytical method of calculating the rotation of the
PA for the case of the planar light trajectories (Loktev et al. 2022). This
assumption is fulfilled in the Schwarzschild metric but does not hold for the
Kerr metric. Nevertheless, we can use this approach as the first-order
approximation for the case of a spinning BH.

The considered geometry is shown in Fig. 1. The disk is assumed to be flat and
located in the equatorial plane of the spinning BH. An element of the
accretion disk surface, located at the tip of the radius vector $\bm{r}$, is
described in Boyer-Lindquist coordinates by an azimuth angle $\varphi$ and the
length $r$. We utilize a 3-dimensional Cartesian coordinate system to describe
the polarization frame’s rotation. The normal to the disk $\hat{n}$ is aligned
with the $z$ axis and the direction to the observer $\hat{o}$ is in the
$x$-$z$ plane making an inclination $i$ to the normal:

$\hat{n}$ $\displaystyle=$ $\displaystyle(0,0,1),$ (15) $\hat{o}$
$\displaystyle=$ $\displaystyle(\sin i,0,\cos i).$ (16)

In these coordinates, the unit vector of the disk element

$\mbox{\boldmath$\hat{r}$}=(\cos\varphi,\sin\varphi,0)$ (17)

makes an angle $\psi$ with the the observer direction

$\cos\psi=\mbox{\boldmath$\hat{r}$}\cdot\mbox{\boldmath$\hat{o}$}=\sin
i\cos\phi.$ (18)

Close to the disk surface, photon escapes along a unit vector

$\mbox{\boldmath$\hat{k}$}_{0}=[\sin\alpha\
\mbox{\boldmath$\hat{o}$}+\sin(\psi-\alpha)\
\mbox{\boldmath$\hat{r}$}]/\sin\psi,$ (19)

where $\alpha$ is the angle between the radius vector and the observer
direction,

$\cos\alpha=\mbox{\boldmath$\hat{r}$}\cdot\mbox{\boldmath$\hat{k}$}_{0},$ (20)

and is related to $\psi$ through the light bending formula (see below). The
photon momentum makes an angle $\zeta$ with the normal

$\cos\zeta=\mbox{\boldmath$\hat{n}$}\cdot\mbox{\boldmath$\hat{k}$}_{0}=\frac{\sin\alpha}{\sin\psi}\cos
i.$ (21)

Here we assume that the fluid velocity is aligned with the azimuthal unit
vector

$\mbox{\boldmath$\hat{\varphi}$}=(-\sin\varphi,\cos\varphi,0)=\mbox{\boldmath$\hat{v}$}.$
(22)

The matter moves with Keplerian speed given by Eq. (1).

The unit vector of photon momentum in the laboratory frame
$\mbox{\boldmath$\hat{k}$}_{0}$ and the one in the comoving frame
$\mbox{\boldmath$\hat{k}$}^{\prime}_{0}$ are related by the Lorentz
transformation:

$\mbox{\boldmath$\hat{k}$}^{\prime}_{0}=\delta\left[\mbox{\boldmath$\hat{k}$}_{0}-\gamma\beta\mbox{\boldmath$\hat{v}$}+(\gamma-1)\mbox{\boldmath$\hat{v}$}(\mbox{\boldmath$\hat{v}$}\cdot\mbox{\boldmath$\hat{k}$}_{0})\right],$
(23)

where $\gamma=1/\sqrt{1-\beta^{2}}$ is the Lorentz factor and

$\delta=\frac{1}{\gamma(1-\beta\cos\xi)}=\gamma(1+\beta\cos\xi^{\prime})$ (24)

is the Doppler factor and

$\displaystyle\cos\xi$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\hat{v}$}\cdot\mbox{\boldmath$\hat{k}$}_{0}=-\frac{\sin\alpha}{\sin\psi}\sin
i\sin\phi,$ (25) $\displaystyle\cos\xi^{\prime}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\hat{v}$}\cdot\mbox{\boldmath$\hat{k}$}^{\prime}_{0}=\frac{\cos\xi-\beta}{1-\beta\cos\xi}.$
(26)

We then also get

$\cos\zeta^{\prime}=\mbox{\boldmath$\hat{n}$}\cdot\mbox{\boldmath$\hat{k}$}^{\prime}_{0}=\delta\cos\zeta.$
(27)

The geodesics are not computed explicitly – we only need the relation between
the angles $\alpha$ and $\psi$, for which we use the following analytical
approximate formula (Poutanen 2020a):

$\cos\alpha=1-(1-u)y\left\\{1+\frac{u^{2}y^{2}}{112}-\frac{{\rm
e}}{100}uy\left[\ln\left(1-\frac{y}{2}\right)+\frac{y}{2}\right]\right\\},$
(28)

where $y=1-\cos\psi$ and $u=1/r$. The formula is not applicable for $r<1$
(within the Schwarzschild radius), hence for the cases with
${r}_{\rm{ISCO}}<1$, i.e., $a\gtrsim 0.943$, artpol cannot be used.

Under the Schwarzschild metric assumption, the geodesics are flat; therefore,
the parallel transport of the Lorentz frame along the geodesics is
unnecessary. The rotation of the polarization basis is a sum of several simple
rotations due to the gravitational light bending (the general relativity
effect) and Lorentz aberration (the special relativity effect),

$\chi^{\text{tot}}=\chi^{\text{GR}}+\chi^{\text{SR}}.$ (29)

Analytical expressions for those have been derived in Loktev et al. (2022):

$\tan\chi^{\text{GR}}=\frac{\cos i\sin\varphi}{\cos\varphi+\tilde{a}\sin i},$
(30)

where $\tilde{a}=(1-\cos\alpha\cos\psi)/(\cos\alpha-\cos\psi)$, and

$\tan\chi^{\text{SR}}=-\beta\
\frac{\cos\alpha\cos\zeta}{\sin^{2}\zeta-\beta\cos\xi}.$ (31)

Expressions for $\chi^{\text{GR}}$ and $\chi^{\text{SR}}$ in vector form can
be found in Loktev et al. (2022).

The PA is measured in the main polarization basis, formed by the disk normal
and the observer vector:

$\mbox{\boldmath$\hat{e}$}_{1}=\frac{\mbox{\boldmath$\hat{n}$}-\cos{i}\
\mbox{\boldmath$\hat{o}$}}{\sin{i}},\qquad\mbox{\boldmath$\hat{e}$}_{2}=\frac{\mbox{\boldmath$\hat{o}$}\times\mbox{\boldmath$\hat{n}$}}{\sin{i}}.$
(32)

The total PA in this basis is given by the sum of the relativistic rotations
and the intrinsic PA $\chi_{0}$,

$\chi=\chi_{0}+\chi^{\text{tot}}=\chi_{0}+\chi^{\text{SR}}+\chi^{\text{GR}}.$
(33)

The PA $\chi_{0}$ in the comoving frame is described in the polarization basis
that is formed by the local normal vector $\hat{n}$ and the photon momentum
vector in that frame $\mbox{\boldmath$\hat{k}$}^{\prime}_{0}$:

$\mbox{\boldmath$\hat{e}$}_{1}^{\prime}=\frac{\mbox{\boldmath$\hat{n}$}-\cos{\zeta^{\prime}}\
\mbox{\boldmath$\hat{k}$}_{0}^{\prime}}{\sin{\zeta^{\prime}}},\qquad\mbox{\boldmath$\hat{e}$}_{2}^{\prime}=\frac{\mbox{\boldmath$\hat{k}$}_{0}^{\prime}\times\mbox{\boldmath$\hat{n}$}}{\sin{\zeta^{\prime}}}.$
(34)

In our case of electron scattering atmosphere, the intrinsic PA
$\chi_{0}=\pi/2$.

The area of an element of the disk surface is expressed as

${\rm d}S=\frac{{R}_{\rm S}^{2}r\,{\rm d}r\,{\rm d}\varphi}{\sqrt{{\cal D}}}.$
(35)

It occupies the solid angle ${\rm d}\Omega$ on the sky and can be computed in
Schwarzschild approximation as

${\rm d}\Omega=\frac{{\rm d}S\cos\zeta}{D^{2}}{\cal L}=\frac{{R}_{\rm
S}^{2}}{D^{2}}\frac{r\,{\rm d}r\,{\rm d}\varphi}{\sqrt{\cal D}}{\cal
L}\cos\zeta,$ (36)

where $D$ is the distance to the source and the lensing factor ${\cal L}$ is
defined for Schwarzschild space-time as (Beloborodov 2002)

${\cal L}=\frac{1}{1-u}\frac{{\rm d}\cos\alpha}{{\rm d}\cos\psi}.$ (37)

It can be computed analytically (Poutanen 2020a)

${\cal L}=1+\frac{3u^{2}y^{2}}{112}-\frac{{\rm
e}}{100}uy\left[2\,\ln\left(1-\frac{y}{2}\right)+y\frac{1-3y/4}{1-y/2}\right]$
(38)

following Eq. (28).

Next, the intensity in the comoving frame ${I}^{\prime}_{E^{\prime}}$ is
related to the observed one as

$\bm{I}_{E}=g^{3}\bm{I}^{\prime}_{E^{\prime}},$ (39)

where $g$ is the redshift factor, which can be expressed in the case of a disk
in the equatorial plane of the BH as (see Eq. C13 in Li et al. 2005)

$g=E/E^{\prime}=\gamma\left[{\cal X}+{\cal Y}\beta+({\cal X}\beta+{\cal
Y})\cos\xi^{\prime}\right],$ (40)

where

$\displaystyle{\cal X}$ $\displaystyle=\sqrt{\cal D/\cal A\phantom{.}},$ (41)
$\displaystyle{\cal Y}$ $\displaystyle=a/\sqrt{4r^{4}\cal A\phantom{.}},$ (42)
$\displaystyle\cal A$ $\displaystyle=1+(r+1)\frac{a^{2}}{4r^{3}}.$ (43)

In the case of $a=0$, the redshift factor is reduced to
$g=\delta\sqrt{1-1/r}$.

The total disk flux is obtained by integration over the surface:

$\displaystyle\bm{F}_{E}$ $\displaystyle\equiv$
$\displaystyle\begin{bmatrix}F_{I}(E)\\\ F_{Q}(E)\\\
F_{U}(E)\end{bmatrix}=\int{\rm d}\Omega\begin{bmatrix}I_{E}\\\ Q_{E}\\\
U_{E}\end{bmatrix}$ (44) $\displaystyle=$ $\displaystyle\frac{{R}_{\rm
S}^{2}}{D^{2}}\int\limits_{{r}_{\rm{ISCO}}}^{r_{\textrm{out}}}\frac{r\,{\rm
d}r}{\sqrt{\cal D}}\int_{0}^{2\pi}{\rm d}\varphi\,g^{3}{\cal L}\cos\zeta\ {\bf
M}(r,\varphi)\bm{I}^{\prime}_{E^{\prime}}(\zeta^{\prime}),$

where ${\bf M}$ is the rotation matrix defined at every surface element of the
disk, ${\rm d}r\,{\rm d}\varphi$, as following

${\bf M}(r,\varphi)=\begin{bmatrix}1&0&0\\\ 0&\cos 2\chi^{\rm tot}&-\sin
2\chi^{\rm tot}\\\ 0&\sin 2\chi^{\rm tot}&\cos 2\chi^{\rm tot}\end{bmatrix}.$
(45)

The PD and PA are then defined from the Stokes parameters of the total flux.
The PD is obtained as

$p(E)=\frac{\sqrt{F^{2}_{Q}(E)+F^{2}_{U}(E)}}{F_{I}(E)},$ (46)

while the observed PA can be computed from either

$\tan{2\chi(E)}=\frac{F_{U}(E)}{F_{Q}(E)}$ (47)

or

$\tan{\chi(E)}=\frac{p(E)F_{I}(E)-F_{Q}(E)}{F_{U}(E)}.$ (48)

Following Loktev et al. (2022), we can use dimensionless energy $x=E/kT_{*}$
(and $x^{\prime}=x/g$) and scale the luminosity to $\sigma_{\rm
SB}T_{*}^{4}{R}_{\rm S}^{2}$ to get the dimensionless luminosity $l_{x}$ in
the following form

$\displaystyle x\bm{l}_{x}$ $\displaystyle=$
$\displaystyle\frac{60}{\pi^{4}}\int\limits_{{r}_{\rm{ISCO}}}^{r_{\textrm{out}}}\\!\\!\frac{r\,{\rm
d}r}{\sqrt{{\cal D}}}\int_{0}^{2\pi}{\rm d}\varphi\,\ {\cal L}\ \cos\zeta\
\frac{g^{4}}{f_{\rm c}^{4}}$ (49) $\displaystyle\times$
$\displaystyle\frac{x^{\prime 4}a_{\rm es}(\zeta^{\prime})}{{\rm
e}^{\displaystyle x^{\prime}t(r,a)/f_{\rm c}}-1}\begin{bmatrix}1\\\ p_{\rm
es}(\zeta^{\prime})\cos(2\chi)\\\ p_{\rm
es}(\zeta^{\prime})\sin(2\chi)\end{bmatrix},$

where $t(r,a)=r^{3/4}f^{-1/4}(r,a)$ is related to the disk temperature (Eq.
5). In this notation, the spectral shape and normalization of $xl_{x}$ are
independent of the BH mass and accretion rate.

Figure 2: Influence of relativistic effects on the PA and the emission angle.
The rotation of the PA (a sum of the GR and SR effects) $\chi^{\rm tot}$
(upper panels) and the emission angle $\zeta^{\prime}$ (lower panels) computed
at $r={r}_{\rm{ISCO}}$ (given in the upper left corner of each panel) for
spins $a=0.2$ (left panels), 0.5 (middle panels), and 0.8 (right panels) and
inclinations $i=30\degr$ (red), $60\degr$ (green), and $80\degr$ (blue lines).
The results using artpol code are shown with the solid lines, while the
results from the ray-tracing in Kerr metric are presented with the dashed
lines. The difference between the two values is shown under each panel.

Image of the disk can be reproduced using Cartesian coordinates $X$ and $Y$
(expressed in units of ${R}_{\rm S}$) on the plane of the sky:

$\begin{bmatrix}X\\\ Y\end{bmatrix}=b\begin{bmatrix}-\sin\Phi\\\
\cos\Phi\end{bmatrix}=\frac{r}{\sqrt{\cal
D}}\frac{\sin\alpha}{\sin\psi}\begin{bmatrix}\sin\varphi\\\ -\cos
i\cos\varphi\end{bmatrix},$ (50)

where $\Phi$ is the position angle of the point where the photon hits the
plane of the sky measured counterclockwise from the projection of the disk
axis on the sky and

$b=\frac{r}{\sqrt{\cal D}}\sin\alpha$ (51)

is the approximation for the impact parameter, which is exact in the
Schwarzschild case (Pechenick et al. 1983; Beloborodov 2002).

###  2.3 Numerical ray-tracing in Kerr metric

Let us now compare our results to numerical ray-tracing calculations to
determine the accuracy of artpol. We construct the image of the BH accretion
disk in the Kerr metric using the arcmancer code (Pihajoki et al. 2018),
designed to integrate the exact geodesic equation numerically. The resulting
trajectories are parametrized by the point at which they intersect the
(observer) image plane.

Figure 3: Same as Fig. 2 for the spin of $a=0.8$ and three ring radii of $r=5$
(left panels), 3 (middle panels), and 2 (right panels).

To obtain an image of the accretion disk, we first define the region of
interest on the image plane to include the relevant region of the disk. The
image plane is positioned at a distance $r_{\rm img}=2500$, where the effects
of general relativity on the photon trajectory are negligible. We then
introduce a grid of nodes and calculate the observed radiation intensity from
the accretion disk at each node. The rectangular grid on the image plane is
orthogonal to the vector pointing toward the observer. At each point on the
grid, we define a local Lorentz frame with the basis vectors along Boyer-
Lindquist coordinate grid. The Lorentz frames constitute the initial
conditions for the geodesic curves with the tangent vectors normal to the
image plane. Essentially, a set of parallel rays are propagated from the
points of the grid towards the BH, backward in time, until they either
intersect the disk surface or reach a distance $>\\!r_{\mathrm{img}}$ (as
measured from the BH). When all the geodesics are calculated, we compute the
local intensity and polarization for each point that intersects the disk
surface, which we assume lies in the equatorial plane of the Kerr BH. We
consider the velocity and energy dissipation profile following Novikov &
Thorne (1973).

The Stokes parameters are parallel transported in the “polarization frame”
along the geodesic in Kerr space-time. The polarization frame
$\mathcal{P}=\\{{v}^{a},{h}^{a}\\}$ is a pair of orthogonal space-like vectors
that are orthogonal to the geodesic itself and the four-velocity of the
observer. Initially, at the image plane, the polarization frame consists of
the corresponding two vectors of the local Lorentz frame. Then, the
polarization frame is parallel-transferred along the geodesic to the disk
surface to obtain the polarization parameters. Then, the frame is projected to
the rest frame of the fluid in the disk (comoving frame). The propagation of
the geodesics, transporting of the Lorentz frames, and the projection of the
screen to the comoving frame are automatically performed by arcmancer.

All the necessary values, namely, the redshift $g$, the rotation of the
polarization basis $\chi^{\rm tot}$, and the emission angle $\zeta^{\prime}$
can be defined in the comoving frame, where the geodesic hits the disk
surface. The gravitational redshift factor $g$ is computed using the curve
tangents of the initial and the endpoints of the light trajectory. For more
details on this numerical procedure, we refer to Sect. 4 of Pihajoki et al.
(2018). We also correct the redshift by the factor
$g_{0}=\sqrt{1-1/r_{\text{img}}}$ to mitigate the redshift between the
observer at $r_{\text{img}}$ and the one at infinity. The total redshift
agrees with the one given by the analytical expression given by Eq. (40).

To compute the overall flux from the disk in the case of numerical ray-
tracing, we sum the intensities over the whole image grid as

$\bm{F}_{E}=\frac{{R}_{\rm S}^{2}}{D^{2}}\iint\limits_{r<r_{\textrm{out}}}{\rm
d}x\ {\rm d}y\ {\bf M}(x,y)\ \bm{I}^{\prime}_{E^{\prime}}(r,\zeta^{\prime}),$
(52)

where ${\rm d}x\,{\rm d}y$ is the size of the pixel at $r_{\text{img}}$,
measured in Schwarzschild radii, ${\bf M}(x,y)$ is the rotation (Mueller)
matrix along the trajectory from the point $(x,y)$ on the image plane to the
point on the disk surface with a coordinate $r<r_{\textrm{out}}$.

The secondary, and higher-order images of the disk and its bottom side are
visible through the gap between the ${r}_{\rm{ISCO}}$ and the BH event
horizon. From there, a segment of the disk surface can be seen multiple times,
depending on how many revolutions its photons make around the BH before
escaping toward the observer. These trajectories are not accounted for in
artpol because we find the contribution of the secondary images to the
observed flux are below 0.5% (see Sect. 3.3 below) for any inclination and
spin values (see also Zhou et al. 2020). Lastly, we note that it is also
possible to compute polarization rotation using Walker–Penrose theorem (e.g.,
see the appendix in Li et al. 2009); however, it does not significantly
decrease the computing time as long as the geodesics are computed numerically
anyway.

Figure 4: Polarization parameters of radiation from the disk at inclinations
$i=30\degr$ (upper row), $60\degr$ (middle row), and $80\degr$ (lower row),
and different spin parameters $a=0.2$ (left column), 0.5 (middle column), and
0.8 (right column). The black lines outline an even polar grid on the disk
with rays spaced by $30\degr$ in azimuth and contours spaced by $5{R}_{\rm S}$
in radius (the corresponding values are denoted in each ring). The sticks
represent the polarization parameters derived for the center of each grid
segment. The polarized disk emission were calculated using the exact numerical
interaction of the geodesics in the Kerr metric (with the arcmancer code;
green lines) and using the approximate analytical formulae in the
Schwarzschild metric (with the artpol code; cyan lines). The coordinates of
the cyan sticks were computed using Eq. (50). The length of the sticks is
proportional to the observed PD from the region, and orientation shows the
observed PA. The 5% polarization stick is shown in the bottom right panel to
scale. The colors reflect the logarithm of the bolometric intensity relative
to the maximal value across all panels. Note that differences between the
exact numerical solution and the approximate artpol method are visible only at
the inner radii for systems with the highest spin values.

##  3 Applications

In this section, we show the applications of the artpol technique to spectra
and polarization signatures of matter near a Kerr BH. We first verify the
accuracy of the artpol for the PA rotation $\chi$ and the emission zenith
angle $\zeta^{\prime}$ (measured in the comoving frame) for the case of the
narrow disk ring located at different distances from the BH. We proceed to the
comparison of polarized images of the accretion disk. Finally, we show the
spectro-polarimetric energy distributions of the accretion disk obtained with
analytical and numerical approaches. In all cases, we assume that the matter
(ring/disk) rotates counterclockwise. Positive spin values correspond to the
prograde rotation (BH spin vector aligned with the orbital momentum of the
accretion disk).

Figure 5: Spectra of the luminosity $xl_{x}$ (upper panels), PD (middle
panels) and PA (bottom panels) for different inclinations $i=30\degr$ (left),
$60\degr$ (middle), and $80\degr$ (right column) and BH spin $a=-1$ (pink
lines), $0$ (red), $0.5$ (orange), $0.8$ (green), and $0.94$ (blue). The solid
lines correspond to the results obtained with the approximate analytical
formulae (artpol), while the dashed lines correspond to the exact numerical
integration of the geodesics (arcmancer). The accretion disk extends from
${r}_{\rm{ISCO}}$ to $r_{\rm{out}}=3000$ for each case. Narrow panels show the
difference between the results given by the two methods. The black dotted
lines denote the results for a Newtonian disk from Loktev et al. (2022).

###  3.1 Narrow ring

The difference between Stokes parameters computed via artpol and numerical ray
tracing is expected to be high at the ${r}_{\rm{ISCO}}$, as the frame-dragging
effects – omitted in artpol – are most important here. In Fig. 2, we show the
total rotation angle $\chi^{\rm tot}$ and the emission zenith angle
$\zeta^{\prime}$ at the ISCO computed by two methods for different spin values
and observer inclinations. For artpol, both PA rotation and $\zeta^{\prime}$
are systematically shifted, over the azimuth, with respect to the numerical
values. The difference is highest at $\varphi\sim 180\degr$ and, in general,
is larger for higher inclinations. The effect is caused by the close BH
approach of the photon trajectories starting from regions of the disk located
behind the BH. Naturally, the difference is larger for higher spin values. For
example, in the case of $a=0.8$ (${r}_{\rm{ISCO}}\approx 1.45$) and
$i=30\degr$ the maximum $|\Delta\chi^{\rm tot}|$ is just $10\degr$ and
$|\Delta\zeta^{\prime}|$ is only $9\degr$. We also note here that the largest
error on $\chi^{\rm tot}$ corresponds to a small emission angle
$\cos\zeta^{\prime}\sim 1$ at which PD is minimal (see Eq. 14); therefore, the
polarized flux is affected very little.

While the accuracy of artpol can be low for certain parts of the narrow ring
of matter at the ISCO, for the case of the zero torque boundary condition
(Novikov & Thorne 1973), the total flux emerging from this ring is zero. The
error is then weighted with the dissipation profile, which achieves maximum at
distances further than ISCO. In Fig. 3 we show $\chi^{\rm tot}$ and
$\zeta^{\prime}$ for $a=0.8$ at different radii $r=2$, 3, and 5, computed with
artpol and numerical ray-tracing. The shapes of the lines are similar, but the
deviations between the two methods are much smaller. For example, the maximum
$|\Delta\chi^{\rm tot}|$ and $|\Delta\zeta^{\prime}|$ for $r=2$ and
$i=30\degr$ are $5\degr$ and $4\degr$, respectively, and those are peak values
only at $\varphi\approx 180\degr$. We note that for $r>5$, artpol gives
$\chi^{\rm tot}$ and $\zeta^{\prime}$ that are nearly indistinguishable from
the numerical results, even for BHs with extreme spins.

###  3.2 Disk imaging

artpol can be used to construct images of the disk (see more details in Loktev
et al. 2022). In Fig. 4, we show the observed intensity and polarization of an
accretion disk for BH spin parameters $a=0.2$, 0.5, and 0.8, and observer
inclinations of $i=30\degr$, 60°, and 80°. The observed intensity for the
innermost parts of the standard disk with the electron scattering dominating
atmosphere is color-mapped. The sticks depict the polarization vector computed
with the ray-tracing and analytical techniques in the Schwarzschild metric.
The difference in polarization parameters can barely be seen only in the
innermost regions for high spin cases. Otherwise, the effect of BH rotation is
so tiny that the sticks visually coincide. The sticks are displaced
predominantly in the azimuthal direction, mainly caused by the frame-dragging
effect, while the radial displacement due to the approximate impact parameter
estimation is negligible. The most distorted region of the disk appears to be
directly behind the black hole, where the surface of the disk appears to be
dragged with the BH spin. This effect is, again, especially pronounced for
higher observer inclinations.

###  3.3 Disk spectra and polarization

Using the artpol technique, we perform fast computations of the spectral
properties, which are essential to extract the spin information in both
reflection spectroscopy and continuum-fitting methods and spectral dependence
of the polarization signatures. Comparison of the computed spectral energy
distribution, PD, and PA for the case of the optically thick, razor-thin
accretion disk to the exact ray-tracing calculations are shown in Fig. 5. The
local model for the disk is described in Sect. 2.1.

In the upper panels of Fig. 5, we show the multicolor blackbody accretion disk
spectra viewed at inclinations $i=30\degr$ (left column), $60\degr$ (middle
column) and $80\degr$ (right column). We consider different BH spins from a
retrograde one $a=-1$ to the prograde $a=0.94$. The solid lines correspond to
the calculations using artpol, and the dashed lines correspond to numerical
ray-tracing calculations. The differences between the methods are shown below
each panel. The deviations increase with the spin value, energy, and
inclination, mainly due to the increased importance of the frame-dragging
effects. These are more important for the innermost radii, contributing to
higher energies and the photon trajectories passing close to the BH, which
appear mostly at higher inclinations. For higher spin values, the ISCO radii
are smaller, and the frame-dragging effect is more pronounced. At the same
time, in the case of the retrograde BH rotation (BH rotates clockwise in our
case), $a=-1$, smaller deviations are caused by the more distant location of
the ISCO (${r}_{\rm{ISCO}}=4.5$). We find that for all the considered spin
values, the deviations in the flux values are smaller than 20% at energies
where the flux is relatively high (i.e., where it drops by less than a factor
of 10 from the peak value).

We can ask here a question whether the difference between the fluxes computed
with the two methods has its origin in the presence of secondary images. Using
arcmancer we computed the accretion disk spectra for two extreme BH spins and
two inclinations. The results shown in Fig. 6 demonstrate that the secondary
images contribute at most 0.5% of the flux for the high spin BH and large
inclination. At lower inclinations or smaller spin values, their contribution
is even smaller.

The middle and lower panels of Fig. 5 show the observed PD and PA. The trend
of increasing deviations towards higher spins, energies, and inclinations is
also evident in these quantities. For all considered cases, however, the
deviations remain small as compared to the typical observational errors
achieved in IXPE observations. The maximum deviations in the predicted PD are
about 0.4% for $a=0.94$ and are much smaller than the expected accuracy of
measurements with IXPE for smaller spins. The error in PA is below 1° for
$a\leq 0.5$, and even for $a=0.94$, it is at most 5° at energies where the
flux is relatively high. This indicates that fitting the Stokes spectra
performed with our analytical technique will work fairly well for all
inclinations and all considered spins.

Figure 6: Contribution of secondary images to the observed accretion disk flux
for two inclinations $i=30\degr$ (left panel) and $80\degr$ (right panel) and
three BH spins $a=-1$ (pink line), $0$ (red), and $0.95$ (blue) as computed
with the numerical ray-tracing code arcmancer. The contribution of the
secondary images for $i=30\degr$ and $a=0.95$ is below $10^{-5}$. The solid
lines correspond to the full luminosity and the dashed lines to the
contributions of the secondary images.

We note the difference between the values obtained for the Newtonian disk
(Shakura 1973) in the Schwarzschild metric in Loktev et al. (2022) and the
zero-spin case considered here. The difference comes solely from the
relativistic temperature profile in Eq. (5), which is different from the
Newtonian temperature profile from Eq. (7) assumed in Loktev et al. (2022).
Notably, the Newtonian disk polarization profiles resemble relativistic
profiles with the spin parameter $a\sim 0.2$ more than $a=0$.

The reduction in calculation time using artpol compared to arcmancer is
consistent with that reported in Loktev et al. (2022), as we are comparing the
same two methods. Specifically, our calculations of one image using the
approximate analytical formulae are performed in $0.01$ s, while the
calculations of one polarized spectrum take about $0.05$ s. This has to be
compared to calculating an image $400\times 400$ pixels using a ray-tracing
code arcmancer, which takes about $1000$ s. Moreover, our analytical method
offers flexibility, allowing for its application to arbitrary disk geometries,
different energy dissipation profiles, and various local spectra and
polarization.

##  4 Summary

We presented the analytical technique artpol for computing the images and
spectro-polarimetric characteristics of relativistic accretion disks around
Kerr BHs. While the velocity and energy dissipation profiles used in the
calculations correspond to the space-time around spinning BHs, the photon
paths have been considered in the Schwarzschild metric. We showed that artpol
technique is highly efficient, allowing for reduction of the computing time by
a factor of $10^{4}$, and remains accurate, in general within 10% in flux,
0.2% in PD, and $1\degr$ in PA, for spin values $a\leq 0.5$. This enables
artpol technique to be used for a broad range of BH spin parameters. The
deviations from the exact results obtained by numerical ray-tracing techniques
systematically increase towards the highest spins and inclinations, reducing
the accuracy down to 20% in flux, 0.5% in PD, and $7\degr$ in PA for
$a=0.8$–0.94 and $i=80\degr$. The systematic discrepancy arises from the
frame-dragging effect, which is omitted in our analytical ray-tracing
calculations.

Applications of the analytical technique include fast computations of static
spectro-polarimetric signatures of accretion disks of X-ray binaries in the
soft/intermediate/very high states (Ratheesh et al. 2023; Rodriguez Cavero et
al. 2023; Rawat et al. 2023a, b; Jana & Chang 2023; Podgorny et al. 2023), and
account for the relativistic effects on spectra and polarimetric signatures of
Comptonization observed in the hard state of X-ray binaries and Seyfert
galaxies (Krawczynski et al. 2022; Marinucci et al. 2022; Gianolli et al.
2023; Ingram et al. 2023; Tagliacozzo et al. 2023). The long-anticipated
detection of the quasi-periodic oscillations from BH X-ray binaries in the
polarized light would open new prospects for artpol in terms of fast
calculation of the phase-resolved characteristics. The method is also well
suited for the warped disks whose spectro-polarimetric signatures have not
been studied in detail so far.

## Acknowledgments

We acknowledge support from the Jenny and Antti Wihuri foundation (VL), the
Academy of Finland grants 355672 (AV), 333112 (JP), and 349144 (VL, AV, JP),
and the German Academic Exchange Service (DAAD) project 57525212 (VFS). VFS
thanks Deutsche Forschungsgemeinschaft (DFG) grant WE 1312/59-1 for financial
support.

## References

  * Beloborodov (2002) Beloborodov, A. M. 2002, ApJ, 566, L85
  * Cárdenas-Avendaño et al. (2023) Cárdenas-Avendaño, A., Lupsasca, A., & Zhu, H. 2023, Phys. Rev. D, 107, 043030
  * Chan et al. (2013) Chan, C.-k., Psaltis, D., & Özel, F. 2013, ApJ, 777, 13
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Radiative transfer (New York: Dover)
  * Chandrasekhar & Breen (1947) Chandrasekhar, S., & Breen, F. H. 1947, ApJ, 105, 435
  * Connors et al. (1980) Connors, P. A., Piran, T., & Stark, R. F. 1980, ApJ, 235, 224
  * Connors & Stark (1977) Connors, P. A., & Stark, R. F. 1977, Nature, 269, 128
  * Dexter (2016) Dexter, J. 2016, MNRAS, 462, 115
  * Dexter & Agol (2009) Dexter, J., & Agol, E. 2009, ApJ, 696, 1616
  * Dovčiak et al. (2008) Dovčiak, M., Muleri, F., Goosmann, R. W., Karas, V., & Matt, G. 2008, MNRAS, 391, 32
  * Draghis et al. (2023) Draghis, P. A., Miller, J. M., Zoghbi, A., et al. 2023, ApJ, 946, 19
  * Gianolli et al. (2023) Gianolli, V. E., Kim, D. E., Bianchi, S., et al. 2023, MNRAS, 523, 4468
  * Ingram et al. (2015) Ingram, A., Maccarone, T. J., Poutanen, J., & Krawczynski, H. 2015, ApJ, 807, 53
  * Ingram et al. (2023) Ingram, A., Ewing, M., Marinucci, A., et al. 2023, MNRAS, submitted, arXiv:2305.13028
  * Jana & Chang (2023) Jana, A., & Chang, H.-K. 2023, MNRAS, submitted, arXiv:2307.14604
  * Kato et al. (2008) Kato, S., Fukue, J., & Mineshige, S. 2008, Black-Hole Accretion Disks — Towards a New Paradigm (Kyoto: Kyoto Univ. Press)
  * Kerr (1963) Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237
  * Krawczynski (2012) Krawczynski, H. 2012, ApJ, 754, 133
  * Krawczynski & Beheshtipour (2022) Krawczynski, H., & Beheshtipour, B. 2022, ApJ, 934, 4
  * Krawczynski et al. (2022) Krawczynski, H., Muleri, F., Dovčiak, M., et al. 2022, Science, 378, 650
  * Li et al. (2009) Li, L.-X., Narayan, R., & McClintock, J. E. 2009, ApJ, 691, 847
  * Li et al. (2005) Li, L.-X., Zimmerman, E. R., Narayan, R., & McClintock, J. E. 2005, ApJS, 157, 335
  * Lightman & Shapiro (1975) Lightman, A. P., & Shapiro, S. L. 1975, ApJ, 198, L73
  * Loktev et al. (2020) Loktev, V., Salmi, T., Nättilä, J., & Poutanen, J. 2020, A&A, 643, A84
  * Loktev et al. (2022) Loktev, V., Veledina, A., & Poutanen, J. 2022, A&A, 660, A25
  * Loskutov & Sobolev (1979) Loskutov, V. M., & Sobolev, V. V. 1979, Astrofizika, 15, 241
  * Loskutov & Sobolev (1981) —. 1981, Astrofizika, 17, 97
  * Marinucci et al. (2022) Marinucci, A., Muleri, F., Dovciak, M., et al. 2022, MNRAS, 516, 5907
  * McClintock et al. (2014) McClintock, J. E., Narayan, R., & Steiner, J. F. 2014, Space Sci. Rev., 183, 295
  * Miller et al. (2009) Miller, J. M., Reynolds, C. S., Fabian, A. C., Miniutti, G., & Gallo, L. C. 2009, The Astrophysical Journal, 697, 900
  * Motta et al. (2022) Motta, S. E., Belloni, T., Stella, L., et al. 2022, MNRAS, 517, 1469
  * Motta et al. (2014) Motta, S. E., Belloni, T. M., Stella, L., Muñoz-Darias, T., & Fender, R. 2014, MNRAS, 437, 2554
  * Narayan et al. (2021) Narayan, R., Palumbo, D. C. M., Johnson, M. D., et al. 2021, ApJ, 912, 35
  * Nättilä & Pihajoki (2018) Nättilä, J., & Pihajoki, P. 2018, A&A, 615, A50
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus) (New York: Gordon & Breach), 343
  * Page & Thorne (1974) Page, D. N., & Thorne, K. S. 1974, ApJ, 191, 499
  * Pechenick et al. (1983) Pechenick, K. R., Ftaclas, C., & Cohen, J. M. 1983, ApJ, 274, 846
  * Pihajoki et al. (2018) Pihajoki, P., Mannerkoski, M., Nättilä, J., & Johansson, P. H. 2018, ApJ, 863, 8
  * Pineault (1977) Pineault, S. 1977, MNRAS, 179, 691
  * Pineault & Roeder (1977a) Pineault, S., & Roeder, R. C. 1977a, ApJ, 212, 541
  * Pineault & Roeder (1977b) —. 1977b, ApJ, 213, 548
  * Podgorny et al. (2023) Podgorny, J., Marra, L., Muleri, F., et al. 2023, MNRAS, submitted, arXiv:2303.12034
  * Poutanen (2020a) Poutanen, J. 2020a, A&A, 640, A24
  * Poutanen (2020b) —. 2020b, A&A, 641, A166
  * Ratheesh et al. (2023) Ratheesh, A., Dovčiak, M., Krawczynski, H., et al. 2023, ApJ, submitted, arXiv:2304.12752
  * Rawat et al. (2023a) Rawat, D., Garg, A., & Méndez, M. 2023a, ApJ, 949, L43
  * Rawat et al. (2023b) —. 2023b, MNRAS, 525, 661
  * Rees (1975) Rees, M. J. 1975, MNRAS, 171, 457
  * Rodriguez Cavero et al. (2023) Rodriguez Cavero, N., Marra, L., Krawczynski, H., et al. 2023, ApJ, submitted, arXiv:2305.10630
  * Shakura (1973) Shakura, N. I. 1973, Soviet Astronomy, 16, 756
  * Shimura & Takahara (1995) Shimura, T., & Takahara, F. 1995, ApJ, 445, 780
  * Sobolev (1949) Sobolev, V. V. 1949, Uch. Zap. Leningrad Univ., 16
  * Sobolev (1963) —. 1963, A treatise on radiative transfer (Princeton: Van Nostrand)
  * Stark & Connors (1977) Stark, R. F., & Connors, P. A. 1977, Nature, 266, 429
  * Tagliacozzo et al. (2023) Tagliacozzo, D., Marinucci, A., Ursini, F., et al. 2023, MNRAS, submitted, arXiv:2305.10213
  * Uttley et al. (2014) Uttley, P., Cackett, E. M., Fabian, A. C., Kara, E., & Wilkins, D. R. 2014, A&A Rev., 22, 72
  * Viironen & Poutanen (2004) Viironen, K., & Poutanen, J. 2004, A&A, 426, 985
  * Walker & Penrose (1970) Walker, M., & Penrose, R. 1970, Communications in Mathematical Physics, 18, 265
  * Weisskopf et al. (2022) Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, JATIS, 8, 026002
  * Yang & Wang (2014) Yang, X.-L., & Wang, J.-C. 2014, A&A, 561, A127
  * Zhang et al. (2019) Zhang, W., Dovčiak, M., & Bursa, M. 2019, ApJ, 875, 148
  * Zhou et al. (2020) Zhou, M., Ayzenberg, D., Bambi, C., & Nampalliwar, S. 2020, Phys. Rev. D, 101, 043010