#  CHIANTI – an atomic database for emission lines - Paper XVII: Version 10.1,
revised ionization and recombination rates and other updates

Kenneth P. Dere  Department of Physics and Astronomy, George Mason University,
4400 University Drive, Fairfax, VA 22030, USA  P. R. Young NASA Goddard Space
Flight Center, Code 671, Greenbelt, MD 20771, USA  Northumbria University,
Newcastle Upon Tyne NE1 8ST, UK  G. Del Zanna DAMTP, Center for Mathematical
Sciences, University of Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK
E. Landi Department of Climate, Space Sciences and Engineering, University of
Michigan, Ann Arbor, MI, 48109

###### Abstract

The CHIANTI atomic database provides sets of assessed data used for simulating
spectral observations of astrophysical plasmas. This article describes updates
that will be released as version 10.1 of the database. A key component of
CHIANTI is the provision of ionization and recombination rates that are used
to compute the ionization balance of a plasma over a range of temperatures.
Parameters for calculating the ionization rates of all stages of ions from H
through Zn were compiled and inserted into the CHIANTI database in 2009. These
were based on all measurements that were available at the time and
supplemented with distorted wave calculations. Since then, there have been a
number of new laboratory measurements for ions that produce spectral lines
that are commonly observed. Parameters have been fit to these new measurements
to provide improved ability to reproduce the ionization cross sections and
rate coefficients, and these are added to the database. CHIANTI 10.1 also
includes new recombination rates for the phosphorus isoelectronic sequence,
and the updated ionization and recombination rates have been used to calculate
a new ionization equilibrium file. In addition, CHIANTI 10.1 has new electron
collision and radiative datasets for eight ions in the nitrogen and oxygen
isoelectronic sequences, and updated energy level and wavelength data for
seven other ions.

atomic data — atomic processes — Sun: UV radiation — Sun: X-rays, gamma rays —
Ultraviolet: general — X-rays: general

111Released on May 23, 2023

##  1 Introduction

CHIANTI is an atomic database and software package used for modeling
optically-thin emission from astrophysical plasmas. It has been an open-source
project since the first release in 1996 (Dere et al., 1997), and the data and
software are available at https://chiantidatabase.org. The most recent
previous release is CHIANTI 10 (Del Zanna et al., 2021), and the present
article describes the updates and new features of the 10.1 release. Young et
al. (2016) and Del Zanna & Young (2020) summarize the contents of the modern
version of the database and provide applications.

The spectral emission from ionized astrophysical plasmas is a key to our
understanding of the processes that form and govern these regions. A key to
the analysis of these spectra is a knowledge of the ionization state. For
plasmas in a collisional ionization equilibrium, such as the solar corona, the
equilibrium is largely controlled by electron collisional ionization and
electron collisional recombination. Collisional ionization consists of both
direct ionization (DI) and excitation followed by autoionization (excitation-
autoionization or EA). Recombination also consists of two different processes,
radiative recombination (RR) and dielectronic recombination (DR). Dere (2007)
analyzed a large number of measured ionization cross sections in order to make
ionization cross sections and their associated ionization rate coefficients
readily available through the CHIANTI atomic database (Dere et al., 1997,
2009; Del Zanna et al., 2021). Since that analysis, a number of important
ionization cross sections have been measured and published. Here, these new
measurements are analyzed in order to improve the calculation of ionization
cross sections and rate coefficients in CHIANTI.

Recombination has been the focus of a project described by Badnell et al.
(2003) to systematically generate new DR and RR data for entire isoelectronic
sequences. Previous versions of CHIANTI have used these rates for sequences
from hydrogen to silicon. New calculations for the phosphorus sequence were
presented by Bleda et al. (2022) and these are discussed in Sect. 2.15. A new
ionization balance is calculated from the updated ionization and recombination
rates and is discussed in Sect. 2.16.

The level balance within an ion is determined by electron collision strengths,
radiative decay rates and energy levels, and updates to these datasets are
described in Sect. 3.

##  2 Revised ionization and recombination rate coefficients

###  2.1 Approach

The approach used by Dere (2007) was to compile a set of parameters for each
ion to be able to reproduce the available ionization cross sections
measurements. The first priority was to use measured cross sections and when
these were not available, cross sections were calculated with the Flexible
Atomic Code (FAC: Gu, 2002). In this paper, new measurements for 13 ions are
examined. In all cases, this leads to new fits to laboratory measurements or
adjustments to existing theoretical cross sections.

Dere (2007) developed a scaling for ionization cross sections similar to that
of Burgess & Tully (1992) for collision strengths. We define a scaled energy
$U$,

$U=1-\frac{\ln f}{\ln(u-1+f)}$ (1)

and a scaled cross section $\Sigma$,

$\Sigma=\frac{u\,\sigma\,I^{2}}{\ln u+1}$ (2)

where $u=E/I$, $E$ is the energy of the incident electron, $I$ is the
ionization potential, $\sigma$ is the ionization cross-section, and $f$ is a
scaling parameter that is selected by the data assessor.

The scaled energy $U$ varies from 0 at $I$ to unity for infinite energy. The
scaling parameter adjusts the placement of the scaled energies between these
two limits. The scaling is chosen to spread the region just above the
ionization potential so that the variation of $\Sigma$ becomes more apparent
and the fit to the measured values easier to evaluate. We refer to this
process as BTI (Burgess-Tully-Ionization) scaling. A spline fit to the
$\Sigma$ values is performed after $f$ has been chosen, and the spline values
are stored in the diparams CHIANTI file. Typically nine or more spline nodes
are used. Since the low energy behavior of $\Sigma$ is very apparent, it is
possible to determine if additional spline nodes are needed. This is often the
case at the energy region just above the ionization potential, that is
important for determining the electron ionization rate coefficient. This is
illustrated in Sect. 2.14.

An example of fitting BTI scaled measured cross sections is shown in Fig. 1
for S xiii. Here the BTI scaled measurements of Hahn et al. (2012a) are shown
together with the current fit to the BTI scaled measurements below 2360 eV
($U$ = 0.55). The fit to the measurements provides the parameters to reproduce
the direct ionization cross sections of the ion in question. The high energy
limit, $U$ = 1, is given by the Bethe cross section as discussed by Dere
(2007). The deviation of the measurements from the fit between $U=0.55$ and
$U=0.60$ shown in Fig. 1 is due to the excitation-autoionization (EA)
component. Dere (2007) computed the EA parameters using FAC and in the present
work adjustments are made to the threshold and/or the amplitude to best match
the new measurements.

The procedure for each ion examined is discussed in more detail later in the
following Sections.

Finally, we note that Dufresne & Del Zanna (2019) and Dufresne et al. (2020)
recently calculated ionization cross sections for all ions of carbon and
oxygen using FAC and AUTOSTRUCTURE (Badnell, 2011), and provided comparisons
to the available experimental data. The calculations were focused on
ionization from metastable levels, and will be made available in a future
CHIANTI release. Cross sections from the present work are retained for carbon
and oxygen to produce a consistent set of ionization data for all elements.
Cross sections for collisional ionization from ground states in those works
were in most cases within 15% of Dere (2007).

Figure 1: BTI scaled ionization cross sections of S xiii. The label Hahn
refers to Hahn et al. (2012a), and present to the results reported here.

###  2.2 The beryllium isoelectronic sequence

####  2.2.1 O V

Crossed beams measurements of the ionization cross sections of C iii, N iv,
and O v made at the Oak Ridge National Laboratory have been presented by Fogle
et al. (2008). The measurements of C iii and N iv show considerable cross
sections below the ionization potential, indicating a beam containing a
significant population in the 2s2p 3P0 metastable level. These measurements
will not be used in this analysis but the O v measurements do not indicate the
presence of significant populations of excited levels.

As in $\S$ 2.2.2, the direct ionization cross sections are represented by a
fit to the BTI scaled cross sections of Fogle et al. (2008) below about 520 eV
where EA processes become significant. The EA cross sections of Dere (2007)
have been increased by a factor of 2.5 to match the measurements but the two
EA thresholds remain at 552 and 620 eV. The O v measurements of Fogle et al.
(2008), Dere (2007) and the current analysis are presented in Fig. 2. The peak
cross sections near about 310 eV show that the current values are reduced by a
factor of about 0.84 from the FAC values of Dere (2007).

We note that earlier experimental data had significant limitations. The
Crandall et al. (1979) and Falk et al. (1983) experimental results were in
agreement with each other, but they were not able to determine the metastable
fraction present. Loch et al. (2003) produced similar experimental values,
although significant readings were present below threshold, indicating a
linear offset in readings or contamination in the beam. All of the above cross
sections were higher than theory by 30-50%. The Fogle et al. (2008)
experimental cross sections were in good agreement with their R-Matrix
calculations at low and high energies, while at the peak they were about 25%
below. The Dufresne et al. (2020) calculations are in good agreement with the
Fogle et al. (2008) ones.

Figure 2: Ionization cross sections of O v. The label Fogle refers to Fogle et
al. (2008), 2007 refers to Dere (2007), and present to the results reported
here.

####  2.2.2 S XIII

Storage ring measurements of S xiii ionization cross sections have been
presented by Hahn et al. (2012a). A major problem with experimental
measurements of ions in the beryllium isoelectronic sequence is that the ion
beams can contain a significant population of ions in the 2s2p 3P0 metastable
level. This results in cross section measurements that show an enhanced cross
section below the ionization potential of the ground level. As the authors
explain, they used a beam of the 33S isotope that has a nuclear spin. “The
resulting hyperfine interaction inducing a mixing of the 3P1 and 3P0 levels,
decreasing the lifetime of the 3P0 level.” The beam of S xiii ions were kept
in the ion storage ring long enough that the metastable levels decay to the
ground state.

Fig. 3 shows the present and previous cross sections together with the values
of Hahn et al. (2012a). The original FAC calculations indicated excitation-
autoionization contributions at energies of 2412 and 2786 eV for excitations
to the 1s 2s2 2p and 1s 2s2 3l levels, respectively. The energies of the EA
contributions are consistent with the measurements but the magnitudes of both
have been decreased by a factor of 0.75 to be consistent with the
measurements. The ratio of the maximum cross section at about 1600 eV is about
1.03 for the present results to the 2007 FAC values and the ratio at a lower
energy of 1000 eV is about 1.10.

Figure 3: Ionization cross sections of S xiii. The label Hahn refers to Hahn
et al. (2012a), 2007 refers to Dere (2007), and present to the results
reported here.

###  2.3 The boron isoelectronic sequence: Mg VIII

Hahn et al. (2010) reports measurements of the ionization cross sections of Mg
viii performed at the heavy ion storage ring TSR at the Max-Planck-Institut
für Kernphysik. They compare their measurements of the FAC cross sections to
Dere (2007) and find that, near the peak of the cross section, their
measurements are below the FAC calculations but within the error bars. At
energies just above threshold near 300 eV, their measurements are above the
FAC values and at the very lowest energies the differences exceed their error
bars. While the earlier FAC calculations were not too far from the new
measurements, it is worthwhile to perform a fit to the Hahn et al. (2010)
measurements. The FAC calculations of Dere (2007) only considered direct
ionization from the 2s and 2p levels and no EA contributions.

The BTI scaled experimental cross sections have been fit in the manner
described in $\S$ 2.2.2 with particular attention paid to the measurements
just above threshold. The Mg viii measurements of Hahn et al. (2010), the FAC
cross sections of Dere (2007) and the current fits to the experimental data
are presented in Fig. 4. As the authors suggest, “the apparent structure in
the experimental results near 600, 1400 and 1600 eV are due to run to run
variations in the background levels.” At peak cross sections near about 700 eV
the current values are reduced by a factor of about 0.92 from the FAC values
of Dere (2007). However, at an energy of about 300 eV the current values are
about a factor of 1.3 above the FAC cross sections.

Figure 4: Ionization cross sections of Mg viii. The label Hahn refers to Hahn
et al. (2010), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.4 The fluorine isoelectronic sequence: Fe XVIII

The ion Fe xviii provides a number of strong lines between 14 and 18 Å. These
make important contributions to observations by Chandra and other missions
that observe at X-ray wavelengths. Hahn et al. (2013) reported measurements of
the cross sections of Fe xiv, Fe xvii and Fe xviii at the heavy ion storage
ring TSR, and Fe xviii is discussed here.

The BTI scaled experimental cross sections of Hahn et al. (2013) have been fit
in the manner described in $\S$ 2.1 with particular attention paid to the
measurements just above threshold. A comparison of the measurements of Hahn et
al. (2013), Dere (2007) and the present results are shown in Fig. 5. Hahn et
al. (2013) suggest that an EA contribution to the cross section could arrive
from excitation of an inner shell 2s electron but there does not seem to be
any evidence of this in the measurements. Neither do the measurements of the
ionization cross sections of Si vi, also in the fluorine sequence, by Thompson
& Gregory (1994) indicate any significant EA contributions. The FAC
calculations of Dere (2007) only considered direct ionization from the 2s and
2p levels and no EA contributions.

At energies of about 2000 eV, somewhat above threshold, the present cross
section is about a factor of 1.13 greater than the 2007 values, but the ratio
becomes unity at 3200 eV and 0.94 at about 5000 eV.

Figure 5: Ionization cross sections of Fe xviii. The label Hahn refers to Hahn
et al. (2013), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.5 The neon isoelectronic sequence: Fe XVII

Fe xvii is a strong contributor to the X-ray spectrum in the wavelength range
10 to 18 Å and at extreme ultraviolet (EUV) wavelengths between 200 and 300 Å.
A number of the EUV lines have been observed in solar spectra between 180 and
400 Å (Warren et al., 2008; Del Zanna & Ishikawa, 2009).

Measurements of the ionization cross section of Fe xvii are described by Hahn
et al. (2013). A comparison of these measurements for Fe xvii, those of Dere
(2007) and the present results are shown in Fig. 6. The FAC calculations of
Dere (2007) only took into account direct ionization from the 2p shell. The
new fits to the measured cross sections are able to reproduce the measurements
quite well.

Figure 6: Ionization cross sections of Fe xvii. The label Hahn refers to Hahn
et al. (2013), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.6 The magnesium isoelectronic sequence: Fe XV

Fe xv produces a strong line at 284 Å in the EUV spectrum. In addition, it
also provides density-sensitive line pair ratios such as $\lambda$234 to
$\lambda$244 (Cowan & Widing, 1973). The ionization cross-section of Fe xv
consists of both DI and EA components and the Dere (2007) cross sections are
based on calculations with FAC. The EA components include transitions from the
n=2 levels to the n=3 and n=4 levels and the various transitions are grouped
into two EA components with excitation thresholds at 793 and 987 eV. These
compare to the ionization potential of 457 eV. Laboratory measurements of the
ionization cross sections of Fe xv are reported by Bernhardt et al. (2014) and
show both DI and EA components although the measured EA components appear more
complex than the two component description. The CHIANTI ionization cross
sections have been revised by fitting the DI component to the Bernhardt et al.
(2014) cross sections below 793 eV. The EA component to the n=3 levels has
been reduced by a factor of 0.92. The measurements of Bernhardt et al. (2014)
togther with the Dere (2007) values and the present results are shown in Fig.
7. At energies just below the first EA threshold, the ratio of the present
cross section to that of Dere (2007) is about 1.22. At energies where the EA
component is present, the ratio is about 0.99 at 1000 eV and 0.95 at 2000 eV.
Bernhardt et al. (2014) also provide parameters to reproduce their
calculations of the ionization rate coefficient. The rates of Bernhardt et al.
(2014) are very close to those of the revised parameters. The ratio of the
present rates to those of Bernhardt et al. (2014) is 0.83 at 105 K, 0.91 at
2.2 $\times$ 106 K (the temperature of maximum ionization of Fe xv) and 0.93
at 3 $\times$ 106 K.

Figure 7: Ionization cross sections of Fe xv. The label Bernhardt refers to
Bernhardt et al. (2014), 2007 refers to Dere (2007), and present to the
results reported here.

###  2.7 The aluminum isoelectronic sequence: Fe XIV

Fe xiv produces the strong “green” line at 5304 Å that is often observed with
solar coronagraphs. It is also responsible for a number of strong lines at EUV
wavelengths. Intensity ratios of some pairs of line are also sensitive
indicators of electron densities between 109 and 1011 cm-3 (Dere et al.,
1979).

As with Fe xv, the Fe xiv cross section includes both DI and EA components,
and both are prominent in both the measurements of Hahn et al. (2013) and the
FAC cross sections of Dere (2007). A comparison of the measurements of Hahn et
al. (2013), Dere (2007) and the present results are shown in Fig. 8. The
present values consist of a fit to the DI measured cross sections below 783 eV
where the EA components become apparent. The FAC cross sections of Dere (2007)
included excitations of a 2s electron to the n=3 and n=4 levels. The new EA
cross sections consist of adjustments to the EA threshold for the n=3
excitations and the magnitudes for both components are multiplied by a factor
of 0.8.

Figure 8: Ionization cross sections of Fe xiv. The label Hahn refers to Hahn
et al. (2013), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.8 The silicon isoelectronic sequence: Fe XIII

Fe xiii produces a number of strong spectral lines at EUV wavelengths. Some of
these line pairs are useful diagnostics of electron densities between 3
$\times$ 108 and 3 $\times$ 1010 cm-3 (Flower & Nussbaumer, 1974).

A comparison of the measurements of Hahn et al. (2011a, 2012b), Dere (2007)
and the present results are shown in Fig. 9. As with Fe xiv, the cross section
consists of strong DI and EA components. The DI fit parameters are arrived at
by using the measurements below the first EA threshold at 775 eV. The EA
components of the FAC cross sections included excitations from the n=2 levels
to the n = 3, 4, and 5 levels. A better fit was obtained for the present cross
section if the excitation to the n=5 levels was removed, the excitation to the
n=4 levels reduced by a factor of 0.8 and the transition to the n=3 levels
unchanged. At 400 eV, just above the threshold, the ratio of the present cross
section to that of the 2007 cross section is about 2.0, at 500 eV the ratio is
about 1.17, at 600 eV the ratio is about 1.05, and at 700 eV the ratio is
about 0.99. Above 780 eV where the EA cross sections occur, the ratio of the
present cross sections to the 2007 values is about 0.8. These ratios are
important for understanding the changes in the ionization rate coefficients
discussed in $\S$ 2.14.

Figure 9: Ionization cross sections of Fe xiii. The label Hahn refers to Hahn
et al. (2011a, 2012b), 2007 refers to Dere (2007), and present to the results
reported here.

###  2.9 The phosphorus isoelectronic sequence: Fe XII

Fe xii also produces a number of strong lines at EUV wavelengths. Some are
strong enough that they comprise the primary component of narrow band EUV
images such as the Extreme-ultraviolet Imaging Telescope (EIT) on SOHO and the
Atmospheric Imaging Assembly (AIA) on the Solar Dynamics Observatory (SDO). A
number of Fe xii line pairs form ratios that are reliable diagnostics of
electron densities between 108 and 1011 cm-3. They have been observed with the
XUV Slitless Spectrograph on Skylab (Dere et al., 1979) and with the Extreme-
ultraviolet Imaging Spectrometer (EIS) on the Hinode satellite (Dere et al.,
2007; Warren & Brooks, 2009).

A comparison of the measurements of Hahn et al. (2011b), Dere (2007) and the
present results are shown in Fig. 10. The DI FAC cross sections of Dere (2007)
took into account ionization from the n=2 levels and the 3s and 3p levels.
These have been replaced with a fit to the Hahn et al. (2011b) measurements
below the first EA threshold at 753 eV. The previous EA cross sections
considered excitations from the n=2 levels to the n=3, 4, and 5 levels. The
comparison in Fig. 10 indicates that measurements are best matched if the
excitations to the n=4 and 5 levels are removed. Otherwise, it has not been
necessary to change the energy of the n=3 threshold or magnitude.

Just above threshold the ratio of the present cross section to that of Dere
(2007) is about 2.5; at 400 eV, the ratio is about 1.4; and at 600 eV, it is
about 1.08. Above 1000 eV where the EA cross sections occur, the ratio is
about 0.8. Again, these ratios are important for understanding the changes in
the ionization rate coefficients discussed in $\S$ 2.14.

Figure 10: Ionization cross sections of Fe xii. The label Hahn refers to Hahn
et al. (2011b), 2007 refers to Dere (2007), and present to the results
reported here.

###  2.10 The sulphur isoelectronic sequence: Fe XI

The ion Fe xi contributes a number of strong lines to the EUV spectrum. It
also provides density sensitive line ratios that are useful in the density
range between 109 and 1011 cm-3 (Dere, 2020). A comparison of the measurements
of Hahn et al. (2012c), Dere (2007) and the present results are shown in Fig.
11. The FAC EA cross sections included excitations from the n=2 to the n=4 and
5 levels. The summed rates have been collected to form three EA cross
sections. The threshold energies have been adjusted to reproduce the measured
cross sections. For the 2007 EA cross sections, the lowest energy threshold
was 882 eV and this has been moved to 721 eV. The other two thresholds have
been adjusted by a similar amount. The magnitudes of the three original EA
cross sections have been scaled by factors of 1.3, 4.0 and 2.0. The ratio of
the present cross sections to the 2007 values have been evaluated. At 300 eV,
just above the threshold, the ratio is about 1.6, at 500 eV the ratio is 1.01
at 700 eV, just before the EA cross sections occur, the ratio is about 0.96.
At energies above 1000 eV, the ratio is about 0.97.

Figure 11: Ionization cross sections of Fe xi. The label Hahn refers to Hahn
et al. (2012c), 2007 refers to Dere (2007), and 2022 to the results reported
here.

###  2.11 The chlorine isoelectronic sequence: Fe X

Fe x produces a number of lines at EUV wavelengths as well as the coronal
“red” line at optical wavelengths. Some EUV line pairs form ratios that are
sensitive to electron densities between 3 $\times$108 and 1011 cm-3.

The ionization cross sections of Fe x have been measured by Hahn et al.
(2012c) and a comparison of these measurements, the FAC cross sections of Dere
(2007) and the present results are shown in Fig. 12. From this figure, it is
clear that both the DI and EA components of the original FAC calculations need
to be modified to accurately reproduce the measured values. The original FAC
calculations took into account DI ionization from the n=2 and n=3 levels. The
EA components included excitation from the n=2 levels to the n=4 and n=5
levels and excitations from the n=3 to the n=4, 5, and 6 levels. The latter
cross sections for excitations from the n=3 levels occur just above the
ionization threshold but are not very apparent in the total cross section. For
the DI component, a fit to the measurements has been obtained to the
measurements at energies below 650 eV. For the EA component, only the cross
sections from the n=2 levels are of importance and it has been necessary to
rescale the magnitudes and shift the energy thresholds. The energy threshold
for the lowest energy EA component has been moved from 862 eV to 720 eV and
the second component by a similar amount. The magnitudes of the two EA
components have been scaled by factors of 1.8 and 2.4. The ratio of the
present to the 2007 cross section is about 1.25 at an energy of 300 eV, just
above the DI threshold. At 500 eV the ratio is 1.06 and at energies above 1000
eV, the ratio is about 1.02.

Figure 12: Ionization cross sections of Fe x. The label Hahn refers to Hahn et
al. (2012c), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.12 The argon isoelectronic sequence: Fe IX

Fe ix produces a very strong line at 171 Å that is produced at about 8
$\times$ 105 K and shows characteristics of the solar upper transition region
and low corona. It comprises the primary component of narrow band EUV images
such as the EIT on SOHO and the AIA on SDO. A comparison of the measurements
of Hahn et al. (2016), Dere (2007) and the present results are shown in Fig.
13. The DI cross sections of Dere (2007) include ionization from the n=3
levels. These have been replaced by a fit to the measurements at energies
below about 700 eV and the result is that the DI cross section is reduced by
about 40 per cent. The EA components of the cross section in the original FAC
calculations included excitations from the n=2 and n=3 levels and these cross
sections were summed into three components with the lowest energy component at
842 eV. Hahn et al. (2016) found the dominant EA component to have a threshold
at about 650 eV and attributed it to excitations from the n=2 level to the 3d
level. This component was not included in the original FAC calculations. This
component is now modeled by a fit to the measurements above 700 eV with a
single EA component. At an energy of 250 eV, just above the ionization
potential of 233 eV, the ratio of the present cross section to the previous
FAC cross section is about 1.5. At 300 eV, the ratio is about 1.03, at 400 eV
the ratio is about 1.04, and at 600 eV, near the peak of the DI component, the
ratio is about 0.72. In the region above 700 eV, the ratio is about 0.85.

Figure 13: Ionization cross sections of Fe ix. The label Hahn refers to Hahn
et al. (2016), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.13 The potassium isoelectronic sequence: Fe VIII

Fe viii produces a number of lines at EUV wavelengths. These are generally not
particularly strong but do provide information on plasmas around 5 $\times$
105 K. Intensity ratios of several lines provide density diagnostics between
106 and 108 cm-3. A comparison of the measurements of Hahn et al. (2015), Dere
(2007) and the present results are shown in Fig. 14. Procedures were followed
during the experimental measurements to ensure that as few ions in metastable
states were present in the beam. The authors state that ions in metastable
states comprised only about 6% of the beam. The original FAC calculations
included DI from the 3p and 3d shells and EA excitations from the 3p63d to the
3p53d nl levels with n=4,5,6. The ionization potential is about 156 eV and the
excitations to the n=4,5 and 6 autoionizing levels occur just above the
ionization potential. The various EA excitations were grouped into five EA
cross sections with energies between 156 and 200 eV. These are evident as
step-like structures in the calculated ionization cross section just above
threshold. With the measurements of Hahn et al. (2015) now available, fits to
the measurements were used to represent the cross sections. These cross
sections are shown in Fig. 14.

Figure 14: Ionization cross sections of Fe viii. The label Hahn refers to Hahn
et al. (2015), 2007 refers to Dere (2007), and present to the results reported
here.

###  2.14 Ionization Rate Coefficients

The ionization rate coefficient, $R(T)$, is obtained by integrating the cross-
section over a Maxwell-Boltzmann distribution, and is given by

$R(T)=\int_{v_{IP}}^{\infty}v\,\sigma(E)\,f(v,T)\,dv$ (3)

where $v$ is the velocity of the electron, $E$ ($=\frac{1}{2}\,m\,v^{2}$) is
the energy of the electron, and $f$ is the Maxwell velocity distribution at a
temperature $T$. The integration is carried out from the velocity of the
electron with an energy equal to the ionization potential IP.

The Fe xiv ionization rate coefficients for the present rates and those of
Dere (2007) are shown in Fig. 15 together with the ratio of the rate
coefficients. The displayed temperature range approximately corresponds to
where the Fe xiv ionization fraction is within 0.01 of its peak value. It can
be seen that the revised rate coefficients do not show a significant
difference from the previous values even though there are considerable
differences in the ionization cross sections (see Fig. 8). The largest
difference is in the EA components that lie well above the ionization
potential. In order to understand this lack of change, the Maxwell-Boltzmann
distribution of electron velocities for a temperature $1.91\times 10^{6}$ K
has been plotted on Fig. 16. The ratios of the new rate coefficients to those
previously in CHIANTI for the ions in the present work are shown in Table 1.
In each case, the rate coefficients are calculated at the temperature of
maximum ionization.

Figure 15: Top: Ionization rate coefficients of Fe xiv. The label 2007 refers to Dere (2007), and 2022 to the present results. Bottom: The ratio of the present rate coefficients to those of 2007. Figure 16: Ionization cross sections of Fe xiv. The label Hahn refers to Hahn et al. (2013), 2007 refers to Dere (2007), present to the results reported here. Also plotted is the Maxwell-Boltzmann (MB) velocity distribution for a temperature of 1.91 MK and $\mathbf{m<v>^{2}/2}$ for the same temperature is displayed as a vertical line. Table 1: The ratio of the ionization rate coefficients reported here to those of Dere (2007) Ion |  $T$ (106 K) | Ratio
---|---|---
S xiii | 2.7 | 1.11
O v | 0.24 | 0.96
Mg viii | 0.79 | 1.09
Fe xviii | 7.9 | 1.09
Fe xvii | 5.6 | 1.06
Fe xv | 2.2 | 0.99
Fe xiv | 1.9 | 0.97
Fe xiii | 1.7 | 1.09
Fe xii | 1.6 | 1.17
Fe xi | 1.3 | 1.09
Fe x | 1.1 | 1.14
Fe ix | 0.79 | 0.72
Fe viii | 0.56 | 1.00

###  2.15 Recombination rate coefficients for the phosphorus isoelectronic
sequence

Recently Bleda et al. (2022) have reported radiative and dielectronic
recombination coefficients for ions in the phosphorus sequence for ions with
nuclear charge Z from 16 through 30 plus a few ions with higher nuclear
charges. These rates refer to the recombination of S-like ions to the P-like
ions. In the case of the recombination of Fe xii to Fe xi it was possible for
the authors to compare their calculations with the measurements of Novotný et
al. (2012). As the authors state, ”Given the complexity of the problem, …, the
overall agreement, …, is quite good.” Bleda et al. (2022) provide fit
parameters for their coefficients that are reproduced in the files released in
the current CHIANTI release.

In general, the changes from the rates in the version 10.0 CHIANTI database
are often not that large. In the case of S ii, the version 10.0 coefficients
were taken from the work of Mazzotta et al. (1998). The agreement between the
two is very good at high temperature but below $2\times 10^{4}$ K the rates of
Bleda et al. (2022) include the significant lower temperature processes.
However, the differences are mostly of importance for photoionized plasmas.
For Ca vi the new DR rates of Bleda et al. (2022) are about 50% lower than the
rates of Shull & van Steenberg (1982) that were included in CHIANTI version
10.0. For Ni xiv Bleda et al. (2022) show considerable difference between
their DR rates and those of Mazzotta et al. (1998). Version 10.0 of CHIANTI
still included the DR rates of Shull & van Steenberg (1982) and these are
quite close to the new rates of Bleda et al. (2022).

###  2.16 A revised ionization balance

With the revised set of ionization and recombination rate coefficients, an
updated ionization balance has been calculated and will be distributed as the
default ionization equilibrium for CHIANTI 10.1. Figure 17 compares the new
ionization fraction curves for the iron ions Fe x–xiv with those from CHIANTI
10. Comparing all of the ionization fraction curves between the new and
previous tabulations, there are 11 ions with $\log\,T_{\mathrm{max}}$ values
shifted by 0.05 dex and seven ions with peak ionization fractions that are
$>10$% different compared to the previous ionization fraction curves. The
latter include V viii,ix, Sc vi,vii, Ca vi and K v. The largest difference is
a 30% decrease for S i.

Figure 17: The revised ionization balance for Fe ions for stages X through
XIV. The full line represents the present calculation and the dashed line the
version 10.0 calculation

##  3 New and updated atomic models

###  3.1 Hydrogen isoelectronic sequence: C VI

An error in the proton rate coefficients used for C vi was identified by C. M.
Gunasekera and has been corrected. The rate coefficients for the transitions
between the three levels of the $n=2$ configuration were computed by Zygelman
& Dalgarno (1987) and added to CHIANTI Version 6 (Dere et al., 2009). Zygelman
& Dalgarno (1987) tabulated their rate coefficients with a scaling of the
atomic number, $Z$, to the power of three. In removing the scaling, the
numbers were inadvertently multiplied by $Z^{3}$ instead of divided, leading
to rate coefficients too large by a factor $4.7\times 10^{4}$.

The proton rates for C vi provide a depopulation mechanism for the $2s$
${}^{2}S_{1/2}$ metastable level in addition to the two-photon decay
mechanism. The erroneous rates resulted in proton excitation being the
dominant depopulation mechanism. With the correct rates, the two-photon
process becomes the dominant process. As a consequence, the two-photon
continuum for C vi is now more than a factor 1000 stronger than previously for
a temperature of 1 MK and an electron number density of $10^{10}$ cm-3.

A second effect is that the erroneous proton rates resulted in the populations
of the $2p$ ${}^{2}P_{1/2,3/2}$ levels being enhanced by 44% and 10%
respectively, for a temperature of 1 MK and an electron number density of
$10^{10}$ cm-3. The Ly-$\alpha$ line for C vi, which is a self-blend of the
transitions from the ${}^{2}P$ levels to $1s$ ${}^{2}S_{1/2}$, is consequently
18% weaker in the new C vi model compared to CHIANTI 10.

###  3.2 Boron isoelectronic sequence: O IV

Errors were found in the energy levels for this ion such that experimental
energies had been inadvertently assigned to the wrong levels. In particular,
the energies for levels 27, 28 and 51 had been assigned to levels 29, 30 and
52, respectively, and vice versa. The energies have been updated and a new set
of wavelengths have been derived from the updated energies. No other changes
were made to the O iv model.

###  3.3 Carbon isoelectronic sequence

####  3.3.1 Na VI

The CHIANTI 10 energies were obtained from NIST which has the data compiled by
Sansonetti (2008). As the $2s2p^{3}$ ${}^{5}S_{2}$ energy was obtained by
isoelectronic fitting performed by Edlén (1985), there is an uncertainty that
propagates to all quintet levels in the ion, even though many give rise to
experimental wavelengths. For the CHIANTI 10 model, none of the quintet levels
were assigned an experimental energy for this reason. For the present release,
the $2s2p^{3}$ ${}^{5}S_{2}$ energy has been assigned, allowing a number of
additional assignments for other quintet levels.

Feldman et al. (1997) identified the ${}^{3}P_{2}$–${}^{5}S_{2}$ transition in
solar spectra at 988.71 Å, which implies a ${}^{5}S_{2}$ energy of
$102\,998\pm 3.0$ cm-1. Sansonetti (2008) rejected this value because Feldman
et al. (1997) did not also identify the ${}^{3}P_{1}$–${}^{5}S_{2}$
transition, which is the stronger of the two transitions. However, this line
is predicted to be at 977.48 Å (or 977.37 Å if using the Edlén energy),
placing it close to C iii 977.02 Å. This is one of the strongest lines in the
solar spectrum making it difficult to identify the Na vi line. Given that the
${}^{5}S_{2}$ energy derived from the Feldman et al. (1997) wavelength is only
12 cm-1 from the Edlén estimated value we add it to CHIANTI since it accounts
for a known solar emission line. As a consequence, all of the quintet level
energies marked with “+x” in the NIST table were adjusted and added to
CHIANTI. There are seven levels in all: $3s$ ${}^{5}P_{J}$ (indices 35–37),
$3p$ ${}^{5}P_{J}$ (indices 59–61) and $4s$ ${}^{5}P_{2}$ (index 189). The
theoretical values of $4s$ ${}^{5}P_{0,1}$ (indices 187 and 188) were adjusted
based on the experimental energy of ${}^{5}P_{2}$.

Finally, the experimental energies for the $2s^{2}2p^{2}$ ${}^{1}D_{2}$ and
${}^{1}S_{0}$ levels were updated to the values from Young et al. (2011).

####  3.3.2 Ca XV

Updated atomic models for the entire carbon isoelectronic sequence were added
to CHIANTI 10. However, it was found that the Ca xv effective collision
strengths were omitted from the update. This has been rectified for CHIANTI
10.1.

###  3.4 Nitrogen isoelectronic sequence

Mao et al. (2020) presented a large-scale scattering calculation using the
$R$-matrix intermediate coupling frame transformation (ICFT) method, to
calculate effective collision strengths of N-like ions from O ii to Zn xxiv.
725 fine-structure levels were included in both the configuration interaction
target and close-coupling collision expansion, for states up to $n=5$. We are
providing new models based on these calculations for the O ii, Si viii, Ar xii
and Ca xiv ions. As shown in Mao et al. (2020), the new rates for these ions
are significantly different than those previously available in CHIANTI, and
they affect some line ratio diagnostics.

For the experimental energies we have used a combination of NIST values and
those that were assessed in the previous CHIANTI versions. Only bound states
have been retained for O II and Si VIII, as the model of the autoionizing
states requires significant additions and these ions do not produce strong
satellite lines.

For O ii, we have replaced the effective collision strengths for the first
five states of the ground configuration with those from Tayal (2007), up to
105 K, as we consider them to be more accurate at low temperatures. We note
that the values of the two calculations at 105 K differ by only a few percent,
as shown in Mao et al. (2020).

The structure calculations carried out with autostructure (AS) for the
scattering work are generally not as accurate as larger-scale ones, and
generally the radiative data for the lower states can be improved, either by
using experimental energies or replacing them with more accurate values. We
have compared the $A$-values with different calculations and found generally
small variations of the order of 10–20% for the important transitions within
the lowest states.

For O ii, Si viii, and Ar xii, we have replaced the $A$-values of the
transitions within the 2s2 2p3, 2s 2p4 configurations with the non-
relativistic multiconfiguration Hartree-Fock (MCHF) calculations of Tachiev &
Froese Fischer (2002). For Ca xiv, we have replaced the $A$-values of the
transitions within the 2s2 2p3, 2s 2p4 configurations with the many-body
perturbation theory (MBPT) calculations by Wang et al. (2016). The AS
radiative data only differed by typically 10–20%, compared to the MCHF or MBPT
values.

###  3.5 Oxygen isoelectronic sequence

Mao et al. (2021) performed a large-scale scattering calculation using the
$R$-matrix ICFT method, to calculate effective collision strengths of O-like
ions from Ne III to Zn XXIII over a wide range of temperatures. The targets
included 630 fine-structure levels up to $nl=5d$ for each ion. As discussed in
Mao et al. (2021), significant differences with the data that are present in
CHIANTI were found for a few ions important for astrophysical applications.
They include Si vii, S ix, Ar xi, and Ca xiii, for which CHIANTI had limited
models and effective collision strengths calculated with the distorted wave
approximation. Lines in the EUV, UV and near infrared have significantly
different predicted intensities with the new atomic data. For this reason, we
have included new models for these ions. Work is in progress to update the
models for the other ions, which will be made available in a future CHIANTI
version.

For the experimental energies we have used a combination of NIST values and
those that were assessed in the previous CHIANTI versions. Only bound states
have been retained, as the model of the autoionizing states requires
significant additions.

For Si vii and S ix we have replaced the A-values of the transitions within
the $n=2$ complex (ground configuration 2s2 2p4, 2s 2p5 and 2p6, ten lowest
states) with those calculated with the Multiconfiguration Hartree-Fock (MCHF)
codes by Tachiev & Froese Fischer (2002), noting that the Mao et al. (2021)
A-values differed by typically 10–20% only.

For Ar xi and Ca xiii we have replaced the AS values with those calculated by
Song et al. (2021) with the Multiconfiguration Dirac-Hartree-Fock (MCDHF)
GRASP2K codes, described in Jönsson et al. (2013). Differences with the AS
values are about 10%.

###  3.6 Silicon isoelectronic sequence: Fe XIII

Zhang et al. (2021) carried out a large-scale MCDHF atomic structure
calculation for Fe xiii, providing theoretical energies that are generally
close to the experimental energies, typically within a few hundreds of cm-1.
Zhang et al. (2021) reviewed previous studies and confirm the new
identifications suggested by Del Zanna (2012) (based on the Del Zanna (2011)
calculations) for the $n=4\to n=3$ transitions, with one exception. We adopt
the Zhang et al. (2021) identification. The $n=4\to n=3$ transitions fall in
the soft X-ray region of 50 Å – 150 Å.

A few weak EUV transitions from the 3s 3p2 3d levels were tentatively
identified by Del Zanna (2011) and assigned experimental wavelengths, for
their potential importance in analysing Hinode EIS observations. However, the
analysis carried out by Zhang et al. (2021) provides evidence for alternative
identifications. Some of them have been adopted for the present version, while
other tentative identifications have not been included. However, we are
providing as ‘best’ energies those from Zhang et al., hence the expected
wavelengths should be close to the experimental ones.

###  3.7 Argon isoelectronic sequence: Fe IX

Ryabtsev et al. (2022) provided a new set of experimental energies for the
$3p^{4}3d^{2}$ and $3p^{5}4f$ configurations and these have been added to
CHIANTI with the modifications described below.

The CHIANTI 10 model had experimental energies for all 12 of the $3p^{5}4f$
fine-structure levels from O’Dwyer et al. (2012), and these have been updated
with the values from Ryabtsev et al. (2022). Agreement is good for nine of the
levels, but larger differences of between 160 and 560 cm-1 are found for the
remaining levels due to differences in line identifications—see Ryabtsev et
al. (2022) for more details.

Comparisons of the $A$-values listed by Ryabtsev et al. (2022) with those in
the CHIANTI model (from Del Zanna et al., 2014) showed that it was necessary
to swap the Ryabtsev et al. (2022) labels for the ${}^{1}F_{3}$ and
${}^{3}F_{3}$ levels, and the ${}^{1}G_{4}$ and ${}^{3}F_{4}$ levels to get
agreement. Both sets of authors found that these levels are strongly mixed,
hence the $LSJ$ labels are not accurate descriptors for them. The Del Zanna et
al. (2014) and Ryabtsev et al. (2022) $A$-values are in very good agreement
after this modification.

The $3p^{4}3d^{2}$ configuration has 111 fine-structure levels, but only five
previously had experimental energies. Sixteen new experimental energies have
been taken from Ryabtsev et al. (2022), and the five existing energies have
been updated, with only small changes. Ryabtsev et al. (2022) provided
experimental energies for nine additional levels, but three of these were not
used as they were listed as ‘questionable.’ Six further energies were not used
either because the identifications were not consistent with predictions from
the CHIANTI atomic model, or because the identifications were uncertain due to
line blending. Full details of how these decisions were made are given in
Young (2023a).

Young & Landi (2009) performed a study of Fe ix lines observed by the Extreme
Ultraviolet Imaging Spectrometer (EIS Culhane et al., 2007), and they
identified seven lines that could be attributed to Fe ix but for which it was
not possible to assign atomic transitions. Two of these lines were
independently identified as Fe ix lines by Del Zanna (2009a), but also without
atomic transitions. Del Zanna et al. (2014) provided transition information
for two lines, and Ryabtsev et al. (2022) identified six of them. The two
works agree for the line at 192.63 Å, but differ for the line at 194.80 Å. The
Ryabtsev et al. (2022) identification is supported by multiple lines from the
same upper level that are measured in laboratory spectra.

The CHIANTI wavelengths in the present version have been derived from the new
energy levels.

###  3.8 Calcium isoelectronic sequence: Fe VII

Experimental energies have been updated using the compilation of Kramida et
al. (2022) that provided energies for previously unassigned levels and updated
energies of known levels. A number of levels are highly mixed, particularly
for the $3p^{5}3d^{3}$ configuration and it was not always possible to match
the Kramida et al. (2022) level with a level in the existing CHIANTI model. In
these cases the Kramida et al. (2022) level was not used. Details on which
level energies were not chosen for the CHIANTI model and how the decisions
were made are described in Young (2023b). In summary, the new model has 17
levels with newly-assigned experimental energies, and 10 levels for which the
updated energy is at least 100 cm-1 different from the previous energy. The
wavelengths in the present version have been updated with the new energies.
However, as pointed out in the previous CHIANTI release paper, some
inconsistencies remain between the predicted intensities of the strongest
lines and those observed by Hinode EIS. Del Zanna (2009b) used an alternative
set of atomic data and identifications to find relatively good agreement with
observations. As different calculations provide very different energies and
rates (even branching ratios), the solution of this complex problem will
require further calculations and further assessments.

##  4 Elemental Abundances

The CHIANTI software uses solar photospheric elemental abundances as the
default when computing synthetic spectra. For CHIANTI versions 9 and 10, the
photospheric abundances were stored in the file
sun_photospheric_2015_scott.abund. This contains the abundances of Asplund et
al. (2009), supplemented with values from Scott et al. (2015a, b) and Grevesse
et al. (2015) for some elements.

For CHIANTI 10.1 the default abundance file—named
sun_photospheric_2021_asplund.abund—contains data from the compilation of
Asplund et al. (2021). These authors provide a comprehensive and detailed
review of the photospheric abundance of all elements up to uranium,
complementing recent re-evaluations and new measurements in the literature for
the abundances of many elements with original results. In both cases, improved
values were obtained mostly from the coupling of photospheric observed spectra
with 3D, NLTE hydrodynamical simulations of the outer convective zone and
atmospheric layers of the Sun, and improved atomic data and collisional rates
for many species. Additional data were taken from laboratory analysis of
meteorites, return samples from the Genesis mission, helioseismological data,
and sunspot observations.

In comparison to the previous default abundance set, the Asplund et al. (2021)
abundances show differences of 10% or more for four elements: lithium
($-19$%), neon ($+35$%), chlorine ($-19$%) and titanium ($+10$%). The most
significant for the spectral modeling of ionized plasmas is the neon
abundance.

The new dataset provides a Ne/O relative abundance ratio of 0.23, larger than
the previous value of 0.17. The new value is in line with recent spectroscopic
re-evaluations of the Ne/O abundance ratio in the transition region (Young,
2018) and in the solar corona during solar minimum (Landi & Testa, 2015). The
Ne/O ratio, and the absolute abundances of both Ne and O are very important
for solving the current discrepancy between helioseismology determinations of
the structure of the solar interior with model predictions, that need accurate
absolute abundance measurements to calculate solar interior opacities
(Christensen-Dalsgaard, 2021).

We finally note that the coronal abundance data sets provided in earlier
CHIANTI versions (e.g., sun_coronal_1992_feldman.abund and
sun_coronal_2012_schmelz.abund) were largely based on solar active region
observations where the relative abundances of low vs. high first ionization
potential (FIP) elements showed an increase (the so-called FIP bias), compared
to their (relative) photospheric values (see, e.g. the reviews by Laming,
2015; Del Zanna & Mason, 2018). The FIP bias was found to have averaged values
around 3–4. Those coronal abundance data sets were obtained by taking into
account those averaged values and applying empirical corrections (the FIP
bias) to the older photospheric abundance data sets. The previous coronal
abundance datasets are therefore not consistent with the new photospheric
abundance file discussed here. We have therefore introduced a new coronal
abundance file sun_coronal_2021_chianti.abund that is derived from the new
photospheric abundance file by multiplying abundances for low-FIP elements
(FIP $\leq$ 10 eV) by a factor $10^{0.5}$. This file is provided as a
representative coronal abundance file that may be useful for investigating how
abundances affect synthetic spectra, for example. Users should be aware that
the magnitude of the FIP bias has been found to vary amongst different solar
structures, and also with temperature.

Abundance files (both photospheric and coronal) previously distributed with
CHIANTI will continue to be available but are moved to the abundance/archive
directory within the database.

##  5 Conclusions

This article has described the most recent updates to the CHIANTI atomic
database that will be distributed as version 10.1. Ionization rates for 13
ions have been updated, as well as the complete set of dielectronic
recombination data for the phosphorus isoelectronic sequence. The core
datasets of energy levels, radiative decay rates and electron excitation rates
have been updated for 15 ions.

Dere (2007) assembled the experimental measurements of ionization cross
sections available at that time. These were supplemented with calculations of
the ionization cross sections with the FAC package (Gu, 2002) where
measurements did not exist or were unreliable. Together with recombination
rates from the scientific literature, version 9.0 of the CHIANTI atomic
database was assembled with a full set of rates covering all ions of all
elements from H through Zn. Since that time, there have been a number of new
experimental measurements of ionization cross sections, as well as new
calculations of recombination rates. The new ionization measurements have been
used to revised the CHIANTI parameterization of these cross sections and these
have been inserted into the latest version of the CHIANTI atomic database
(Version 10.1). In most cases, the differences between the revised rates and
the previous CHIANTI rates is not that great. The parameters for recombination
of Bleda et al. (2022) have also been inserted into the latest version of
CHIANTI. In addition, a new ionization equilibrium has been calculated from
these revised rates.

A number of atomic models have been updated. For the nitrogen isoelectronic
sequence, new models for the O ii, Si viii, Ar xii and Ca xiv ions are now
based on the recent calculations of Mao et al. (2020). For the oxygen
isoelectronic sequence, the calculations of Mao et al. (2021) have allowed us
to improve the previously limited models for Si vii, S ix, Ar xi, and Ca xiii.

All of these are now included in the latest release of the CHIANTI atomic
database version 10.1.

We thank Dr. Stefan Schippers for supplying the measurements of Bleda et al.
(2022) in a machine-readable format. KPD and PRY work has been supported by
NASA grants 80NSSC21K0110 and 80NSSC21K1785. GDZ acknowledges support from
STFC (UK) via the consolidated grants to the atomic astrophysics group (AAG)
at DAMTP, University of Cambridge (ST/P000665/1. and ST/T000481/1). EL has
been supported by NASA grants 80NSSC22K0750 and 80NSSC20K0185.

## References

  * Asplund et al. (2021) Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141, doi: 10.1051/0004-6361/202140445
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Badnell (2011) Badnell, N. R. 2011, Comput. Phys. Commun., 182, 1528, doi: 10.1016/j.cpc.2011.03.023
  * Badnell et al. (2003) Badnell, N. R., O’Mullane, M. G., Summers, H. P., et al. 2003, A&A, 406, 1151, doi: 10.1051/0004-6361:20030816
  * Bernhardt et al. (2014) Bernhardt, D., Becker, A., Grieser, M., et al. 2014, Phys. Rev. A, 90, 012702, doi: 10.1103/PhysRevA.90.012702
  * Bleda et al. (2022) Bleda, E. A., Altun, Z., & Badnell, N. R. 2022, A&A, 668, A72, doi: 10.1051/0004-6361/202244043
  * Burgess & Tully (1992) Burgess, A., & Tully, J. A. 1992, A&A, 254, 436
  * Christensen-Dalsgaard (2021) Christensen-Dalsgaard, J. 2021, Living Reviews in Solar Physics, 18, 2, doi: 10.1007/s41116-020-00028-3
  * Cowan & Widing (1973) Cowan, R. D., & Widing, K. G. 1973, ApJ, 180, 285, doi: 10.1086/151962
  * Crandall et al. (1979) Crandall, D. H., Phaneuf, R. A., Hasselquist, B. E., & Gregory, D. C. 1979, Journal of Physics B Atomic Molecular Physics, 12, L249, doi: 10.1088/0022-3700/12/7/008
  * Culhane et al. (2007) Culhane, J. L., Harra, L. K., James, A. M., et al. 2007, Sol. Phys., 243, 19, doi: 10.1007/s01007-007-0293-1
  * Del Zanna (2009a) Del Zanna, G. 2009a, A&A, 508, 501, doi: 10.1051/0004-6361/200913082
  * Del Zanna (2009b) —. 2009b, A&A, 508, 501, doi: 10.1051/0004-6361/200913082
  * Del Zanna (2011) —. 2011, A&A, 533, A12, doi: 10.1051/0004-6361/201117121
  * Del Zanna (2012) —. 2012, A&A, 546, A97, doi: 10.1051/0004-6361/201219923
  * Del Zanna et al. (2021) Del Zanna, G., Dere, K. P., Young, P. R., & Landi, E. 2021, ApJ, 909, 38, doi: 10.3847/1538-4357/abd8ce
  * Del Zanna & Ishikawa (2009) Del Zanna, G., & Ishikawa, Y. 2009, A&A, 508, 1517, doi: 10.1051/0004-6361/200911729
  * Del Zanna & Mason (2018) Del Zanna, G., & Mason, H. E. 2018, Living Reviews in Solar Physics, 15, 5, doi: 10.1007/s41116-018-0015-3
  * Del Zanna et al. (2014) Del Zanna, G., Storey, P. J., Badnell, N. R., & Mason, H. E. 2014, A&A, 565, A77, doi: 10.1051/0004-6361/201323297
  * Del Zanna & Young (2020) Del Zanna, G., & Young, P. R. 2020, Atoms, 8, 46, doi: 10.3390/atoms8030046
  * Dere (2007) Dere, K. P. 2007, A&A, 466, 771, doi: 10.1051/0004-6361:20066728
  * Dere (2020) —. 2020, MNRAS, 496, 2334, doi: 10.1093/mnras/staa1645
  * Dere et al. (2007) Dere, K. P., Doschek, G. A., Mariska, J. T., et al. 2007, PASJ, 59, S721, doi: 10.1093/pasj/59.sp3.S721
  * Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149, doi: 10.1051/aas:1997368
  * Dere et al. (2009) Dere, K. P., Landi, E., Young, P. R., et al. 2009, A&A, 498, 915, doi: 10.1051/0004-6361/200911712
  * Dere et al. (1979) Dere, K. P., Mason, H. E., Widing, K. G., & Bhatia, A. K. 1979, ApJS, 40, 341, doi: 10.1086/190590
  * Dufresne & Del Zanna (2019) Dufresne, R. P., & Del Zanna, G. 2019, A&A, 626, A123, doi: 10.1051/0004-6361/201935133
  * Dufresne et al. (2020) Dufresne, R. P., Del Zanna, G., & Badnell, N. R. 2020, MNRAS, 497, 1443, doi: 10.1093/mnras/staa2005
  * Edlén (1985) Edlén, B. 1985, Phys. Scr, 31, 345, doi: 10.1088/0031-8949/31/5/008
  * Falk et al. (1983) Falk, R. A., Stefani, G., Camilloni, R., et al. 1983, Phys. Rev. A, 28, 91, doi: 10.1103/PhysRevA.28.91
  * Feldman et al. (1997) Feldman, U., Behring, W. E., Curdt, W., et al. 1997, ApJS, 113, 195, doi: 10.1086/313048
  * Flower & Nussbaumer (1974) Flower, D. R., & Nussbaumer, H. 1974, A&A, 31, 353
  * Fogle et al. (2008) Fogle, M., Bahati, E. M., Bannister, M. E., et al. 2008, ApJS, 175, 543, doi: 10.1086/525256
  * Grevesse et al. (2015) Grevesse, N., Scott, P., Asplund, M., & Sauval, A. J. 2015, A&A, 573, A27, doi: 10.1051/0004-6361/201424111
  * Gu (2002) Gu, M. F. 2002, ApJ, 579, L103
  * Hahn et al. (2010) Hahn, M., Bernhardt, D., Lestinsky, M., et al. 2010, ApJ, 712, 1166, doi: 10.1088/0004-637X/712/2/1166
  * Hahn et al. (2011a) Hahn, M., Grieser, M., Krantz, C., et al. 2011a, ApJ, 735, 105, doi: 10.1088/0004-637X/735/2/105
  * Hahn et al. (2011b) Hahn, M., Bernhardt, D., Grieser, M., et al. 2011b, ApJ, 729, 76, doi: 10.1088/0004-637X/729/1/76
  * Hahn et al. (2012a) —. 2012a, Phys. Rev. A, 85, 042713, doi: 10.1103/PhysRevA.85.042713
  * Hahn et al. (2012b) Hahn, M., Grieser, M., Krantz, C., et al. 2012b, ApJ, 761, 79, doi: 10.1088/0004-637X/761/1/79
  * Hahn et al. (2012c) Hahn, M., Becker, A., Grieser, M., et al. 2012c, ApJ, 760, 80, doi: 10.1088/0004-637X/760/1/80
  * Hahn et al. (2013) Hahn, M., Becker, A., Bernhardt, D., et al. 2013, ApJ, 767, 47, doi: 10.1088/0004-637X/767/1/47
  * Hahn et al. (2015) —. 2015, ApJ, 813, 16, doi: 10.1088/0004-637X/813/1/16
  * Hahn et al. (2016) —. 2016, Journal of Physics B Atomic Molecular Physics, 49, 084006, doi: 10.1088/0953-4075/49/8/084006
  * Jönsson et al. (2013) Jönsson, P., Gaigalas, G., Bieroń, J., Fischer, C. F., & Grant, I. 2013, Computer Physics Communications, 184, 2197, doi: https://doi.org/10.1016/j.cpc.2013.02.016
  * Kramida et al. (2022) Kramida, A., Ryabtsev, A. N., & Young, P. R. 2022, ApJS, 258, 37, doi: 10.3847/1538-4365/ac3a7e
  * Laming (2015) Laming, J. M. 2015, Living Reviews in Solar Physics, 12, 2, doi: 10.1007/lrsp-2015-2
  * Landi & Testa (2015) Landi, E., & Testa, P. 2015, ApJ, 800, 110, doi: 10.1088/0004-637X/800/2/110
  * Loch et al. (2003) Loch, S., Colgan, J., Pindzola, M., et al. 2003, Physical Review A, 67, 042714
  * Mao et al. (2020) Mao, J., Badnell, N. R., & Del Zanna, G. 2020, A&A, 643, A95, doi: 10.1051/0004-6361/202039195
  * Mao et al. (2021) —. 2021, A&A, 653, A81, doi: 10.1051/0004-6361/202141464
  * Mazzotta et al. (1998) Mazzotta, P., Mazzitelli, G., Colafrancesco, S., & Vittorio, N. 1998, A&AS, 133, 403, doi: 10.1051/aas:1998330
  * Novotný et al. (2012) Novotný, O., Badnell, N. R., Bernhardt, D., et al. 2012, ApJ, 753, 57, doi: 10.1088/0004-637X/753/1/57
  * O’Dwyer et al. (2012) O’Dwyer, B., Del Zanna, G., Badnell, N. R., Mason, H. E., & Storey, P. J. 2012, A&A, 537, A22, doi: 10.1051/0004-6361/201117926
  * Ryabtsev et al. (2022) Ryabtsev, A. N., Kononov, E. Y., & Young, P. R. 2022, ApJ, 936, 60, doi: 10.3847/1538-4357/ac7d51
  * Sansonetti (2008) Sansonetti, J. E. 2008, Journal of Physical and Chemical Reference Data, 37, 1659, doi: 10.1063/1.2943652
  * Scott et al. (2015a) Scott, P., Asplund, M., Grevesse, N., Bergemann, M., & Sauval, A. J. 2015a, A&A, 573, A26, doi: 10.1051/0004-6361/201424110
  * Scott et al. (2015b) Scott, P., Grevesse, N., Asplund, M., et al. 2015b, A&A, 573, A25, doi: 10.1051/0004-6361/201424109
  * Shull & van Steenberg (1982) Shull, J. M., & van Steenberg, M. 1982, ApJS, 48, 95, doi: 10.1086/190769
  * Song et al. (2021) Song, C. X., Zhang, C. Y., Wang, K., et al. 2021, Atomic Data and Nuclear Data Tables, 138, 101377, doi: 10.1016/j.adt.2020.101377
  * Tachiev & Froese Fischer (2002) Tachiev, G. I., & Froese Fischer, C. 2002, A&A, 385, 716, doi: 10.1051/0004-6361:20011816
  * Tayal (2007) Tayal, S. S. 2007, ApJS, 171, 331, doi: 10.1086/513107
  * Thompson & Gregory (1994) Thompson, J. S., & Gregory, D. C. 1994, Phys. Rev. A, 50, 1377, doi: 10.1103/PhysRevA.50.1377
  * Wang et al. (2016) Wang, K., Si, R., Dang, W., et al. 2016, ApJS, 223, 3, doi: 10.3847/0067-0049/223/1/3
  * Warren & Brooks (2009) Warren, H. P., & Brooks, D. H. 2009, ApJ, 700, 762, doi: 10.1088/0004-637X/700/1/762
  * Warren et al. (2008) Warren, H. P., Feldman, U., & Brown, C. M. 2008, ApJ, 685, 1277, doi: 10.1086/591075
  * Young (2023a) Young, P. 2023a, Updating Fe IX energy levels for CHIANTI 10.1, doi: 10.5281/zenodo.7803672
  * Young (2023b) —. 2023b, Updating Fe VII energy levels for CHIANTI 10.1, doi: 10.5281/zenodo.7799540
  * Young (2018) Young, P. R. 2018, ApJ, 855, 15, doi: 10.3847/1538-4357/aaab48
  * Young et al. (2016) Young, P. R., Dere, K. P., Landi, E., Del Zanna, G., & Mason, H. E. 2016, Journal of Physics B Atomic Molecular Physics, 49, 074009, doi: 10.1088/0953-4075/49/7/074009
  * Young et al. (2011) Young, P. R., Feldman, U., & Lobel, A. 2011, ApJS, 196, 23, doi: 10.1088/0067-0049/196/2/23
  * Young & Landi (2009) Young, P. R., & Landi, E. 2009, ApJ, 707, 173, doi: 10.1088/0004-637X/707/1/173
  * Zhang et al. (2021) Zhang, X. H., Del Zanna, G., Wang, K., et al. 2021, ApJS, 257, 56, doi: 10.3847/1538-4365/ac2a3f
  * Zygelman & Dalgarno (1987) Zygelman, B., & Dalgarno, A. 1987, Phys. Rev. A, 35, 4085, doi: 10.1103/PhysRevA.35.4085


# Quasinormal modes in higher-derivative gravity: Testing the black hole
parametrization and sensitivity of overtones

R. A. Konoplya<EMAIL_ADDRESS>Research Centre for Theoretical
Physics and Astrophysics,
Institute of Physics, Silesian University in Opava,
Bezručovo náměstí 13, CZ-74601 Opava, Czech Republic

###### Abstract

The fundamental quasinormal modes of black holes in higher-derivative gravity
given by the Einstein-Weyl action are known to be moderately corrected by the
Weyl term. Here we will show that the first several overtones are highly
sensitive to even a relatively small Weyl correction, which might be important
when representing the earlier stage of the black hole ringdown. In addition,
we have solved the problem related to the analytical parametrized
approximation of the numerical black hole solution in the Einstein-Weyl
theory: In some range of parameters the approximation for the metric developed
up to the third order leads to the unusual highly nonmonotonic behavior of the
frequencies. We have shown that this problem can be solved via the extension
of the parametrization of the metric to higher orders until reaching the
regime when the frequencies do not change with further increasing of the
order.

###### pacs:

04.50.Kd,04.70.-s

##  I Introduction

Proper oscillation frequencies of black holes called quasinormal modes
Kokkotas:1999bd ; Berti:2009kk ; Konoplya:2011qq  have been the focus of
theoretical and experimental studies during the past few decades. The
observation of gravitational waves from the merger of two black holes by the
LIGO and VIRGO Collaborations Abbott:2016blz  and a series of subsequent
observations LIGOScientific:2020ibl ; KAGRA:2013rdx  agree with General
Relativity. Nevertheless, the angular momenta and masses of resultant black
holes are determined with a large uncertainty of tens of percent. This leaves
a wide window for alternative theories of gravity Konoplya:2016pmh ;
Yunes:2016jcc ; Berti:2018vdi . That is why quasinormal modes of the four-
dimensional asymptotically flat black holes in various alternative theories of
gravity have been intensively studied.

Among various alternative theories of gravity those with higher curvature
corrections are important, because they are predicted by the low-energy limit
of string theory. One such theory is represented by the Einstein action
corrected by the Weyl term. Then, the general vacuum four-dimensional theory
up the second order in curvature is described by the following Lagrangian:

$\displaystyle L$ $\displaystyle=$ $\displaystyle\sqrt{-g}(\gamma R-\alpha
C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}+\beta R^{2}),$ (1)

where $g$ is determinant of the matrix of a metric tensor; $\alpha$, $\beta$
and $\gamma$ are constants; and $C_{\mu\nu\rho\sigma}$ is the Weyl tensor. For
spherically symmetric and asymptotically flat solutions we can choose
$\gamma=1$ and $\beta=0$ without loss of generality, so that only the coupling
constant $\alpha$ matters for nonrotating black holes. The Schwarzschild
solution is a solution of the Einstein-Weyl equations, yet in Lu:2015cqa , it
was shown that a non-Schwarzschild asymptotically flat solution representing a
black hole is possible and a numerical solution for it was obtained. Using the
general parametrization for asymptotically flat black hole spacetimes in
arbitrary metric theories of gravity Rezzolla:2014mua ; Konoplya:2016jvv ;
Younsi:2016azx , the numerical black hole solution was approximated by an
analytical expressions for the metric functions in Ref. Kokkotas:2017zwt .
This theory and the corresponding black hole solutions were studied in a
number of works Einstein-Weyl:2018pfe ; QuadraticGravity:2016kip ;
Zinhailo:2019rwd ; Zinhailo:2018ska ; Konoplya:2019ppy ; Wu:2019uvq .

The asymptotically flat non-Schwarzschild solution in the Einstein-Weyl theory
has an interesting property: when the Weyl coupling constant increases
approaching the extremal limit (beyond which only the Schwarzschild solution
survives), the black hole mass decreases, tending to zero and having the
entire time, a nonvanishing radius. Therefore, the near-extremal non-
Schwarzschild solution may be a “light black hole,” which is unlikely to be
created in the course of the gravitational collapse of stars.

The low-lying quasinormal modes of non-Schwarzschild black holes in the
Einstein-Weyl theory have been studied in a few works. The eikonal formula for
the quasinormal modes of test fields was obtained in Ref. Kokkotas:2017zwt .
Massless scalar field perturbations were considered in Ref. Cai:2015fia , but
as was shown in Refs. Zinhailo:2019rwd ; Zinhailo:2018ska  the analysis based
on the Wentzel-Kramers-Brillouin (WKB) approach and time-domain integration of
the numerical solution apparently suffers from very large numerical
inaccuracy, making the final results doubtful. Indeed, Table I in Ref.
Cai:2015fia  shows that the difference between the WKB and time-domain
integration results reaches approximately $70\%$. The problem consists in the
insufficient precision of the numerical metric function when further taking
higher derivatives and integrating the wave equation. The obtained frequencies
have random, noisy-like dependence on the coupling parameter $\alpha$ (as was
shown in Fig. 1 of Ref. Zinhailo:2018ska ), which turns to be far from the
much more accurate values found later in Refs. Zinhailo:2019rwd ;
Zinhailo:2018ska . In Refs. Zinhailo:2019rwd ; Zinhailo:2018ska  the problem
was solved using the analytical approximation for the metric function
Kokkotas:2017zwt . There, the scalar, electromagnetic and Dirac low-lying
quasinormal modes were studied in detail. However, despite the quasinormal
frequencies this time depended smoothly upon the coupling parameter $\alpha$,
one technical problem remained: at some value of $\alpha$ a highly
nonmonotonic dependence, a kind of "splash," was observed for the Dirac field
perturbations at lower multipole numbers $\ell$ (see Fig. 5(b) in Ref.
Zinhailo:2019rwd ). In the present paper we will show that such strange
behavior is not any kind of numerical inaccuracy of the calculation method for
quasinormal modes, but is related to bad accuracy of the parametrization at
those values of $\alpha$. The higher-order parametrization remedies this
problem.

However, our main interest here is related to the first several overtones
which were not studied previously, because only the Wentzel-Kramers-Brillouin
or time-domain integration methods were used for the analysis of the spectrum
of the Einstein-Weyls black holes Cai:2015fia ; Zinhailo:2019rwd ;
Zinhailo:2018ska . It may seem that the fundamental mode provides sufficient
information about the signal; however, recently it has been shown that while
the lowest mode may not feel the deformations near the event horizon of the
black holes, the first few overtones may be highly sensitive to the tiny
changes of the near-horizon geometry Konoplya:2022pbc ; Konoplya:2022hll . An
important recent finding was made in Ref. Giesler:2019uxc  (and later studied
in Refs. Oshita:2021iyn ; Forteza:2021wfq ; Oshita:2022pkc ) showing that the
first several overtones must be taken into account in order to model the
ringdown phase, obtained by accurate numerical relativity simulations in the
beginning of the quasinormal ringing and not only at the last stage. This
indicates also that the actual quasinormal ringing starts much earlier than
expected. Once the overtones are taken into account, the modeling and linear
profile of the ringdown are in the full concordance and allow us to extract
the angular mass and momentum of the black hole. Even though it seems that the
current LIGO/VIRGO observational data do not allow reliable detection of the
overtones in the gravitational-wave signals Capano:2021etf ; Cotesta:2022pci ,
there are indications that in the future LISA experiments the overtones could
be significantly excited and cannot be neglected at the early ringdown phase
Oshita:2022yry .

The paper is organized as follows: Section II is devoted to the basic equation
for the metric and the wave euqation. In Sec. III we briefly describe the
boundary condition for quasinormal modes, explain the method used for their
calculations and review the obtained results. In the last section we discuss a
couple of open questions.

##  II Black hole metric and wave equation

Figure 1: The effective potential for scalar field perturbations $\ell=0$:
$p=1.14$ (top, red), $p=1$ (green), $p=0.876$ (bottom, blue).

The general spherically symmetric metric has the form

$ds^{2}=-A(r,p)dt^{2}+\dfrac{dr^{2}}{B(r,p)}+r^{2}(\sin^{2}\theta
d\phi^{2}+d\theta^{2}),$ (2)

where the dimensionless parameter $p$ is

$p=\dfrac{r_{0}}{\sqrt{2\alpha}},$ (3)

and $r_{0}$ is the black hole radius.

The metric functions $A(r,p)$ and $B(r,p)$ were obtained numerically in Refs.
Lu:2015cqa ; Lu:2015psa . The analytical approximation for the black hole
metric functions $A(r,p)$ and $B(r,p)$ was found in Ref. Kokkotas:2017zwt
and, as these metric functions have rather lengthy form we do not write them
here explicitly and refer the reader to Ref. Kokkotas:2017zwt  instead.

The Schwarzschild solution is the exact solution of the Einstein-Weyl
equations for all $p$, but the second solution also describes an
asymptotically flat black hole existing at some interval $(p_{min},p_{max})$
Lu:2015cqa ; Lu:2015psa :

$\displaystyle p_{min}\approx\dfrac{1054}{1203}\approx 0.876,$ (4a)
$\displaystyle p_{max}\approx 1.14.$ (4b)

The general covariant equation for a minimally coupled scalar field has the
form:

$\dfrac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\Phi)=0,$
(5)

while for an electromagnetic field the general covariant Maxwell equations
read as follows:

$\dfrac{1}{\sqrt{-g}}\partial_{\mu}(F_{\rho\sigma}g^{\rho\nu}g^{\sigma\mu}\sqrt{-g})=0.$
(6)

Here, $A_{\mu}$ is a vector potential and
$F_{\rho\sigma}=\partial_{\rho}A_{\sigma}-\partial_{\sigma}A_{\rho}$. The
function $\Phi$ can be represented in terms of spherical harmonics and radial
part in the following way+:

$\Phi(t,r,y,\phi)=e^{\pm i\omega t}R_{\omega\ell}(r)Y_{\ell(\theta,\phi)}.$
(7)

Here, $Y_{l}(\theta,\phi)$ are the standard spherical harmonics, and $\ell$ is
the multipole number. Then, we rewrite the radial component
$R_{\omega\ell}(r)$ via the new function,

$\Psi_{\omega l}(r)=rR_{\omega l}(r).$ (8)

Thus, taking into account the replacement (7), Eqs (5) and (6) can be
represented in the following general form:

$\dfrac{d^{2}\Psi}{dr_{*}^{2}}+(\omega^{2}-V(r,p))\Psi=0,$ (9)

where the “tortoise coordinate” $r_{*}$ is

$dr_{*}=\dfrac{dr}{\sqrt{A(r,p)B(r,p)}}.$ (10)

The corresponding effective potentials for the scalar $V_{s}$ and
electromagnetic $V_{e}$ perturbations can be written in the following forms
Konoplya:2006rv :

$V_{s}(r,p)=\dfrac{\ell(\ell+1)A(r,p)}{r^{2}}+\ $
$\dfrac{B(r,p)\dfrac{\partial A(r,p)}{\partial r}+A(r,p)\dfrac{\partial
B(r,p)}{\partial r}}{2r},$ (11)
$V_{e}(r,p)=\dfrac{\ell(\ell+1)A(r,p)}{r^{2}}.$ (12)

The effective potentials have the form of a positive definite potential
barrier which decays at the event horizon and at infinity (see an example for
the $\ell=0$ scalar field perturbations in Fig. 1).

Figure 2: Real and imaginary parts of quasinormal frequencies for $\ell=0$
scalar perturbations as a function of $p$ found by the Frobenius method;
$r_{+}=1$.

Figure 3: Real and imaginary parts of quasinormal frequencies for $\ell=1$
electromagnetic perturbations as a function of $p$ found by the Frobenius
method; $r_{+}=1$

##  III Quasinormal modes: sensitivity of overtones

Quasinormal modes satisfy the following boundary conditions

$\Phi(r_{*})\propto e^{\pm i\omega r_{*}},\qquad r_{*}\to\pm\infty,$ (13)

corresponding to purely ingoing waves at the event horizon and purely outgoing
waves at the cosmological horizon. Thus, no waves are coming from either of
the boundaries, which means that a response in perturbations is detected once
the source of perturbations stops acting, or, in other words, proper
oscillation frequencies are observed.

In the previous literature quasinormal modes were studied with the help of
either the third-order WKB method Cai:2015fia  or the sixth-order WKB method
with the Prony approximants and time-domain integration Zinhailo:2019rwd ;
Zinhailo:2018ska . While the time-domain integration and the sixth-order WKB
method usually provides sufficiently good accuracy, there are some
deficiencies: the WKB series converges only asymptotically and cannot
_guarantee_ accuracy, while time-domain integration, especially for $\ell=0$
scalar perturbations allows us to extract values of the frequencies with only
a few digits. In this respect, the Frobenius method which was applied to
finding quasinormal modes by E. Leaver Leaver:1985ax  is based on the
_convergent_ procedure. Therefore, it guarantees results with as many digits
as one wishes. More importantly, while both WKB and time-domain integration
fail when searching for overtones with $n>\ell$, the Leaver method allows one
to find higher overtones accurately as well. That is why, in order to find the
accurate values of QNMs we will apply the Leaver method Leaver:1985ax .

###  III.1 Frobenius method

The wavelike equation (9) has regular singularities at the event horizon
$r=r_{+}$ and an irregular singularity at spatial infinity $r=\infty$.
Further, we will use the new function,

$\Psi(r)=e^{i\omega
r}r^{\lambda}\left(\frac{r-r_{+}}{r}\right)^{-i\omega/f^{\prime}(r_{+})}y(r),$
(14)

where $\lambda$ is defined in such a way that $y(r)$ is regular at $r=\infty$
once $\Psi(r)$ corresponds to the purely outgoing wave at spatial infinity. We
notice that $\Psi(r)$ behaves as a purely ingoing wave at the event horizon,
if $y(r)$ is regular at $r=r_{+}$. Therefore, we represent $y(r)$ in terms of
the following Frobenius series:

$y(r)=\sum_{k=0}^{\infty}a_{k}\left(\frac{r-r_{+}}{r}\right)^{k}.$ (15)

Then, we find that the coefficients $a_{k}$ satisfy some recurrence relation,
which can be reduced to the three-term recurrence relation via Gaussian
elimination (see Ref. Konoplya:2011qq  for details). Finally, using the
recurrence relation coefficients, we find the equation with an infinite
continued fraction with respect to $\omega$, which is satisfied once the
series in Eq. (15) converges at $r=\infty$ i.e., when $\Psi(r)$ obeys the
quasinormal boundary conditions. In order to calculate the infinite continued
fraction, we also use the so called Nollert improvement Nollert:1993zz , which
was generalized in Zhidenko:2006rs  for an arbitrary number of terms of the
recurrence relation.

###  III.2 WKB method

In the frequency domain we will use the semianalytic WKB approach applied by
Will and Schutz Schutz:1985zz  for finding quasinormal modes. The Will-Schutz
formula was extended to higher orders in Refs. Iyer:1986np ; Konoplya:2003ii ;
Matyjasek:2017psv  and made even more accurate when using the Padé
approximants in Refs. Matyjasek:2017psv ; Hatsuda:2019eoj . The general WKB
formula has the form Konoplya:2019hlu ,

$\displaystyle\omega^{2}$ $\displaystyle=$ $\displaystyle V_{0}+A_{2}({\cal
K}^{2})+A_{4}({\cal K}^{2})+A_{6}({\cal K}^{2})+\ldots$ $\displaystyle-$
$\displaystyle i{\cal K}\sqrt{-2V_{2}}\left(1+A_{3}({\cal K}^{2})+A_{5}({\cal
K}^{2})+A_{7}({\cal K}^{2})+\ldots\right),$

where ${\cal K}=n+1/2$ is a half-integer. The corrections $A_{k}({\cal
K}^{2})$ of the order $k$ to the eikonal formula are polynomials of ${\cal
K}^{2}$ with rational coefficients and depend on the values of higher
derivatives of the potential $V(r)$ in its maximum. In order to increase the
accuracy of the WKB formula, we will follow the procedure of Matyjasek and
Opala Matyjasek:2017psv  and use the Padé approximants. Here, for the
fundamental mode $n=0$, we will use the 13th-order WKB method with
$\tilde{m}=7$, where $\tilde{m}$ is defined in Refs. Matyjasek:2017psv ;
Konoplya:2019hlu , because this choice provides the best accuracy in the
Schwarzschild limit and there is hope that this will be the case for more
general metrics. We will see that the 13th-order WKB formula gives much better
agreement with the accurate Frobenius method than the sixth and lower orders
used previously in the literature for this case.

###  III.3 Numerical and semianalytical calculations of quasinormal modes

As was mentioned earlier, the mass of the black holes diminishes when $p$ is
increased, achieving zero at some extremal value of $p$. Therefore, such light
near extremal black holes apparently could not be created as a result of the
collapse of stars. Therefore, first of all, the astrophysically motivated
range of the coupling parameter $p$, that is relatively near the Schwarzschild
limit, will be considered in this section.

From the accurate Frobenius data shown in Figs. 2 and 3 for a fixed value of
the radius of the event horizon $r_{+}=1$ we can see that for both scalar and
electromagnetic perturbations at the lowest multipole number $\ell$ the real
oscillation frequency of the fundamental mode $n=0$ is increased only by a few
percent, when the coupling $p$ runs from $0.876$ until $0.92$. The damping
rate is changed even less: by a fraction of one percent. Thus, the fundamental
mode is almost indistinguishable from the Schwarzschild one in this range of
$p$.

Though, as the mass changes with $p$, the black hole might be distinguishable
from the Schwarzschild solution if one takes into account the inspiral phase.
At the same time, we can see that higher overtones are much more sensitive to
variation of the coupling $p$. Thus, the $n=6$ mode changes already by a few
hundred percent, which is 2 orders higher a rate of change than that of the
fundamental mode. This outburst of overtones takes place mainly for the real
oscillation frequency, while the damping rate is increased at an almost
equidistant spacing for large $n$, so that it changes roughly proportionally
to its value in the Schwarzschild limit.

From Figs. 4 and 5 we can see that near the Schwarzschild limit, the WKB
method leads relatively close results at the 7th and 13th orders and is close
to the sixth order results obtained in Ref. Zinhailo:2018ska . With the help
of the Frobenius method in the limit $p\approx 0.876$ we reproduce the
Schwarzschild values of quasinormal modes with great accuracy for all
overtones under consideration, while for the fundamental mode we reproduce the
results of Ref. Zinhailo:2018ska  within the expected accuracy of the WKB
method Schutz:1985zz ; Iyer:1986np ; Konoplya:2003ii ; Matyjasek:2017psv ;
Konoplya:2019hlu  at the sixth order used there.

Summarizing the accuracy provided by different methods and in order to
distinguish the inaccuracy due to the method from that coming from the
parametrization of the metric we conclude that while the Frobenius method is
accurate for any values of $p$ and overtone number $n$, the WKB method can
only be safely used for calculation of the fundamental mode and, and the 13th
order is necessary to achieve sufficient accuracy for large values of $p$.

##  IV Higher-order parametrization for the metric

Figure 4: Real and imaginary parts of the fundamental quasinormal mode
$\ell=n=0$ for scalar perturbations as a function of $p$ calculated with the
help of the 13th WKB order (lower curve for $Re\omega$ and upper curve for
$Im\omega$) and the Frobenius method; $\tilde{m}=7$, $r_{+}=1$. The
nonmonotonic behavior near $p\approx 1.06$ is due to an insufficient number of
orders (four) of the parametrization for the metric function.

Figure 5: Real and imaginary parts of the fundamental quasinormal mode $\ell=1$, $n=0$ for electromagnetic perturbations as a function of $p$ calculated with the help of the 13th WKB order with $\tilde{m}=7$, and with the Frobenius method and 7th WKB order ($\tilde{m}=7$); $r_{+}=1$. The results of the 13th WKB order and accurate Frobenius data are very close (upper curve for $Re\omega$ and lower one for $Im\omega$) and cannot be distinguished on the plot. The nonmonotonic behavior near $p\approx 1.06$ is due to an insufficient number of orders (four) of the parametrization for the metric function. $p$ | 4th order | 5th order | 6th order | 7th order
---|---|---|---|---
1.060 | $0.272020-0.383228i$ | $0.275324-0.382319i$ | $0.276697-0.380141i$ | $0.276618-0.380274i$
1.061 | $0.272979-0.384515i$ | $0.275658-0.383442i$ | $0.277033-0.381243i$ | $0.276953-0.381377i$
1.062 | $0.276117-0.385789i$ | $0.275995-0.384565i$ | $0.277371-0.382345i$ | $0.277292-0.382479i$
1.063 | $0.266280-0.381901i$ | $0.276324-0.385698i$ | $0.277701-0.383457i$ | $0.277625-0.383592i$
1.064 | $0.270399-0.386619i$ | $0.276664-0.386822i$ | $0.278042-0.384560i$ | $0.277963-0.384696i$
1.065 | $0.271456-0.388172i$ | $0.277000-0.387952i$ | $0.278042-0.384560i$ | $0.278302-0.385806i$
1.066 | $0.272082-0.389472i$ | $0.277338-0.389083i$ | $0.278718-0.386779i$ | $0.278641-0.386916i$
Table 1: The fundamental ($\ell=n=0$) mode for the scalar field calculated by the Frobenius method for the metrics which are truncation of the parametrization at different orders; $r_{0}=1$. We can see that the quasinormal modes converge with the increasing of the order and the difference between the 6th and 7th order results is a small fraction of $1\%$. $p$ | 4th order | 5th order | 6th order | 7th order
---|---|---|---|---
1.060 | $0.592591-0.300992i$ | $0.592580-0.301161i$ | $0.592699-0.301224i$ | $0.592689-0.301232i$
1.061 | $0.593138-0.301681i$ | $0.593014-0.301842i$ | $0.593134-0.301906i$ | $0.593124-0.301915i$
1.062 | $0.594020-0.302348i$ | $0.593448-0.302524i$ | $0.593568-0.302588i$ | $0.593559-0.302597i$
1.063 | $0.592245-0.302414i$ | $0.593881-0.303208i$ | $0.594001-0.303273i$ | $0.593992-0.303281i$
1.064 | $0.593756-0.303607i$ | $0.594312-0.303889i$ | $0.594434-0.303955i$ | $0.594424-0.303964i$
1.065 | $0.594346-0.304331i$ | $0.594743-0.304573i$ | $0.594865-0.304639i$ | $0.594855-0.304648i$
1.066 | $0.594839-0.305028i$ | $0.5951726-0.305256i$ | $0.595295-0.305323i$ | $0.594993-0.305655i$
Table 2: The fundamental ($\ell=1$, $n=0$) mode for the electromagnetic field calculated by the Frobenius method for the metrics which are truncation of the parametrization at different orders; $r_{0}=1$. We can see that the quasinormal modes converge with the increasing of the order and the difference between the 6th and 7th order results is a small fraction of one percent. $p$ | $a_{5}$ | $b_{5}$ | $a_{6}$ | $b_{6}$ | $a_{7}$ | $b_{7}$
---|---|---|---|---|---|---
1.060 | $-1/3743$ | $-603/1802$ | $-301/2336$ | $-563/1673$ | $1413/1060$ | $11/636$
1.061 | $-1/5867$ | $-564/1685$ | $-533/4174$ | $-1565/4649$ | $371/274$ | $19/1099$
1.062 | $-1/13431$ | $-984/2939$ | $-165/1304$ | $-461/1369$ | $337/245$ | $7/405$
1.063 | $1/48062$ | $-571/1705$ | $-85/678$ | $-351/1042$ | $4/3$ | $1/57$
1.064 | $1/8669$ | $-878/2621$ | $-155/1248$ | $-650/1929$ | $1535/1081$ | $12/695$
1.065 | $1/4780$ | $-641/1913$ | $-475/3861$ | $-1112/3299$ | $1013/702$ | $19/1101$
1.066 | $1/3308$ | $-5491/16383$ | $-343/2815$ | $-701/2079$ | $5447/3714$ | $39/2261$
Table 3: Higher coefficients of the metric parametrization obtained with the
help of the automatic procedure developed in Kokkotas:2017zwt ; $r_{0}=1$. The
lower coefficients are found as functions of $p$ in Kokkotas:2017zwt .

Following the parametrization procedure of Ref. Rezzolla:2014mua  the new
functions $A$ and $B$ are defined as follows:

$\displaystyle A(r)$ $\displaystyle\equiv$ $\displaystyle xh(x)\,,$ (17)
$\displaystyle\frac{A(r)}{B(r)}$ $\displaystyle\equiv$ $\displaystyle
f(x)^{2}\,,$ (18)

where $x$ is the dimensionless compact coordinate

$x\equiv 1-\frac{r_{0}}{r}\,.$ (19)

We write the above two functions in the following form:

$\displaystyle h(x)$ $\displaystyle=$ $\displaystyle
1-\epsilon(1-x)+(a_{0}-\epsilon)(1-x)^{2}+{\tilde{h}}(x)(1-x)^{3},$
$\displaystyle f(x)$ $\displaystyle=$ $\displaystyle
1+b_{0}(1-x)+{\tilde{f}}(x)(1-x)^{2},$ (20)

where ${\tilde{h}}(x)$ and ${\tilde{f}}(x)$ are expressed in terms of the
continued fractions, in order to describe the metric near the event horizon
$x=0$:

$\displaystyle{\tilde{h}}(x)=\frac{a_{1}}{\displaystyle 1+\frac{\displaystyle
a_{2}x}{\displaystyle 1+\frac{\displaystyle a_{3}x}{\displaystyle
1+\frac{\displaystyle a_{4}x}{\displaystyle 1+\ldots}}}}\,,$
$\displaystyle{\tilde{f}}(x)=\frac{b_{1}}{\displaystyle 1+\frac{\displaystyle
b_{2}x}{\displaystyle 1+\frac{\displaystyle b_{3}x}{\displaystyle
1+\frac{\displaystyle b_{4}x}{\displaystyle 1+\ldots}}}}\,.$ (21)

At the event horizon, one has
${\tilde{h}}(0)=a_{1},\quad{\tilde{f}}(0)=b_{1}$.

From Figs. 4 and 5 we see that at a fixed radius of the event horizon, the
real oscillation frequency grows monotonically and damping rate decreases once
$p$ is growing, except for a small region near $p\approx.1.06$. There, a kind
of “splash ” occurs, which is not usual for quasinormal spectra of black
holes, because normally the frequencies depend on physical parameters in a
smooth way.

Quasinormal modes obtained here in Figs. 2 and 3 (for values $p$ near the
Schwarzschild limit) as well as in Zinhailo:2018ska ; Zinhailo:2018ska  were
found for the fourth-order approximation Kokkotas:2017zwt  of the numerical
solution Lu:2015cqa . This means that $a_{5}=a_{6}...=0$, and
$b_{5}=b_{6}...=0$. In the regime of relatively small deviation from the
Schwarzschild geometry (as in figs. 2 and 3 here) the fourth-order
approximation is sufficient that is, the quasinormal modes do not change
significantly if the expansion is truncated at higher orders. In other words,
the convergence regime is achieved.

However, as we can see from Tables I and II, at higher $p$ in the range where
the nonmonotonic dependence of $\omega$ on $p$ occurs, this is not so, and the
difference between quasinormal modes calculated for the forth- and fifth-order
metrics may reach a few percent. At the same time, when higher-order
approximations are used according to the general automatic procedure in Ref.
Kokkotas:2017zwt , the quasinormal modes converge quickly. For example, the
difference between quasinormal modes for the 5-th and 6-th order metrics is
about one percent or less, while the difference between the sixth- and
seventh-order is already a small fraction of $1\%$. We can also see that the
convergence selectively depends on $p$ in this range of parameters, being the
worst for $p=1.065$, but better for smaller and larger $p$. The coefficients
of the metric parametrization $\epsilon$, $a_{i}$ and $b_{i}$ for $i\leq 4$
are found in general form, that is, as functions of $p$ and $r_{0}$ in Ref.
Kokkotas:2017zwt , while higher-order coefficients can be seen in Table III
here for particular values of $p$ considered in Tables I and II. Using the
above data in Kokkotas:2017zwt  and Table III one can obtain quasinormal modes
for other spin fields and multipole numbers $\ell$.

##  V Conclusions

The present paper is revising the earlier studies of quasinormal modes of non-
Schwarzschild asymptotically flat black holes in the Einstein-Weyl gravity in
two aspects: accurate calculation of overtones which cannot be found with the
WKB or time-domain integration methods and testing the convergence of the
parametrization for the metric. We have found that:

  * •

Here we have accurately calculated quasinormal modes of the Einstein-Weyl
black hole by two alternative methods which achieve the regime of convergence
and are very well in concordance with each other. This way we have shown that
the discrepancy between Refs. Zinhailo:2019rwd ; Zinhailo:2018ska  and Ref.
Cai:2015fia  is resolved in the favor of Ref. Zinhailo:2019rwd ;
Zinhailo:2018ska  at small coupling, while at larger coupling the results of
Zinhailo:2019rwd ; Zinhailo:2018ska  are corrected. In other words, the
accurate quasinormal modes obtained here are in a very good agreement with
Refs. Zinhailo:2019rwd ; Zinhailo:2018ska  in the regime of small coupling and
significantly differ from Cai:2015fia .

  * •

The first several overtones are very sensitive to the change of the coupling
constant. While the fundamental mode is changed only within a couple of
percent, the higher overtones can change by hundreds of percent. This may be
important for reproduction of the quasinormal ringing at earlier times
Giesler:2019uxc , for the so-called mode instability Jaramillo:2020tuu  and
for testing the geometry of the event horizon Konoplya:2022pbc . In Ref.
Jaramillo:2021tmt  it is shown that the overtone outburst may already occur
within the near-future detection range while in Ref. Oshita:2022yry  the
possibility of the detection of overtones with LISA is argued.

  * •

The unusual nonmonotonic dependence of the quasinormal frequencies upon the
coupling parameter $p$ which was found in earlier publications
Zinhailo:2019rwd ; Zinhailo:2018ska  is the result of the insufficient number
of terms at which the metric parametrization was truncated, when approximating
the numerical solution by an analytical expression. Extension to higher orders
allows one to achieve convergence when computing quasinormal modes.

Our work could be extended in a number of ways, first of all, by consideration
of the Dirac and gravitational perturbations. The latter would be a highly
challenging problem, because the linearization of perturbation equations in
the higher-derivative gravity will most probably lead to a system of chained
differential equations of higher than the second order. Nevertheless, the
phenomenon of sensitivity of the overtones found here is connected with the
near-horizon deformation of the geometry of the non-Schwarzschild black hole
Konoplya:2022pbc ; Konoplya:2022hll . Therefore, qualitatively the same
behavior of overtones must take place for perturbations of other spin fields
as well.

###### Acknowledgements.

The author would like to thank the support of the Czech Science Foundation
(GAČR) under Grant No. 19-03950S and A. Zhidenko for most useful discussions.

## References

  * (1) K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999) [gr-qc/9909058].
  * (2) E. Berti, V. Cardoso and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009) [arXiv:0905.2975 [gr-qc]].
  * (3) R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011) [arXiv:1102.4014 [gr-qc]].
  * (4) B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 6, 061102 (2016) [arXiv:1602.03837 [gr-qc]].
  * (5) R. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. X 11, 021053 (2021) [arXiv:2010.14527 [gr-qc]].
  * (6) B. P. Abbott et al. [KAGRA, LIGO Scientific, Virgo and VIRGO], Living Rev. Rel. 21, no.1, 3 (2018) [arXiv:1304.0670 [gr-qc]].
  * (7) R. Konoplya and A. Zhidenko, Phys. Lett. B 756, 350 (2016) [arXiv:1602.04738 [gr-qc]].
  * (8) N. Yunes, K. Yagi and F. Pretorius, Phys. Rev. D 94, no. 8, 084002 (2016) [arXiv:1603.08955 [gr-qc]].
  * (9) E. Berti, K. Yagi, H. Yang and N. Yunes, Gen. Rel. Grav. 50, no. 5, 49 (2018) [arXiv:1801.03587 [gr-qc]].
  * (10) H. Lu, A. Perkins, C. N. Pope and K. S. Stelle, Phys. Rev. Lett. 114, no.17, 171601 (2015) [arXiv:1502.01028 [hep-th]].
  * (11) L. Rezzolla and A. Zhidenko, Phys. Rev. D 90, no.8, 084009 (2014) [arXiv:1407.3086 [gr-qc]].
  * (12) R. Konoplya, L. Rezzolla and A. Zhidenko, Phys. Rev. D 93, no.6, 064015 (2016) [arXiv:1602.02378 [gr-qc]].
  * (13) Z. Younsi, A. Zhidenko, L. Rezzolla, R. Konoplya and Y. Mizuno, Phys. Rev. D 94, no.8, 084025 (2016) [arXiv:1607.05767 [gr-qc]].
  * (14) K. Kokkotas, R. A. Konoplya and A. Zhidenko, Phys. Rev. D 96, no.6, 064007 (2017) [arXiv:1705.09875 [gr-qc]].
  * (15) J. Podolsky, R. Svarc, V. Pravda and A. Pravdova, Phys. Rev. D 98, no. 2, 021502 (2018) [arXiv:1806.08209 [gr-qc]]; H. Lü, A. Perkins, C. N. Pope and K. S. Stelle, Phys. Rev. D 96, no. 4, 046006 (2017) [arXiv:1704.05493 [hep-th]]; B. Holdom and J. Ren, Phys. Rev. D 95, no. 8, 084034 (2017) [arXiv:1612.04889 [gr-qc]].
  * (16) K. Lin, W. L. Qian, A. B. Pavan and E. Abdalla, EPL 114, no. 6, 60006 (2016) [arXiv:1607.04473 [gr-qc]].
  * (17) A. F. Zinhailo, Eur. Phys. J. C 79, no.11, 912 (2019) doi:10.1140/epjc/s10052-019-7425-9 [arXiv:1909.12664 [gr-qc]].
  * (18) A. F. Zinhailo, Eur. Phys. J. C 78, no.12, 992 (2018) [arXiv:1809.03913 [gr-qc]].
  * (19) R. A. Konoplya and A. Zhidenko, [arXiv:2209.00679 [gr-qc]].
  * (20) R. A. Konoplya, A. F. Zinhailo, J. Kunz, Z. Stuchlik and A. Zhidenko, JCAP 10 (2022) 091 [arXiv:2206.14714 [gr-qc]].
  * (21) M. Giesler, M. Isi, M. A. Scheel and S. Teukolsky, Phys. Rev. X 9, no.4, 041060 (2019) [arXiv:1903.08284 [gr-qc]].
  * (22) N. Oshita, Phys. Rev. D 104, no.12, 12 (2021) [arXiv:2109.09757 [gr-qc]].
  * (23) X. J. Forteza and P. Mourier, Phys. Rev. D 104, no.12, 124072 (2021) [arXiv:2107.11829 [gr-qc]].
  * (24) N. Oshita, [arXiv:2208.02923 [gr-qc]].
  * (25) C. D. Capano, M. Cabero, J. Westerweck, J. Abedi, S. Kastha, A. H. Nitz, Y. F. Wang, A. B. Nielsen and B. Krishnan, [arXiv:2105.05238 [gr-qc]].
  * (26) R. Cotesta, G. Carullo, E. Berti and V. Cardoso, Phys. Rev. Lett. 129, no.11, 111102 (2022) [arXiv:2201.00822 [gr-qc]].
  * (27) N. Oshita and D. Tsuna, [arXiv:2210.14049 [gr-qc]].
  * (28) J. L. Jaramillo, R. Panosso Macedo and L. Al Sheikh, Phys. Rev. X 11, no.3, 031003 (2021) [arXiv:2004.06434 [gr-qc]].
  * (29) E. W. Leaver, Proc. Roy. Soc. Lond. A 402, 285-298 (1985)
  * (30) H. P. Nollert, Phys. Rev. D 47, 5253-5258 (1993)
  * (31) A. Zhidenko, Phys. Rev. D 74, 064017 (2006) [arXiv:gr-qc/0607133 [gr-qc]].
  * (32) R. A. Konoplya and A. F. Zinhailo, Phys. Rev. D 99, no.10, 104060 (2019) [arXiv:1904.05341 [gr-qc]].
  * (33) C. Wu, D. C. Zou and M. Zhang, Nucl. Phys. B 952, 114942 (2020) [arXiv:1904.10193 [gr-qc]].
  * (34) Y. F. Cai, G. Cheng, J. Liu, M. Wang and H. Zhang, JHEP 1601, 108 (2016) [arXiv:1508.04776 [hep-th]].
  * (35) H. Lu, A. Perkins, C. N. Pope and K. S. Stelle, Phys. Rev. D 92, no. 12, 124019 (2015) [arXiv:1508.00010 [hep-th]].
  * (36) R. A. Konoplya and A. Zhidenko, Phys. Lett. B 644, 186-191 (2007) [arXiv:gr-qc/0605082 [gr-qc]].
  * (37) B. F. Schutz and C. M. Will, Astrophys. J. 291, L33 (1985).
  * (38) S. Iyer and C. M. Will, Phys. Rev. D 35, 3621 (1987).
  * (39) R. A. Konoplya, Phys. Rev. D 68, 024018 (2003) [gr-qc/0303052].
  * (40) J. Matyjasek and M. Opala, Phys. Rev. D 96, no. 2, 024011 (2017) [arXiv:1704.00361 [gr-qc]].
  * (41) Y. Hatsuda, Phys. Rev. D 101, no.2, 024008 (2020) [arXiv:1906.07232 [gr-qc]].
  * (42) R. A. Konoplya, A. Zhidenko and A. F. Zinhailo, Class. Quant. Grav. 36, 155002 (2019) [arXiv:1904.10333 [gr-qc]].
  * (43) J. L. Jaramillo, R. Panosso Macedo and L. A. Sheikh, Phys. Rev. Lett. 128, no.21, 211102 (2022) doi:10.1103/PhysRevLett.128.211102 [arXiv:2105.03451 [gr-qc]].
  * (44) N. Oshita and D. Tsuna, [arXiv:2210.14049 [gr-qc]].