July 2022 IPMU22-0022
revised version

Formation of primordial black holes after Starobinsky inflation

Daniel Frolovsky a, Sergei V. Ketov a,b,c,d,# and Sultan Saburov a

a Interdisciplinary Research Laboratory, Tomsk State University

36 Lenin Avenue, Tomsk 634050, Russia

b Department of Physics, Tokyo Metropolitan University

1-1 Minami-ohsawa, Hachioji-shi, Tokyo 192-0397, Japan

c Research School of High-Energy Physics, Tomsk Polytechnic University

2a Lenin Avenue, Tomsk 634028, Russia

d Kavli Institute for the Physics and Mathematics of the Universe (WPI)

The University of Tokyo Institutes for Advanced Study,

Kashiwa 277-8583, Japan

#<EMAIL_ADDRESS>
Abstract

We adapted the Appleby-Battye-Starobinsky model of $F(R)$ gravity towards
describing double cosmological inflation and formation of primordial black
holes with masses up to $10^{19}$ g in the single-field model. We found that
it is possible to get an enhancement of the power spectrum of scalar curvature
perturbations to the level beyond the Hawking (black hole evaporation) limit
of $10^{15}$ g, so that the primordial black holes resulting from
gravitational collapse of those large primordial perturbations can survive in
the present universe and form part of cold dark matter. Our results agree with
the current measurements of cosmic microwave background radiation within
$3\sigma$ but require fine-tuning of the parameters.

##  1 Introduction

One of the current candidates for cold dark matter (CDM) is given by
primordial black holes (PBH) presumably formed in the early Universe during or
after inflation [1, 2, 3, 4, 5, 6, 7, 8].

Cosmological inflation in the early Universe can be described by the
Starobinsky model [9], see also Ref. [10] for a recent review. The Starobinsky
model perfectly fits current observations of the cosmic microwave background
(CMB) radiation, relies on gravitational interactions, while its only
(inflaton mass) parameter $M$ is fixed by the CMB amplitude. Starobinsky’s
inflaton (scalaron) has the clear origin as the physical degree of freedom of
the higher-derivative gravity, and can be identified with the Nambu-Goldstone
boson related to spontaneous breaking of the scale invariance [11, 12]. The
current observational constraints on inflation are available in Refs. [13, 14,
15].

The Starobinsky model is based on the non-perturbative $(R+R^{2}/6M^{2})$
gravity and does not lead to PBH production, so the question arises whether it
can be modified in the context of more general $F(R)$ gravity models in order
to generate PBH. The principal answer to the same question in the context of
the dual scalar-tensor (or quintessence) gravity models is known to be in
affirmative by demanding the existence of a near-inflection point in the
inflaton scalar potential by using the double inflation scenario leading to an
enhancement of the power spectrum of scalar perturbations [16, 17, 18]. In the
literature, the single-field models of PBH production are unrelated to the
Starobinsky inflation and lead to a reduction of the effective number of
e-folds (relevant to CMB) and the related decrease in the value of the tilt
$n_{s}$ of CMB scalar perturbations, see e.g., Refs. [16, 19, 20]. For
instance, when insisting on the best CMB fit (i.e. the present central
observational value) and robust predictions, one gets the typical PBH masses
of the order $10^{8}$ g [21] that is well below the Hawking black hole
evaporation limit of $10^{15}$ g, so that those light PBH do not survive in
the present universe and cannot be part of CDM.

In this letter we revisit those issues in the context of Starobinsky’s
inflation by relaxing the agreement with the CMB measurements up to $3\sigma$,
while simultaneously maximizing the total number of e-folds for double
inflation. The double inflation with a near inflection point is achieved by a
careful choice of the $F(R)$-gravity function that asymptotically leads to
Starobinsky’s $R^{2}$ term in the large (scalar) curvature limit.

Our paper is organized as follows. In Sec. 2 we introduce our $F(R)$ gravity
model and the corresponding scalar potential. Inflationary dynamics is
investigated in Sec. 3. In Sec. 4 we study large scalar perturbations and
demonstrate the existence of a significant peak in the power spectrum leading
to PBH formation. Sec. 5 is our Conclusion.

##  2 The model

A generic $F(R)$ gravity model leads to singularities in cosmological
evolution and does not agree with observations. A viable cosmological model
has to obey the no-ghost or stability conditions, must have a well-defined
Newtonian limit, satisfy the Solar system tests, and keep successes of
inflation and the standard cosmology. All that makes a choice of new viable
$F(R)$ gravity functions non-trivial.

Those consistency issues were addressed by Appleby, Battye and Starobinsky in
Ref. [22] where they proposed a viable $F(R)$ gravity model with the action

$S={{M_{\rm Pl}^{2}}\over{2}}\int d^{4}x\sqrt{-g}F(R)~{}~{},$ (1)

simultaneously describing the Starobinsky inflation and the present dark
energy, by using the following $F$-function:

$F(R)=\epsilon{\raise-3.0pt\hbox{${\hskip
0.75pt}\\!_{AB}$}}g\ln\left[{{\cosh\left({\textstyle{R\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{\epsilon_{AB}}$}}}}-b\right)}\over{\cosh(b)}}\right]+(1-g)R+{{R^{2}}\over{6M^{2}}}~{}.$
(2)

Here $M_{\rm Pl}\sim 10^{18}$ GeV denotes the (reduced) Planck mass, and the
$M$ stands for the inflaton (Starobinsky’s scalaron) mass $M\sim 10^{-5}M_{\rm
Pl}\sim 10^{13}$ GeV. The $\epsilon{\raise-3.0pt\hbox{${\hskip
0.75pt}\\!_{AB}$}}$ is defined by

$\epsilon_{\raise-3.0pt\hbox{${\hskip
0.75pt}\\!_{AB}$}}={{R_{0}}\over{2g\ln(1+e^{2b})}}~{}~{},$ (3)

so that the model has three free parameters: the dimensional parameter
$R_{0}$, and the two dimensionless parameters $g$ and $b$. In Ref. [22], the
parameter $R_{0}$ represented the vacuum value of the scalar curvature with
$\sqrt{R_{0}}\sim 10^{-33}$ eV, in order to describe the present dark energy,
whereas the remaining parameters $(g,b)$ determined the shape of the inflaton
scalar potential. The allowed range of the parameters $(b,g)$ was given in
Fig. 4 of Ref. [22], which specifies (above the curve) the parameter values
needed to get a metastable de Sitter vacuum.

We use the same function (2) with the important difference in the physical
scale and physical interpretation: instead of describing the dark energy we
want to describe double inflation for the purpose of PBH production, by
assigning the much higher value of $R_{0}$ just under the inflationary scale,
$\sqrt{R_{0}}<H_{\rm inf.}\sim 10^{14}$ GeV. The power spectrum is well
constrained on the CMB scale $H_{\rm inf.}$ by observational data, but it is
not on smaller scales.

The standard conversion from the $F(R)$ gravity description to the (dual)
scalar-tensor gravity description in terms of the canonical inflaton scalar
field $\phi$ with the potential $V$ is given by 111See e.g., Ref. [23] for
more details.

$\phi=\sqrt{{\textstyle{3\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}}M_{\rm Pl}^{2}\ln\left({{2}\over{M_{\rm
Pl}^{2}}}F^{\prime}\right)~{},\quad V(R)=M_{\rm
Pl}^{4}{{F^{\prime}R-F}\over{4(F^{\prime})^{2}}}~{},$ (4)

where the primes denote the derivatives with respect to the argument (scalar
curvature $R$).

In the case of the $F$-function (2), we find the potential

$V(R)={{y^{2}}\over{2}}\Bigg{[}{{R^{2}}\over{6M^{2}}}+gR\tanh\Bigg{(}{{R}\over{\epsilon_{AB}}}-b\Bigg{)}-g\epsilon{\raise-3.0pt\hbox{${\hskip
0.75pt}\\!_{AB}$}}\ln\Bigg{(}{{\cosh({{R}\over{\epsilon_{AB}}}-b)}\over{\cosh(b)}}\Bigg{)}\Bigg{]}~{}~{},$
(5)

where we have used the notation [23]

$y\equiv\exp\left(-\sqrt{{\textstyle{2\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{3}$}}}}}\phi/M_{\rm
Pl}\right)=\Bigg{[}1-g+{{R}\over{3M^{2}}}+g\tanh\Bigg{(}{{R}\over{\epsilon_{AB}}}-b\Bigg{)}\Bigg{]}^{-1}~{}~{}.$
(6)

The quantity $y$ is small during the Starobinsky inflation, so that it can be
used as an expansion variable. The scalar potential $V(\phi)$ cannot be
obtained analytically because it requires finding the inverse function
$R(\phi)$ from Eq. (6). However, it is possible numerically.

The profile of the canonical inflaton potential for the parameter values
$g=0.41$, $b=2.89$ and $R_{0}=0.1M^{2}$ is given in Fig. 1. There are three
physically different regimes: (I) the large curvature regime (relevant to CMB)
where the potential is the same as that in the Starobinsky model governed by
the last term in Eq. (2), (II) the lower (still high) curvature regime where
the scalar potential has a shallow meta-stable de Sitter minimum and a small
bump on the left of it, (III) the low curvature regime near the Minkowski
vacuum, relevant for reheating after inflation. Actually, the free parameters
$(g,b)$ were chosen to get the desired profile of the scalar potential in Fig.
1. A significant change in the parameter values may destroy this profile. The
regimes (I+II) realize double inflation.

Figure 1:  The shape of the inflaton potential $V(\phi)$ (left). Zooming the
potential $V(\phi)$ for the selected values of the inflaton field $\phi$
(right). The parameters are given by $g=0.41$, $b=2.89$ and $R_{0}=0.1M^{2}$.

##  3 Double inflation

In this Section we use the standard formulae for describing inflation with a
given (canonical) inflaton scalar potential, see e.g., Ref. [24] for a review.

The equations of motion are given by

$\ddot{\phi}+3H\dot{\phi}+V^{\prime}(\phi)=0,\qquad
3H^{2}={{1}\over{M^{2}_{\rm
Pl}}}\left[{\textstyle{1\over\vphantom{2}\smash{\raise
0.60275pt\hbox{$\scriptstyle{2}$}}}}\dot{\phi}^{2}+V(\phi)\right]~{},\qquad\dot{H}=-{{1}\over{2M^{2}_{\rm
Pl}}}\dot{\phi}^{2}~{}~{},$ (7)

where $H(t)$ is Hubble function, and the dots denote the derivatives with
respect to time.

The duration of inflation is measured by the e-foldings number

$N=\int^{t}_{t_{\rm end}}H(t)dt\approx{{1}\over{M_{\rm
Pl}^{2}}}\int_{\phi}^{\phi_{in}}{{V(\phi)}\over{V^{\prime}(\phi)}}d\phi~{}~{}.$
(8)

The slow roll parameters in terms of the scalar potential are given by

$\epsilon(\phi)={{M_{\rm
Pl}^{2}}\over{2}}\left[{{V^{\prime}(\phi)}\over{V(\phi)}}\right]^{2},\qquad\eta(\phi)=M_{\rm
Pl}^{2}\left|{{V^{\prime\prime}(\phi)}\over{V(\phi)}}\right|~{}~{},$ (9)

and they are to be evaluated at the horizon crossing ($\phi_{\rm in}$). We
also use the slow-roll parameters in terms of the Hubble flow functions,

$\epsilon(t)=-{{\dot{H}}\over{H^{2}}}~{}~{},\qquad\eta(t)={{\dot{\epsilon}}\over{H\epsilon}}~{}~{},$
(10)

with the time $t_{c}$ at the horizon crossing. As the running arguments we use
scalar curvature $R$, cosmic time $t$, inflaton field $\phi$ and e-folds $N$
that are all related to each other.

The CMB tilt of scalar perturbations $(n_{s})$ and the CMB tensor-to-scalar
ratio $(r)$ are given by

$n_{s}=1+2\eta(\phi_{\rm in})-6\epsilon(\phi_{\rm
in})=1-2\epsilon(t_{c})-\eta(t_{c})~{}~{},\qquad
r=16\epsilon(\phi_{in})=16\epsilon(t_{c})~{}.$ (11)

Our strategy is to numerically solve the equations of motion by using the
initial conditions adapted to the original Starobinsky model, by choosing
$\phi_{\rm in}\approx 5.17M_{\rm Pl}$ and $\dot{\phi}_{\rm in}=0$. Since the
Starobinsky inflationary solution is an attractor [9], the dependence upon
local changes of the initial conditions should be irrelevant, except for a
duration of inflation. Our numerical solutions to the Hubble function and the
first slow roll parameter are given in Fig. 2.

Figure 2: On the left: the Hubble function $H(t)$. On the right: the first
slow roll parameter $\epsilon(t)$. The initial conditions are $\phi_{\rm
in}=5.17M_{\rm Pl}$ and $\dot{\phi}_{\rm in}=0$. The end of each stage of
inflation is defined by $\epsilon(t_{\rm end})=1$.

As is clear from Fig. 2, there is the double inflation indeed, with the two
plateaus leading to slow-roll of the inflaton. Between the two stages of slow-
roll inflation, there is a short phase of the so-called ”ultra-slow-roll”
where the first slow-roll parameter becomes very small, and the slow-roll
approximation is broken. Actually, the ultra-slow-roll regime implies fast
rolling of the inflaton.

The end of inflation ($\epsilon=1)$ is achieved at $\phi_{\rm end}\approx
0.9M_{\rm Pl}$ that comes earlier than $\eta=1$ at $\phi_{\rm end}\approx
0.64M_{\rm Pl}$. We demand the total e-foldings number, given by a sum of two
slow-roll stages of inflation, $N_{\rm total}=N_{e}+\Delta N$, to be as large
as possible, $N_{\rm total}\approx 60\div 65$.

##  4 Power spectrum of perturbations and PBH masses

The simple analytic formula for the power spectrum of scalar (curvature)
perturbations [16]

$P_{R}={{H^{2}}\over{8M^{2}_{\rm Pl}\pi^{2}\epsilon}}$ (12)

gives a good approximation as long as the Hubble flow slow-roll parameters
(10) are much less than one. The exact power spectrum should be derived by
solving the Mukhanov-Sasaki (MS) equation [19]. As is clear from our plot on
the right side of Fig. 2, Eq. (12) does not apply for a small relevant part of
the power spectrum with $tM\approx 110\pm 5$.

An ultra-slow-roll regime (also seen on the right side of Fig. 2) leads to an
enhancement of the power spectrum, i.e. the appearance of a peak related to
large perturbations. The height of the peak should be at least $10^{6}$ times
higher than the reference CMB value $P_{R}\sim 10^{-9}$, in order to form the
PBH surviving until the present times.

The scalar tilt $n_{s}$ is related to the power spectrum via the relation
$n_{s}={{d\ln P_{R}}\over{d\ln k}}$ in terms of the scale variable
$k=aH=\dot{a}$. We computed the power spectrum numerically, from Eq. (12) by
using Mathematica and the original code for numerical calculations on a
computer. Our results are given by Fig. 3. The height of the peak in Fig. 3
corresponds to $P_{R}(k)\sim 10^{-5}$ that is the three orders of magnitude
less than its actual value $P_{R}(k)\sim 10^{-2}$.

Figure 3: The power spectrum $P_{R}(k)$ of scalar (curvature) perturbations in
our model: numerically calculated from the Mukhanov-Sasaki equation (in red)
and from Eq. (12) in the slow-roll approximation (in black).

The masses of generated PBH can be estimated from the peak as follows [25]:

$M_{\rm PBH}\simeq{{M_{\rm Pl}^{2}}\over{H(t_{\rm peak})}}\exp\left[2(N_{\rm
total}-N_{\rm peak})+\int_{t_{\rm peak}}^{t_{\rm
total}}\epsilon(t)H(t)dt\right]~{}~{}.$ (13)

The right-hand-side of this equation is most sensitive to the difference
$(N_{\rm total}-N_{\rm peak})$, while the integral is not very sensitive to
the detailed shape of the power spectrum and contributes merely the sub-
leading correction. We numerically solved the MS equation, obtained the exact
power spectrum and compared it with the result of the slow-roll approximation
based on Eq. (12), see Fig. 3. The exact (MS-based) peak is the three orders
of magnitude higher.

We summarize our findings in the Table where we give the values of the key
observables $n_{s}$, $r$ and $M_{\rm PBH}$ in our model, with the fine-tuned
parameters $b=2.89$ and $g=0.41$.  $\phi_{in}/M_{\rm Pl}$ $n_{s}$ $\Delta N$
$r$ $N_{\rm total}$ $M_{\rm PBH}$ 5.099 0.95657 21 0.00532 65.03 $2.17\cdot
10^{19}$ g 5.12 0.95737 20 0.005139 65.09 $5.72\cdot 10^{18}$ g 5.146 0.95831
19 0.004922 65.03 $8.41\cdot 10^{17}$ g 5.17 0.95915 18 0.004730 65.02
$1.43\cdot 10^{17}$ g 5.095 0.95646 16 0.0053565 60.01 $3.89\cdot 10^{15}$ g
5.12 0.95738 15 0.005139 60.01 $7.38\cdot 10^{14}$ g

It follows from the Table that the value of the tensor-to scalar ratio stays
well inside the current observational bound, $r<0.036$ ($2\sigma$). The value
of the index $n_{s}$ of scalar perturbations agrees with the PLANCK
measurements [13, 14, 15],

$n_{s}=0.9649\pm 0.0042\quad(1\sigma)~{},$ (14)

within $2\sigma\div 3\sigma$.

In order to form PBH with masses beyond the Hawking evaporation limit of
$10^{15}$ g, so that those PBH can survive in the present universe and form
CDM, it is crucial to have the duration $\Delta N$ of the second stage of
inflation for 18 e-folds at least or higher, up to 21. The $M_{\rm PBH}$ grows
with increasing $\Delta N$, but consistency (within $3\sigma$) with the
measured value of $n_{s}$ is lost beyond $\Delta N=22$. It is also clear from
the Table that the total duration of inflation should be as long as 65
e-folds.

##  5 Conclusion

In this letter we adapted the $F(R)$ gravity model [22] for double inflation
and PBH production. The significant enhancement of the power spectrum due to
large scalar perturbations, namely, the large peak with the height about
$10^{7}$ times larger the CMB reference level is apparently possible. The PBH
resulting from gravitational collapse of primordial curvature perturbations
can have masses $10^{17}\div 10^{19}$ g, so that they can also survive in the
present universe and may form part of CDM. These PBH masses are very light and
smaller than those obtainable in the multi-field models of inflation, see
e.g., Ref. [24].

To get agreement with the CMB measurements within $3\sigma$, we tuned the
parameters $(b,g,R_{0})$. The $n_{s}$ becomes smaller, and $r$ becomes larger,
when compared to the original Starobinsky model.

When compared to the results of Ref. [16], our peak in Fig. 3 is equally high
but is narrower. When compared to the results of Ref. [21], our PBH masses
(13) in the Table are significantly larger.

It is possible to further generalize the master $F$-function by including a
phenomenological description of dark energy, simply by using the Appleby-
Battye-Starobinsky ansatz (2) twice, with more parameters $(b_{\rm vac},g_{\rm
vac})$ and $\sqrt{R_{\rm vac}}\sim 10^{-33}$ eV.

Fine tuning of the parameters in our model is needed to get a significant
enhancement of the power spectrum of scalar perturbations leading to the PBH
with masses beyond the Hawking evaporation limit and, hence, the possible PBH
DM at present. A significant change of the parameters would lead to a
significant reduction of the power spectrum enhancement with the masses of
emerged PBH below the Hawking limit and no PBH DM. However, the remnants of
those PBH may still form DM [4].

Our approach is entirely phenomenological and classical. Considerations of
quantum corrections require another framework and are beyond the scope of this
investigation (see, however, Refs. [26, 27]).

## Acknowledgements

DF and SS were supported by Tomsk State University. DF was also supported by
the Foundation for the Advancement of Theoretical Physics and Mathematics
”BASIS”. SVK was supported by Tokyo Metropolitan University, the Japanese
Society for Promotion of Science under the grant No. 22K03624, the World
Premier International Research Center Initiative (MEXT, Japan), and the Tomsk
Polytechnic University development program Priority-2030-NIP/EB-004-0000-2022.
SS was also supported by the Huawei Technologies.

The authors are grateful to Cristiano Germani, Maxim Khlopov and Ranjan Laha
for correspondence.

## References

  * [1] I. Novikov and Y. Zeldovic, “Cosmology,” Ann. Rev. Astron. Astrophys. 5 (1967) 627–649.
  * [2] S. Hawking, “Gravitationally collapsed objects of very low mass,” Mon. Not. Roy. Astron. Soc. 152 (1971) 75.
  * [3] A. Dolgov and J. Silk, “Baryon isocurvature fluctuations at small scales and baryonic dark matter,” Phys. Rev. D 47 (1993) 4244–4255.
  * [4] J. D. Barrow, E. J. Copeland, and A. R. Liddle, “The Cosmology of black hole relics,” Phys. Rev. D 46 (1992) 645–657.
  * [5] B. J. Carr, “Primordial black holes as a probe of cosmology and high energy physics,” Lect. Notes Phys. 631 (2003) 301–321, arXiv:astro-ph/0310838.
  * [6] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, “Primordial black holes—perspectives in gravitational wave astronomy,” Class. Quant. Grav. 35 no. 6, (2018) 063001, arXiv:1801.05235 [astro-ph.CO].
  * [7] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “Constraints on primordial black holes,” Rept. Prog. Phys. 84 no. 11, (2021) 116902, arXiv:2002.12778 [astro-ph.CO].
  * [8] B. Carr and F. Kuhnel, “Primordial Black Holes as Dark Matter: Recent Developments,” Ann. Rev. Nucl. Part. Sci. 70 (2020) 355–394, arXiv:2006.02838 [astro-ph.CO].
  * [9] A. A. Starobinsky, “A new type of isotropic cosmological models without singularity,” Phys. Lett. B 91 no. 1, (1980) 99 – 102.
  * [10] S. V. Ketov, “On the equivalence of Starobinsky and Higgs inflationary models in gravity and supergravity,” J. Phys. A 53 no. 8, (2020) 084001, arXiv:1911.01008 [hep-th].
  * [11] S. V. Ketov and A. A. Starobinsky, “Embedding $(R+R^{2})$-Inflation into Supergravity,” Phys. Rev. D 83 (2011) 063512, arXiv:1011.0240 [hep-th].
  * [12] S. V. Ketov and A. A. Starobinsky, “Inflation and non-minimal scalar-curvature coupling in gravity and supergravity,” JCAP 08 (2012) 022, arXiv:1203.0805 [hep-th].
  * [13] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO].
  * [14] BICEP, Keck Collaboration, P. A. R. Ade et al., “Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season,” Phys. Rev. Lett. 127 no. 15, (2021) 151301, arXiv:2110.00483 [astro-ph.CO].
  * [15] M. Tristram et al., “Improved limits on the tensor-to-scalar ratio using BICEP and Planck,” arXiv:2112.07961 [astro-ph.CO].
  * [16] J. Garcia-Bellido and E. Ruiz Morales, “Primordial black holes from single field models of inflation,” Phys. Dark Univ. 18 (2017) 47–54, arXiv:1702.03901 [astro-ph.CO].
  * [17] C. Germani and T. Prokopec, “On primordial black holes from an inflection point,” Phys. Dark Univ. 18 (2017) 6–10, arXiv:1706.04226 [astro-ph.CO].
  * [18] C. Germani and I. Musco, “Abundance of Primordial Black Holes Depends on the Shape of the Inflationary Power Spectrum,” Phys. Rev. Lett. 122 no. 14, (2019) 141302, arXiv:1805.04087 [astro-ph.CO].
  * [19] I. Dalianis, A. Kehagias, and G. Tringas, “Primordial black holes from $\alpha$-attractors,” JCAP 01 (2019) 037, arXiv:1805.09483 [astro-ph.CO].
  * [20] H. V. Ragavendra, P. Saha, L. Sriramkumar, and J. Silk, “Primordial black holes and secondary gravitational waves from ultraslow roll and punctuated inflation,” Phys. Rev. D 103 no. 8, (2021) 083510, arXiv:2008.12202 [astro-ph.CO].
  * [21] L. Iacconi, H. Assadullahi, M. Fasiello, and D. Wands, “Revisiting small-scale fluctuations in $\alpha$-attractor models of inflation,” arXiv:2112.05092 [astro-ph.CO].
  * [22] S. A. Appleby, R. A. Battye, and A. A. Starobinsky, “Curing singularities in cosmological evolution of F(R) gravity,” JCAP 06 (2010) 005, arXiv:0909.1737 [astro-ph.CO].
  * [23] V. R. Ivanov, S. V. Ketov, E. O. Pozdeeva, and S. Y. Vernov, “Analytic extensions of Starobinsky model of inflation,” JCAP 03 no. 03, (2022) 058, arXiv:2111.09058 [gr-qc].
  * [24] S. V. Ketov, “Multi-Field versus Single-Field in the Supergravity Models of Inflation and Primordial Black Holes,” Universe 7 no. 5, (2021) 115.
  * [25] S. Pi, Y.-l. Zhang, Q.-G. Huang, and M. Sasaki, “Scalaron from $R^{2}$-gravity as a heavy field,” JCAP 05 (2018) 042, arXiv:1712.09896 [astro-ph.CO].
  * [26] L.-H. Liu, T. Prokopec, and A. A. Starobinsky, “Inflation in an effective gravitational model and asymptotic safety,” Phys. Rev. D 98 no. 4, (2018) 043505, arXiv:1806.05407 [gr-qc].
  * [27] S. V. Ketov, “Starobinsky-Bel-Robinson gravity,” arXiv:2205.13172 [gr-qc].


# Spectral properties of Compton inverse radiation: Application of Compton
beams

Eugene Bulyak1 and Junji Urakawa2 1 NSC KIPT, Kharkov, Ukraine  2 KEK,
Tsukuba, Ibaraki, Japan<EMAIL_ADDRESS>
###### Abstract

Compton inverse radiation emitted due to backscattering of laser pulses off
the relativistic electrons possesses high spectral density and high energy of
photons – in hard x–ray up to gamma–ray energies – because of short wavelength
of laser radiation as compared with the classical electromagnetic devices such
as undulators.

In this report, the possibility of such radiation to monochromatization by
means of collimation is studied. Two approaches have been considered for the
description of the spectral–angular density of Compton radiation based on the
classical field theory and on the quantum electrodynamics. As is shown, both
descriptions produce similar total spectra. On the contrary, angular
distribution of the radiation is different: the classical approach predicted a
more narrow radiation cone. Also proposed and estimated is a method of the
‘electronic’ monochromatization based on the electronic subtraction of the two
images produced by the electron beams with slightly different energies. A
‘proof–of–principle’ experiment of this method is proposed for the LUXC
facility of KEK (Japan).

##  1 Introduction

Sources of Compton radiation, in which photons of intense laser pulse
scattered off from relativistic electrons, are able to produce bright x–ray
beams with narrow bandwidth. The process of Compton scattering can be treated
as two–particle elastic scattering: one of the particle is represented by the
electron, another one by the photon. Since the laser photon has negligibly
small energy as compared with electron’s, the recoil of electron is negligibly
small – the photon is scattered off within a narrow cone along electron’s
trajectory.

Such sources have a substantial potential for applications in different areas
of medicine, biology, physics, etc. This potential is emerged not only from
its brightness, but tunability and ‘quasimonochromaticity’ of the spectrum.

The report highlights spectral properties of Compton sources and presents a
potential application of these sources for x–ray angiography, which possesses
a substantial advantage over the conventional methods.

##  2 Spectrum of collimated Compton inverse radiation

The spectrum of inverse Compton radiation is determined by two relations, (i)
kinematic dependence of the energy of a scattered off quantum,
$E_{\mathrm{x}}$ upon the crossing angle $\phi$ between the laser photon and
the electron, energy of both the electron and the photon, and the scattering
angle $\psi$, and (ii) the differential cross section – dependence of
probability of scattering the photon at angle $\psi$. Within the small–angle
approximation, $\psi\ll 1$ (the crossing angle $\phi=0$ corresponds to the
head–on collision), neglected the electron recoil, these relations read ([1]):

$E_{\mathrm{x}}=\frac{2\gamma^{2}(1+\cos\phi)E_{\mathrm{las}}}{1+\gamma^{2}\psi^{2}}\;;\qquad\mathrm{d}\sigma=8\pi
r^{2}_{0}\frac{\psi\gamma^{2}\left(1+\gamma^{4}\psi^{4}\right)}{\left(1+\gamma^{2}\psi^{2}\right)^{4}}\mathrm{d}\psi\;,$

where $\gamma=E_{\mathrm{e}}/m_{\mathrm{e}}c^{2}$ is the Lorentz factor of the
electron, $r_{0}$ is the classical electron radius.

Convolution of the cross section with distribution functions of both the
electron bunch and the photon pulse taking into account the kinematic relation
will produce a real collimated spectrum of the radiation. As it follows from
the kinematics, the width of spectrum is determined by the collimation angles
$\psi_{i}\leq\psi\leq\psi_{f}$, as is depicted in Fig.2. The spectral–angular
density of the collimated Compton radiation for an ideal case (monoenergetic
electrons and laser photons with parallel trajectories) is presented in Fig.2.

Figure 1: Scheme of the Compton radiation collimation.

Figure 2: Spectral-angular density of the Compton radiation.

For the vast majority of practical sources the inner collimation angle is
zero, $\psi_{i}=0$. We will consider below this particular case of an iris
collimator.

In general, the spectrum is dependent on the following factors: the collimator
opening angle, the energy spread of electrons and photons and the angular
spread of their trajectories at the interaction point (IP). In practical
cases, the angular spread of photon ‘trajectories’ produces least impact on
the spectrum, as it can be seen from the kinematics: Its effect is
proportional to the spread squared because IP is usually set up at the laser
waist, while the electron spread impact is $\gamma^{2}$ times larger. Also the
spread of electron energy within the bunch usually much higher than that of
photons. Therefore below we will consider effects of the electron bunch phase
volume – spreads of energy and trajectories – upon the collimated spectrum of
scattered off laser quanta.

###  2.1 Spread of electrons energy

The energy spread dilutes the spectral–angular density curve in Fig.2 along
the energy axis. Since the energy spread in bunches of Compton sources is
usually small – from a few percents down to fraction of a percent, – the
corresponding partial energy spread of x–rays is equal to doubled the
electrons’: For the normal (Gaussian) distribution of electrons’ energy in the
bunch, the ‘pin-hole’ collimated x–ray beam has the doubled reduced
dispersion: $\sigma_{x}/E_{x}\approx 2\sigma_{e}/E_{e}$. (Detailed study on
the spectrum for head-on collision is presented in [2].) Analytical collimated
spectra for finite energy spread of electrons with parallel trajectories are
presented in Fig.4.

Figure 3:  X-ray spectra for the electron energy spreads
$\sigma_{e}\gamma=(0,0.01,0.02,0.03)$ into collimating range
$\psi\gamma=(0,0.1)$.

Figure 4:  X-ray energy spectra for the electron trajectories spreads
$\sigma_{e}\gamma=(0.01,0.1,0.5,1.0)$ into collimating range
$\psi\gamma=(0,0.1)$.

###  2.2 Angular spread of electrons trajectories

The angular spread of electrons’ trajectories induces much more widening into
the spectrum since the radiation is emitted within narrow cone along the
trajectory. (It should be noted that a quantum emitted at the angle $1/\gamma$
to the electron trajectory has a half of the energy of that emitted along the
trajectory.) Also tight focusing of the bunches at the interaction point aimed
to gain the yield of x rays, increases the angular spread.

A collimated spectrum with account for the angular spread has asymmetric and
more complicated shape than that due to the energy spread. It is specified
with steep high–energy cutoff (stemmed from the energy conservation law) and
long low–energy tail, see Fig.4.

Examples of simulated spectra for different collimation opening angles is
presented in Fig.5 (the parameters for simulation were taken similar to those
of LUCX facility, see [3]).

Figure 5:  Spectra for collimation $\psi\gamma=0.2,0.4,0.6,0.8,1.0$ (from the
narrowest to widest profiles). Angular horizontal and vertical spreads in the
bunch are $\gamma\sigma_{x,y}=0.08,0.32$.

From the study on the collimated x-ray spectra generated by Compton sources
the points are followed:

  * •

The spectrum width is limited by the high energy cutoff, no photons with
higher energy.

  * •

The spectrum has the maximum close to the high-energy cutoff.

  * •

High-energy cutoff of the spectrum is steep.

  * •

The collimator opening angle should not be more narrow than the (approximately
doubled) angular spread of electron trajectories at IP. More narrow
collimation will reduce the spectral maximum but not the width.

##  3 Application of Compton sources for angiography

Angiography, a medical x–ray based imaging technique to visualize the inside
of blood vessels and organs of the body, uses a radio-opaque contrast agent
injected into the blood vessel. The difference between x–ray images without
and with the contrast agent in the blood produces a picture of the blood flow.

Iodine–based radiocontrast agents are most commonly used in angiography with
maximum contrast energy of x–rays just above the K–edge of iodine, 33.17 keV.
Hence, the x–rays within a band of $\sim$30–35 keV most effectively produce
the images. Other energies are redundant and even detrimental since they cause
additional radiation load.

The Compton x–ray sources with the tunable energy peak possess certain
advantages over the conventional bremsstrahlung sources. In Fig.6, there an
ideal Compton spectrum and a bremsstrahlung one are presented that produce
equal fluxes within 30–35 keV energy range.

Figure 6: Bremstrahlung (70 keV electrons and tungsten radiator) and Compton
spectra. the vertical lines limit a useful x-ray energy band, the gray areas
beyond show redundant x–ray radiation load from the bremsstrahlung spectrum.

As it can be seen from the picture, the conventional source produces much more
background x–ray quanta (the gray areas beyond 30–35 keV range).

###  3.1 Simulation of Compton angiography

For verification of Compton sources applicability for angiography, a digital
model has been created. The models simulates the image of the tissue composed
from the skeletal muscle, the bone cortical, blood and iodine components
(x–ray attenuation coefficients were taken from NIST [5]). Distribution of the
energy of x–ray quanta is simulated with Monte Carlo method according to the
collimate ideal Compton spectrum. Distribution of the impinging quanta over
the tissue face is uniform, random.

A run of simulation on the model reveals applicability of Compton x–ray source
for angiography as is presented in Fig.7 (left). The model input data was:
x–ray energy ranges 23.3–35 keV; number of quanta in the range $2\times
10^{7}$ (total number over the spectrum $2\times 10^{7}$); the sensor mesh
$100\times 100$ pixels (1 mm2 pixels on $10\times 10\,\mbox{cm}^{2}$ tissue);
the surface densities of tissue’s substances are as follows (in g/cm2): muscle
5, bone 0.5, blood 0.5, iodine 0.0125. The muscle component is uniform density
over the tissue, the bone–like substance is added at the bottom ($2<y<4$ cm),
the blood one at the top of the tissue ($7<y<8$ cm), iodine is placed in the
left half of the blood tape ($x<5$ cm).

Figure 7:  X–ray image with $E_{x}^{max}=35\,\mbox{keV}$ (left) and a
subtracted ‘30 keV’ minus ‘35 keV’ image (right).

To compare the left half of ‘35 keV’ image in Fig.7 (blood with the contrast
agent) with the right half (no contrast agent), one can detect location of the
iodine. Total flux density for the simulated case is $4\times 10^{5}$ quanta
per cm2.

###  3.2 Advanced method of Compton angiography

Specific properties of the Compton x–ray radiation spectrum – steep high-
energy cutoff with the maximum near it, and strong dependence of the maximal
energy upon the energy of electrons, – enable us to propose an advanced method
of Compton angiography.

The method consists in following. Dissimilar to the conventional angiography
where location of the contrast agent is revealed as a result of subtracting
the image without the agent from that with it, the both images are taken with
presence of the contrast agent but one of them at the maximum energy of the
spectrum tuned just below the K–edge of the agent, another – above it.
Difference between these images will display the contrast agent location.

Simulations of the method have validated it. In Fig.7 (right panel) a
subtracted image is presented, the image of the described above tissue at the
maximum energy of Compton x–rays of 35 keV (left panel) was subtracted from
the image with the maximum energy 30 keV. Since opaqueness of the agent
increases with increasing the energy above K–edge while the other components
become more transparent, the location of agent in the difference image is
lighter than the others (see Fig.7, right panel).

Thus, the advanced method is capable to reveal the location of the contrast
agent much faster than the conventional Compton angiography method with about
the same flux of x–ray quanta. Speed up of the imaging is due to the fact that
injection of the contrast agent into the blood vessel will take much more time
than switching of the accelerator energy (we believe it will take a fraction
of second).

###  3.3 ‘Proof–of–principle’ experiment proposal for LUCX

The ‘proof–of–principle’ experiment must verify the basic suggestion for
Compton x–ray source: opaqueness of a contrast agent increases with increasing
the maximal energy of spectra above the K-edge of contrast agent. Also ability
of Compton sources to gather sufficient statistics revealing the location of
the contrast agent should be proved.

Such an experiment is proposed to conduct at LUCX facility [3, 4]. Since the
energy of Compton x–ray photons in this accelerator with YAG laser (1 eV) can
not overlap the K–line of iodine, we propose to use bromine as a contrast
agent with K-edge at 13.4737 keV, see [5].

We have simulated a ‘proof–of–principle’ experiment for LUCX. Results of the
simulation is depicted in Fig.8. The sample of $20\times 20$ cells was
irradiated by $2\times 10^{6}$ quanta of each energy,
$E_{x}^{(max)}=13,16\,\mbox{keV}$ (the images ‘A’ and ‘B’, respectively). The
Compton spectrum was collimated to $1/\gamma$ opening angle. The sample
thickness (muscle skeletal) was 0.5 g/cm2 with additional 0.15 g/cm2 in top
left corner, bromine thickness was 10 mg/cm2 (bottom left corner).

Figure 8:  ‘13 keV’ image (A, left), ‘16 keV’ image (B, right) and the
subtracted density (A-B, centre).

As is follows from the images, one cannot tell from single run (A and B
images) which of the shadows induced by bromine or by an additional thickness.
The bromine shadow is thinner when the spectrum does not overlap the K–peak,
‘A’ image. The subtracted image, ‘A-B’ clearly indicates which of shadows
belong to the additional thickness or to the bromine – the contrast agent is
more transparent for lower energy inverse to the regular behavior.

##  4 Summary and conclusion

X–ray radiation generated by the Compton sources has the spectrum with the
steep high–energy cutoff, which is independent on the collimation opening
angle. The specific shape of the spectrum with tunable maximum is able to
substantially reduce radiation load for the angiography procedure.

We propose a novel angiography procedure consisting in deriving a subtracted
image from the two images having been made with different maximum energy of
the Compton spectrum, one with the energy below the K–edge of a contrast
agent, another – with the maximum energy overlapped it.

The carried out simulations proved the suggestion of advantage the Compton
subtracted scheme for angiography. Also we proposed a ‘proof–of–principle’
experiment that can be conducted at LUCX facility of KEK employing bromine as
a radiocontrast agent.

Authors wishing to acknowledge Prof. A. Dovbnya, Drs. P. Gladkikh and V.
Skomorokhov for their assistance and helpful discussions. Work is supported by
the Photon and Quantum Basic Research Coordinated Development Program by the
Ministry of Education, Culture, Sports, Science and Technology ‘Fundamental
Technology Development for High Brightness X-ray Source and the Imaging by
Compact Accelerator.’

## References

## References

  * [1] Bulyak E and Skomorokhov V 2004 Proc. Eur. Part. Accel. Conf. 2004 report THPKF063
  * [2] Sun C, Li J, Rusev G, Tonchev A P and Wu Y K 2009 Phys. Rev. ST-AB 12 062801
  * [3] Fukuda M, Araki S, Aryshev A, Honda Y, Terunuma N, Urakawa J, Sakaue K and Washio M 2013 Proc. NAPAC’13 report TUPMA01
  * [4] Sakaue K, Washio M, Fukuda M, Honda Y, Terunuma N, Urakawa J 2013 Proc. IPAC 2013 report WEPWA017
  * [5] http://www.nist.gov/pml/data/xraycoef/


#  $IRIS^{2+}$: A comprehensive database of stratified thermodynamic models in
the low solar atmosphere.

Alberto Sainz Dalda Lockheed Martin Solar & Astrophysics Laboratory, 3251
Hanover Street, Palo Alto, CA 94304, USA  Bay Area Environmental Research
Institute, NASA Research Park, Moffett Field, CA 94035, USA.  Aaryan Agrawal
Henry M. Gunn High School, 780 Arastradero Rd, Palo Alto, CA 94306  Bart De
Pontieu Lockheed Martin Solar & Astrophysics Laboratory, 3251 Hanover Street,
Palo Alto, CA 94304, USA  Rosseland Center for Solar Physics, University of
Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway  Institute of Theoretical
Astrophysics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway
M. Gošić Lockheed Martin Solar & Astrophysics Laboratory, 3251 Hanover Street,
Palo Alto, CA 94304, USA  Bay Area Environmental Research Institute, NASA
Research Park, Moffett Field, CA 94035, USA.

###### Abstract

Our knowledge of the low solar atmosphere, i.e., the photosphere and
chromosphere, is based on the knowledge gained from the observations and the
theoretical and numerical modeling of these layers. In this sense, the
thermodynamical and magnetic semi-empirical models of the solar atmosphere
have significantly contributed to the advance in the understanding of the
physics of the Sun. In the past, many of these models have been used as a
reference that helps us to, e.g., constrain the theoretical and numerical
modeling, or to verify the goodness of physical parameters obtained from the
inversion of the spectral lines. Nevertheless, semi-empirical models are quite
limited by the assumptions that are inherent to the approach and do not
necessarily provide an accurate view of the instantaneous and local
thermodynamic conditions in the solar atmosphere. In this work, we provide an
extensive collection of thermodynamic model atmospheres for active regions
(ARs) obtained from the simultaneous inversion of 6 lines sensitive to changes
in the thermodynamical conditions in the chromosphere, and another 6 lines
sensitive to changes in the thermodynamical conditions in the photospere.
These inversions were made using 320 representative profiles (RP) obtained by
clustering the profiles in the umbra, penumbra, pore-like, plage, and
surrounding quiet-sun in 126 active regions. Due to the simultaneous inversion
of the selected lines, the resulting representative model atmosphere (RMA)
samples the thermodynamics from the bottom of the photosphere to the top of
the chromosphere. As a result, this database, named $IRIS^{2+}$ and formed by
$40,320$ RP-RMA pairs, represents the most comprehensive collection of
stratified-in-optical-depth thermodynamic models of the low solar atmosphere.

Sun — photosphere — chromosphere — radiative transfer

††software: $IRIS^{2}$

##  1 Introduction

Modern astrophysics started when we were able to derive physical information -
beyond morphological and phenomenological behavior - from the objects we
observed in the sky. We can establish the last half of the $19^{th}$ century
as the birth of what we now call astrophysics. Two critical techniques were
needed: spectroscopy and photography, i.e.: the capability to record the light
coming from the Sun (or any other object in the sky) as a function of
wavelength. The comparison of data thus obtained with similar data obtained in
the laboratory has led to great progress in understanding the physical
conditions on stars. Key to the interpretation of these new measurements was
the development of models that led to a greater understanding of matter, i.e.,
the atomic model, and light, more precisely, electromagnetic radiation. A
quintessential example of this synergy between experiments, techniques and
theory is the discovery of the magnetic field in the Sun(spots), inferred
through the Zeeman effect in the spectral lines observed on sunspots by Hale
(1908). Since then, many advances, both theoretical and technological have
improved our knowledge about the matter, the electromagnetic radiation, and
the physics of the objects observed in the sky.

The study of the Sun has played a prominent role in the knowledge we now have
about stars, and in the development of the techniques and ideas needed to gain
that knowledge. Thanks to its proximity, we are able to resolve many
structures on its surface (photosphere), and the rest of the solar atmosphere,
i.e., the chromosphere (Rutten, 2007; Cauzzi et al., 2008; Wedemeyer et al.,
2016; Carlsson et al., 2019), transition region, and the corona. Current
instrumentation allows us to recover in great detail the physical conditions
of solar structures as small as $100~{}km$, with a temporal cadence of a few
seconds, and with a high sensitivity both in the magnetic field
($\delta|B|<10~{}G$) and in other derived thermodynamic parameters (e.g.,
$\delta T<50~{}K$).

There are two main ways to gain knowledge about the solar atmosphere. Forward
modeling involves numerical models of the thermodynamic (or
magnetothermodynamic) structure and comparing the synthetic spectral lines or
profiles obtained by solving the radiative transfer equation (RTE) with solar
observations. This method is not unique as it starts with an initial guess at
the thermodynamic parameters and assumptions about the dominant physical
processes. The second method recovers a model atmosphere from iteratively
fitting synthetic profiles - obtained by solving the RTE from a initial guess
model atmosphere and slight subsequent modifications - until a best fit is
found between the observed profile and the synthetic profile. This approach is
also not unique as it sensitively depends on the initial guess and could
suffer from common issues with least-square fitting such as local minima in
the loss function. Because of that, we (should) say that the recovered model
atmosphere is compatible with the observed profile, and therefore that these
physical conditions may occur in the solar atmosphere. This last tool is
called the inversion of spectro(polarimetric) data111The reader interested in
a detailed review of this method can be found in del Toro Iniesta & Ruiz Cobo
(2016).

While neither method provides a unique solution, both can provide valuable
insight into the physical conditions and processes that dominate the solar
atmosphere. For both methods it is important to have an accurate knowledge of
the model of the solar atmosphere. In the forward modeling, this knowledge
comes from the physics recovered from previous observations and by the
theoretical considerations of the problem (e.g., hydrostatic equilibrium,
ambipolar difusion, Khomenko et al. 2014; Martínez-Sykora et al. 2020), which
help us constrain the physical conditions considered in the model. For the
inversion, the feasibility of the initial guess model used to start the
inversion is critical, as well as the theoretical considerations made while
solving the RTE (e.g., statistical equilibrium, non-local thermodynamic
equilibrium, partial frequency distribution, e.g., Leenaarts et al. 2013), and
the atom model used to synthesize the spectral lines. In fact, the theoretical
considerations should be similar in both cases, since both methods try to
reproduce the observables we detect, i.e., the spectral line profiles. Both
methods suffer from computational and algorithmic limitations, whether it is
the inclusion of the effects of multi-fluid interactions (forward modeling,
see Wargnier et al. 2022) or the calculation of 3D radiative transfer in
chromospheric lines (inversions, see Sukhorukov & Leenaarts 2017; Bjørgen et
al. 2019).

In this paper, we focus on providing a collection of thermodynamic model
atmospheres for active regions (ARs). We consider an active region as the non-
exclusive combination of umbra, pores or pore-like structures, penumbra, plage
and the neighboring quiet Sun. There are many solar atmosphere models
available222An excellent set of solar atmosphere models has been collected by
Dr. Basilio Ruiz Cobo. They are available at
https://github.com/BasilioRuiz/SIR-code/tree/master/models., most of them are
particularly devoted to one solar feature, e.g., the HSRA quiet-sun model
(Gingerich et al., 1971), a hot sunspot model (Collados et al., 1994), a
penumbra model (del Toro Iniesta et al., 1994), or a plage model (Model1004
Fontenla et al., 2009). Some of these models are considered as the standard
model for their corresponding features. In the case of the inversions, these
models are mostly used as initial guess models. The results presented in this
paper complement these models by providing a set of $40,320$ representative
model atmosphere (RMA) and their corresponding RPs, which are labeled with the
feature where they were observed, e.g. “plage”, their location on the solar
disk, and the recording date. Therefore, a more diverse set of models and
profile are given, and they can be used to constrain numerical models, as
initial model guesses for inversions, or for the synthesis of profiles at
different wavelengths.

In Section 2 we detail the main characteristics of the IRIS data that we have
used in our study, which spectral lines we have selected for the inversions,
and where in the solar atmosphere these lines are sensitive to changes in the
thermodynamic conditions. We explain how we have selected the different areas
in the AR data, i.e. umbra, pore-like, penumbra, plage, and surrounding quiet
Sun. We justify the need to cluster these data and how we treated the spectral
lines to be simultaneously inverted. In Section 3 we show the results of the
inversions, i.e. the fit between the observed profiles and the synthetic
profile provide by the inversion, as well as the model recovered from that
inversion. As a result, we present a database of synthetic RPs and RMAs in the
low solar atmosphere. This data base is called $IRIS^{2+}$, since it
represents an extension - in terms of the optical depth where the models are
now more accurate - to $IRIS^{2}$. In Section 4 we discuss the confidence and
the usability of $IRIS^{2+}$ database.

##  2 Data

The most critical part of this project is likely to select a set of data that
represents properly the variety of active regions in the Sun. As we have
already mentioned, by active region we understand the non-exclusive
combination of plage, umbra, pore, and/or penumbra with surrounding quiet Sun
(which may or may not be affected by the neighboring active region). Thus, an
IRIS data set containing spectra from a region that combines only quiet sun
and a plage region is considered in this selection as an active region. This
relaxed definition does not diminish the goal of this study: to have a
comprehensive collection of RPs and RMAs that characterizes the thermodynamics
of the main elements of the ARs. To this aim, we have selected 126 IRIS data
sets located at different positions on the solar disk, observed with different
exposure times, and containing different active regions observed in the time
range from July 2013 (the beginning of the IRIS mission) to 2021. Figure 1
shows the distribution of the selected ARs with respect the observing date
(left panel), the location on the solar disk given by $\mu=cos(\theta)$, with
$\theta$ the heliocentric angle (center panel), and the exposure times (right
panel). The size of the field of view (FoV) was $128\times 130~{}arcsec^{2}$
for 104 data sets, $128\times 175~{}arcsec^{2}$ for 14 data sets, $140\times
175~{}arcsec^{2}$ for 3 data sets, $112\times 124~{}arcsec^{2}$ for 3 data
sets, $64\times 175~{}arcsec^{2}$ for 1 data set, and $64\times
124~{}arcsec^{2}$ for 1 data set.

Figure 1: Distribution of the selected data in the $IRIS^{2+}$ database with
respect the observation date (left), the position in the solar disk given by
$\mu=cos(\theta)$ \- with $\theta$ the heliocentric angle (middle), and the
exposure time (right).

###  2.1 Selection of spectral lines

IRIS | Fe I | Ni II | C I | Fe I | Mg II h | Mg II UV3&2 | Mg II k | Fe I | Mg II UV1 | Ti II
---|---|---|---|---|---|---|---|---|---|---
linelist | 2827.33 | 2815.18 | 2810.58 | 2809.15 | 2803.53 | 2798.85&.72 | 2796.35 | 2793.22 | 2791.60 | 2785.46
v36_00 & v38_00 |  | ✓ |  | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |
v36_01 & v38_01 |  | ✓ |  |  | ✓ | ✓ | ✓ |  |  |
v36_02 & v38_02 |  | ✓ |  |  |  | ✓ | ✓ |  |  |
v36_03 & v38_03 | ✓ | ✓ |  |  | ✓ | ✓ | ✓ | ✓ | ✓ |
v36_04 & v38_04 | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
v40_00 |  |  |  |  | ✓ | ✓ | ✓ |  |  | ✓
v40_02 |  |  |  |  | ✓ |  | ✓ |  |  |
v40_09 | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
Table 1: Availability of the spectral lines (in Å) considered in the
$IRIS^{2+}$ for each observation IRIS linelist. The C II 1334 & 1335 Å lines,
not included in the table, are always included in any IRIS linelist,
therefore, they are always considered by $IRIS^{2+}$.

Another important issue is the selection of the photospheric spectral lines.
We have tried different combinations of photospheric lines located in the near
ultra-violet (NUV) spectral range observed by IRIS. The bottom panel in Figure
2 shows the full spectrum taken by IRIS in its NUV channel. In the top panel
of this figure, an average profile of the quiet sun observed at the center of
the solar disk is shown. The selected photospheric lines are indicated by
vertical lines (orange in the bottom, black in the top panel), while the
chromospheric lines, i.e. the Mg II h&k lines and the Mg II UV triplet lines,
are similarly indicated in violet. The full spectrum along the slit that is
shown in the bottom panel of this figure is intended to help the comparison
with the IRIS linelist images333The IRIS linelist of an IRIS observation is
shown in the ”Raster” column that shows up after an observation has been
selected in the right panel at https://iris.lmsal.com/search/. that illustrate
the combination of spectral ranges selected in an IRIS observation. Because of
these linelists (which are used for telemetry reasons, i.e., to minimize the
data volume that is downlinked), we have to take into account the availability
of the selected lines in the different IRIS observation setups or linelists.
Table 1 shows which selected lines are available in the IRIS linelists. Note
that in our study, in addition to the lines observed in the NUV channel, we
also consider the C II 1334 & 1335 Å lines recorded in the far ultraviolet
(FUV) channel.

Figure 2: Top: The photospheric (orange) and chromospheric (violet) selected
lines in the spectral range observed by IRIS around the Mg II h&k lines. The
values in parenthesis indicate the approximate height (in $Mm$) for an optical
depth of 1, and the uncertainty associated with the fact that the formation
occurs in a corrugated layer, as obtained by Pereira et al. (2013). Bottom:
The spectrum as seen on the detector is shown for a better comparison with the
IRIS linelist files, e.g. v36_01, which show all the spectral ranges selected
for a given observation. Note that wavelength increases towards the left in
this panel.

The selection of the photospheric lines is intended to recover as accurately
as possible the variation of thermodynamics from the bottom to the top of the
photosphere. This information is derived from the response function (RF, Mein,
1971; Landi Degl’Innocenti & Landi Degl’Innocenti, 1981) of the intensity
profiles to a variation in a given physical parameter, and from the height at
which the optical depth is above 1, as obtained by Pereira et al. (2013). For
all the photospheric lines in the panels of Figure 2 \- except Fe I 2827.33 Å-
the height is provided in $Mm$. including the range of heights that maps
velocities in the model to observed line shifts, associated with the fact that
the line formation occurs in a corrugated layer. The height estimates are
based on the synthesis of these lines using a modified version of the RH code
(Uitenbroek, 2001; Pereira & Uitenbroek, 2015) on a 3D radiative MHD numerical
simulation obtained by the Bifrost code (Gudiksen et al., 2011) . The RFs of
the intensity (for a multi-line RP in the database) to the variation of the
temperature is shown in the left panel in Figure 3. Note that the optical
depth range where a spectral line is sensitive to the variation of a physical
parameter is slightly different for different solar features (e.g., umbra,
penumbra, filament, plage), and it is also different for the various physical
parameters (temperature, electron density, LOS velocity, microturbulence). The
right panel in Figure 3 shows the optical depth ranges where the selected
lines in the $IRIS^{2+}$ are roughly sensitive to changes in the temperature.
As mentioned above, these ranges summarize the general behavior of the
selected lines for most of the observed solar features in the database. They
should be understood as a reference of where the lines are mostly sensitive to
changes in the thermodynamics.

Figure 3: Left: Response function (RF) of the intensity to variations of $T$
of the selected lines in the $IRIS^{2}$ database. The spectral lines (white)
are shown as a reference. Right: the optical depth ranges where these lines
are on average sensitive to changes of $T$ for the solar features included in
the $IRIS^{2+}$ database. Note that these ranges may slightly vary for each
solar feature and for the other physical parameters in the database. The
values under the lines are given in Å. The values between parenthesis are
given in $Mm$. See Figure 2 for details.

###  2.2 Multi-line representative profiles

Now that we have selected a set of spectral lines that is sensitive to the
thermodynamics from the top of the chromosphere to the bottom of the
photosphere, the next steps is to prepare them for clustering. The need for
clustering our data is due to the computational resources needed to recover
the physical information encoded in the spectral lines, specifically when we
are working with the C II 1334 & 1335 Å and the Mg II h&k lines. As we have
already mentioned the inversion of the latter requires consideration of non-
LTE and partial frequency redistribution of the scattered photons. Therefore,
we follow the same approach introduced by Sainz Dalda et al. (2019), which was
extended by Woods et al. (2021) for the multi-line case of C II 1334 & 1335 Å
lines and Mg II h&k (including Mg II UV2&3 lines), and more recently by Sainz
Dalda & De Pontieu (2022a) for the multi-line case of C II 1334 & 1335 Å
lines, Mg II h&k and all the Mg II UV triplet lines. Now, in addition to the
set of chromospheric lines used by the latter authors we include 6
photospheric lines. Because we are now simultaneously considering lines
sensitive to the thermodynamics at different optical depths we may have RPs
which share the same chomospheric physical conditions but different
photosperic thermodynamic conditions, or the other way around. In other words,
identical chromospheric profiles may be associated with photospheric profiles
that are quite different, or vice versa. We must consider this degeneracy
between the chromospheric set of lines and the photospheric set of lines when
clustering them. Thus, compared to inversions that are based on fewer lines,
we need to consider a larger number of clusters and RPs to capture properly
that degeneracy. We have determined that 320 clusters per data set are able to
represent properly most of the data sets considered in our study within the
hardware limitations we have.

In this investigation, since we are interested in obtaining thermodynamic
models for the ARs we have stratified the data before clustering. That means,
we have first identified those locations where the profiles belong to the
umbra (or pore), the penumbra, the plage, and the quiet Sun. These four
regions are determined by using different intensity thresholds in the
reconstructed intensity map at $2810.58~{}$Åof the data set to be clustered.
The 3 first panels in the top row of Figure 4 show the intensity
spectroheliogram or map for an AR at Mg II k, Ti 2 2785.46 Å and at 2832.04 Å
(photospheric continuum) respectively. The fourth panel shows the areas
identified as umbra, penumbra, plage, and quiet sun. Once these areas are
determined, we impose the number of clusters for each region to be clustered.
The number of clusters for each region is based on the average area covered by
these regions and the variability of different physical conditions in these
areas. For instance, the quiet Sun occupies more area in most of the data
selected compared to the penumbra. On the other hand, the physical conditions
are likely more variable in the latter, although that is not necessary true in
the chromosphere. We have attempted to find a trade-off between the hardware
limitations (to consider at most 320 RPs per data set) and to recover a
meaningful representation of the physics of the selected areas in each data
set. The number of clusters for the umbra, penumbra, plage and quiet-sun are
30, 50, 80, 160 respectively. If there is no umbra, i.e., no umbra has been
detected under the threshold level established, then the features detected are
usually pores, orphan-penumbras, or naked-sunspots. In this case, the number
of clusters for the penumbra are incremented up to 80, and the associated RPs
are labeled as ’pore-like’. If there is neither umbra nor penumbra, the number
of clusters for the plage and quiet Sun are 144 and 176 respectively. In
total, the $IRIS^{2+}$ database has $2,280$ RPs associated with the umbra,
$800$ RPs with the pore-like features, $3,800$ RPs with the penumbra, $12,640$
with plage, and $20,800$ RPs with quiet Sun. To cluster the original IRIS
profiles we created a joint-profile by cropping the selected lines and
concatenating them together. These new profiles are scaled between 0 and 1,
then they are clustered with the number of clusters determined by the criteria
mentioned above. After the clusters are defined, we calculated the RPs
(centroids) as the average of the original profiles, i.e. un-cropped, and
associated with the elements (labels) of each cluster. At this point, the RPs
preserve the original spectral sampling and the original spectral range of the
selected lines. They are ready for the inversion.

###  2.3 Inversion of the Multi-Line Representative Profiles

We have used the multi-atom STockholm inversion Code (STiC de la Cruz
Rodríguez et al., 2016, 2019) to invert the RPs of our selected data.
STiC444STiC is available at https://github.com/jaimedelacruz/stic/. is an MPI-
parallel non-LTE inversion code that utilises a modified version of RH
(Uitenbroek, 2001) to solve the atomic population densities assuming
statistical equilibrium and plane-parallel geometry and it allows including
partial frequency redistribution effects of scattered photons (Leenaarts et
al., 2012). The radiative transport equation is solved using cubic Bezier
solvers (de la Cruz Rodríguez & Piskunov, 2013). The inversion engine of STiC
includes an equation of state extracted from the SME code (Piskunov & Valenti,
2017).

An important step in the multi-line inversion is needed to properly handle the
strong intensity differences between the various lines. This is particularly
important for the C II 1334 & 1335 Å lines since their intensity is much than
the Mg II h&k, the Mg II UV triplet lines, and the nearby photospheric lines.
Therefore, for the inversion we must scale them with respect to the other
lines. We have also considered this between the Mg II UV2&3 lines with respect
to the Mg II h&k lines, and the photospheric lines with respect to the Mg II
h&k lines. The strategy is as follows: we calculate the most frequent value of
the intensity in a spectroheliogram map in the Mg II h&k $k_{3}$, and we scale
or weight the C II 1334 & 1335 Å lines to that value. That means, we try to
scale the C II 1334 & 1335 Å lines to the same relative intensity level as the
Mg II h&k lines. Note that a somewhat similar strategy was used by Sainz Dalda
& De Pontieu (2022a). However, in that case, because the authors were
investigating the thermodynamics of the chromosphere during the maximum of a
X1.0-class flare, the ratio between the maximum intensity in the Mg II h&k
lines and the C II 1334 & 1335 Å lines varies from $\approx 100$ in the non-
flaring region to just $\approx 10$ in the ribbons, both in the same FoV or
dataset. Therefore, the authors decided to stratify the data considering
different populations of profiles with different intensity. As a consequence,
different weights for the C II 1334 & 1335 Å lines with respect to the Mg II
h&k lines were used for different regions in a data set.The need for this
scaling (or weighting) is due the different behavior of the lines in different
solar features, with the most extreme differences occurring between the non-
flaring and flaring regions. The data sets considered in this study do not
have such a strong variation of intensity between the C II 1334 & 1335 Å and
Mg II h&k lines in different regions within the same field of view. Therefore,
only one weight is used for the C II 1334 & 1335 Å lines for the full data
set. For the Mg II UV1 line, the Mg II UV2&3 lines, and the photospheric lines
we used a slightly different strategy. For these lines we obtained the weights
with respect to the Mg II h&k lines in an empirical fashion, after we verified
that using the strategy of the most frequent value of the Mg II h&k $k_{3}$
map did not provide a significant advantage. In other words, we found that for
these lines, on average, the relative intensity with respect to that of the Mg
II h&k lines for different data sets has a small standard deviation. Thus, for
the sake of simplicity, we have used the same combination of weights (with
respect to Mg II h&k) for these lines. For a dataset, STiC only accepts a
wavelength-dependent noise and weight for all the profiles. For simplicity,
the noise is the same averaged noise (standard deviation in a spectral range)
at all the wavelengths, while the different weights are given for different
spectral positions (usuallly ranges or windows). That means, we use an unique
value of the $w_{i}/\sigma$ per data set in the merit function used to
quantify the quality of the fit between the observed $I(\lambda_{i})^{obs}$,
and the inverted profile $I(\lambda_{i},\mathbf{M})^{syn}$, that is:

$\chi^{2}=\frac{1}{\nu}\sum_{i=0}^{q}{(I(\lambda_{i})^{obs}-I(\lambda_{i},\mathbf{M})^{syn})^{2}\frac{w_{i}^{2}}{\sigma_{i}^{2}}}$
(1)

with $i=0,...,q$ the sampled wavelengths, $w_{i}$ their weights, $\sigma_{i}$
the uncertainties of the observation (e.g. photon noise), and $\nu$ the number
of spectral samples. Note that we assume, for simplicity, that $\sigma_{i}$ is
the same for all $i$, as mentioned above. We have used 3 cycles with different
number of nodes for the physical variables in the model. These values are
detailed in the table 2. We used the FALC quiet-sun model (Fontenla et al.,
1993) as the initial guess model for all the inversions.

No. Cycle | 1 | 2 | 3
---|---|---|---
$T$ | 4 | 7 | 9
$v_{turb}$ | 2 | 4 | 8
$v_{los}$ | 2 | 4 | 8
Table 2: Number of nodes in each cycle for the thermodynamics variables
considered during the inversions.

##  3 Results

Figures 4 and 5 show the inversion of a multi-line RP in the penumbra and the
outer part of the plage of the NOAA 12681, respectively, both on 2017-09-26 at
05:09:50UT. In addition to the observed RP (in fuchsia) and its corresponding
inverted, synthetic RP (best fit obtained by STiC, in black), the temperature
($T$), line-of-sight velocity ($v_{los}$), the velocity of turbulent motions
or micro-turbulence ($v_{turb}$), and the logarithm of the electron density
($n_{e}$) are shown in the last row of these figures. The rest position of the
lines is indicated with a vertical dashed line. The inset plots in the fourth
row show the core of the Mg II h&k and Mg II UV2&3 lines.

As a result of the strategy presented in this article, we have an indexed,
labeled, relational data base. Thus, each inverted synthetic RP has an
associated representative model atmosphere (RMA), and the following meta-data:
the IRIS filename and the solar feature (“umbra”, “penumbra”, “pore-like”,
“plage”, and “quiet-sun”) where the RP was clustered, the position on the Sun
of the original observed profiles associated with the RP’s cluster, the
exposure time, and the observation time. Because information is attached to
the RP-RMA pair in $IRIS^{2+}$ one can search for RMAs for a given solar
feature (e.g. ’plage’), at a given time in the solar cycle (e.g. ’2014-01-01’
$<$ date_obs $<$ ’2014-05-31’), and at given position on the Sun. Figure 6
shows the case of a search of the RMAs for all the solar features located at
$0.65<\mu<0.70$. As one would expect, the RMAs show, for a given physical
parameter, some variation both within the same solar feature (e.g., the
$v_{turb}$ for the RMAs in the plage), and between solar features (e.g umbra
versus penumbra). In this example search, only 2 umbras (30 RMA/umbra) and 2
penumbras (50 RMA/penumbra) were observed and included in the data base. In
fact, these 2 umbras and 2 penumbras observations belong to the same AR, and
for that reason the dispersion between the RMAs for the different physical
variables is rather small in these solar features. Several plages and "quiet-
sun" belonging to different ARs are however found in this search, and for that
reason their associated RMAs show a larger variability. We remind the reader
that the term "quiet-sun" here refers to regions within an active region that
are not plage, umbra, pores, or penumbra. It is not evident that such regions
are, in fact, identical to quiet Sun regions that are not part of an active
region. This is because the overlying large-scale magnetic field topology and
the associated canopy can impact the chromosphere of these regions, which are
at the photospheric level perhaps more quiet-Sun like.

The quality of the inversion is not always easy to evaluate from the
$\chi^{2}$ value in the profiles of $IRIS^{2+}$. This is because of the
behavior of the Euclidean distance when applied to problems with high
dimensions, which is the case of our profiles. In many cases, a bad fit in the
core of the Mg II h&k lines penalizes too much the value of $\chi^{2}$ (see
section 3 of Sainz Dalda & De Pontieu 2022b for a more detailed discussion).
In addition, we are now using multiple lines. That increases the possibility
of having a bad $\chi^{2}$ value when one or a few lines are badly fitted, but
the rest are well fitted. Of course, when all the lines are well fitted, the
$\chi^{2}$ is good (i.e., low). The question is what to do when several lines
are well fitted and others are not (i.e., with a resulting high value for
$\chi^{2}$). As we will discuss later, that is not a big problem, since we
have verified that in a large majority the profiles are well fitted. Even for
the minority of RPs where some lines are badly fitted, there are usually
several other lines that are well fitted.

Figure 4: First row: Spectroheliogram map for Mg II h&k $k_{3}$, the Ti I
2785.46 Å, the photospheric continuum at 2832.04 Å, and a map of the various
types of regions that were identified (see text). The other rows are results
from an inversion for the location marked with a white cross in the “Selected
areas” map in the right of the top row. From the second to the sixth row, from
left to right, the observed (dotted-solid in fuchsia) and the inverted profile
(black) for the C II 1334 & 1335 Å lines, Ti I 2785.46 Å, Mg II UV1 Fe I
2793.22 Å, Mg II h&k lines - including Mg II UV2&3-, Fe I 2809.15 Å, C I
2810.58 Å, Ni I 2815.18 Å, and Fe I 2827.33 Å. The bottom row shows the
temperature ($T$, orange), the logarithm of the electron density ($n_{e}$,
violet), velocity of turbulent motions or micro-turbulent ($v_{turb}$, green)
and the line-of-sight ($v_{los}$, blue). Both the inverted profile and the
representative model atmosphere (RMA) are included in the $IRIS^{2+}$, and
they are both labeled as penumbra in the database. Figure 5: Same as Figure 4
for a multi-line RP and RMA in the outer plage (marked with a white cross in
the “Selected areas" panel). Both the inverted profile and eh representative
model atmosphere (RMA) are included in the $IRIS^{2+}$, and they are both
labeled as ’plage’ in the database. Figure 6: Using $IRIS^{2+}$ data base to
find and visualize the different models of umbra, penumbra, plage and quiet-
sun at a given interval of $\mu=cos(\theta)$, with $\theta$ the heliocentric
angle.

##  4 Discussion

In the example of the $IRIS^{2+}$ data base shown in Figure 4, the inverted
profiles fit the observed profiles very well. Therefore, we can be confident
with the representative model atmosphere (RMA) obtained, and because the fit
is very good for all the lines, the uncertainty of the physical values should
be rather small at all the optical depth ranges shown in the right panel of
Figure 3, i.e., where the RFs are large. What if one or more lines are not
fitting well the observed profile? A very picky reader might consider that is
the case, for instance, for the lines Fe I 2793.22 Å, Fe I 2809.15 Å, C I
2810.58 Å, and Ni I 2815.18 Å. In this case, the reader should consider the
values at the optical depths $-3<$ $log(\tau)$$<-1$ with some caution, since
only Fe I 2827.33 Åis fitting the observed profile very well. Which impact do
the lines that do not fit the observations well have in the $IRIS^{2+}$?
Strictly speaking, none. The $IRIS^{2+}$ database is a relational data base
between model atmospheres and their synthetic profiles. The physics encoded in
the latter are strictly due to the values in the model atmosphere and the
considerations made to solve the radiative transfer equation. Practically
speaking, it has an impact if the fits were consistently bad - which is not
the case.

The $IRIS^{2}$ can be used in different ways: as a set of reference model
atmospheres and profiles for the various types of features within ARs; as
constraints for numerical models; as a look-up table for speeding inversions
problems; or as a source atmosphere to synthesize spectral lines both in the
chromosphere and the photosphere, including full-Stokes polarimetric profiles
by considering magnetic field information. In all these cases, we want to have
a better representation of the thermodynamics in the low solar atmosphere, and
in that sense the better the fits are, the better that representation is.
However, those synthetics profiles that show some poorer fits with respect to
the observed RP very likely still occur in the solar atmosphere. Therefore,
since a very large number of the inversions calculated to build the
$IRIS^{2+}$ data base are showing very good fits, we consider this database
likely the most accurate comprehensive set of stratified-depth models and
profiles in the chromosphere and the photosphere.

##  5 Conclusions

Thanks to the extraordinary multi-line capabilities of IRIS, we have created a
database of 40,320 synthetic representative profiles and their corresponding
representative model atmospheres that convey the essential thermodynamic
information of active regions from the bottom of the photosphere to the top of
the chromosphere.

Thanks to the state-of-the-art multi-line multi-atom STiC inversion code, we
have been able to recover the physical information encoded in 6 chromospheric
lines and 6 photospheric lines. Using the well-known k-means technique we have
clustered multi-line IRIS spectral data. This helps us to overcome the task of
inverting all the profiles considered in this work, which would require very
significant computational resources ($\approx 242,000~{}CPU-hours$).

$IRIS^{2+}$ is an unique data base, obtained systematically from an
unprecedented combination of a large number spectral lines sensitive to the
thermodynamics in the low solar atmosphere. We will continue maximizing the
observational capabilities of IRIS, considering more lines and applying
improved inversion methods and advance machine learning techniques.
$IRIS^{2+}$ has a high potential to be used in different ways to gain
knowledge in the low solar atmosphere. $IRIS^{2+}$, as well as $IRIS^{2}$
already is, will be available to the public in different formats, both for IDL
and Python. We encourage the community to exploit $IRIS^{2+}$ to solve the
physics of the ARs.

IRIS is a NASA small explorer mission developed and operated by LMSAL with
mission operations executed at NASA Ames Research center and major
contributions to downlink communications funded by ESA and the Norwegian Space
Agency. This work was supported by NASA contract NNG09FA40C (IRIS). A. Agrawal
was supported by the Lockheed-Martin-Palo Alto Unified School District
internship program. The authors are grateful to Mats Carlsson and Tiago
Pereira for providing the C atom model, and to Jaime de la Cruz Rodríguez for
his helpful comments on the inversion of photospheric lines. .

## References

  * Bjørgen et al. (2019) Bjørgen, J. P., Leenaarts, J., Rempel, M., et al. 2019, A&A, 631, A33
  * Carlsson et al. (2019) Carlsson, M., De Pontieu, B., & Hansteen, V. H. 2019, ARA&A, 57, 189
  * Cauzzi et al. (2008) Cauzzi, G., Reardon, K. P., Uitenbroek, H., et al. 2008, A&A, 480, 515
  * Collados et al. (1994) Collados, M., Martinez Pillet, V., Ruiz Cobo, B., del Toro Iniesta, J. C., & Vazquez, M. 1994, A&A, 291, 622
  * de la Cruz Rodríguez et al. (2016) de la Cruz Rodríguez, J., Leenaarts, J., & Asensio Ramos, A. 2016, ApJ, 830, L30
  * de la Cruz Rodríguez et al. (2019) de la Cruz Rodríguez, J., Leenaarts, J., Danilovic, S., & Uitenbroek, H. 2019, A&A, 623, A74
  * de la Cruz Rodríguez & Piskunov (2013) de la Cruz Rodríguez, J., & Piskunov, N. 2013, ApJ, 764, 33
  * del Toro Iniesta & Ruiz Cobo (2016) del Toro Iniesta, J. C., & Ruiz Cobo, B. 2016, Living Reviews in Solar Physics, 13, 4
  * del Toro Iniesta et al. (1994) del Toro Iniesta, J. C., Tarbell, T. D., & Ruiz Cobo, B. 1994, ApJ, 436, 400
  * Fontenla et al. (1993) Fontenla, J. M., Avrett, E. H., & Loeser, R. 1993, ApJ, 406, 319
  * Fontenla et al. (2009) Fontenla, J. M., Curdt, W., Haberreiter, M., Harder, J., & Tian, H. 2009, ApJ, 707, 482
  * Gingerich et al. (1971) Gingerich, O., Noyes, R. W., Kalkofen, W., & Cuny, Y. 1971, Sol. Phys., 18, 347
  * Gudiksen et al. (2011) Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154
  * Hale (1908) Hale, G. E. 1908, PASP, 20, 220
  * Khomenko et al. (2014) Khomenko, E., Collados, M., Díaz, A., & Vitas, N. 2014, Physics of Plasmas, 21, 092901
  * Landi Degl’Innocenti & Landi Degl’Innocenti (1981) Landi Degl’Innocenti, E., & Landi Degl’Innocenti, M. 1981, Nuovo Cimento B Serie, 62B, 1
  * Leenaarts et al. (2012) Leenaarts, J., Pereira, T., & Uitenbroek, H. 2012, A&A, 543, A109
  * Leenaarts et al. (2013) Leenaarts, J., Pereira, T. M. D., Carlsson, M., Uitenbroek, H., & De Pontieu, B. 2013, ApJ, 772, 89
  * Martínez-Sykora et al. (2020) Martínez-Sykora, J., Leenaarts, J., De Pontieu, B., et al. 2020, ApJ, 889, 95
  * Mein (1971) Mein, P. 1971, Sol. Phys., 20, 3
  * Pereira et al. (2013) Pereira, T. M. D., Leenaarts, J., De Pontieu, B., Carlsson, M., & Uitenbroek, H. 2013, ApJ, 778, 143
  * Pereira & Uitenbroek (2015) Pereira, T. M. D., & Uitenbroek, H. 2015, A&A, 574, A3
  * Piskunov & Valenti (2017) Piskunov, N., & Valenti, J. A. 2017, A&A, 597, A16
  * Rutten (2007) Rutten, R. J. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 368, The Physics of Chromospheric Plasmas, ed. P. Heinzel, I. Dorotovič, & R. J. Rutten, 27
  * Sainz Dalda et al. (2019) Sainz Dalda, A., de la Cruz Rodríguez, J., De Pontieu, B., & Gošić, M. 2019, ApJ, 875, L18
  * Sainz Dalda & De Pontieu (2022a) Sainz Dalda, A., & De Pontieu, B. 2022a, submitted
  * Sainz Dalda & De Pontieu (2022b) —. 2022b, submitted
  * Sukhorukov & Leenaarts (2017) Sukhorukov, A. V., & Leenaarts, J. 2017, A&A, 597, A46
  * Uitenbroek (2001) Uitenbroek, H. 2001, ApJ, 557, 389
  * Wargnier et al. (2022) Wargnier, Q. M., Martinez-Sykora, J., Hansteen, V. H., & De Pontieu, B. 2022, arXiv e-prints, arXiv:2211.02157
  * Wedemeyer et al. (2016) Wedemeyer, S., Bastian, T., Brajša, R., et al. 2016, Space Sci. Rev., 200, 1
  * Woods et al. (2021) Woods, M. M., Sainz Dalda, A., & De Pontieu, B. 2021, ApJ, 922, 137


# Ultrathin Al1-xScxN for low-voltage driven ferroelectric-based devices

Georg Schönweger, , Md Redwanul Islam, , Niklas Wolff, Adrian Petraru, Lorenz
Kienle, Hermann Kohlstedt, Simon Fichtner

###### Abstract

Thickness scaling of ferroelectricity in Al1-xScxN is a determining factor for
its potential application in neuromorphic computing and memory devices. In
this letter, we report on ultrathin (10 nm) Al0.72Sc0.28N films that are
ferroelectrically switchable at room temperature. All-epitaxial
Al0.72Sc0.28N/Pt heterostructures are grown by magnetron sputtering onto
GaN/sapphire substrates followed by an in situ Pt capping approach to avoid
oxidation of the Al0.72Sc0.28N film surface. Structural characterization by
X-ray diffraction and transmission electron microscopy reveals the established
epitaxy. The thus obtained high-quality interfaces in combination with the in
situ capping is expected to facilitate ferroelectric switching of Al1-xScxN in
the ultrathin regime. The analysis of the relative permittivity and coercive
field dependence on the Al0.72Sc0.28N film thicknesses in the range of 100 nm
down to 10 nm indicates only moderate scaling effects, suggesting that the
critical thickness for ferroelectricity is not yet approached. Furthermore,
the deposited layer stack demonstrates the possibility of including ultrathin
ferroelectric Al1-xScxN into all-epitaxial GaN-based devices using sputter
deposition techniques. Thus, our work highlights the integration and scaling
potential of all-epitaxial ultrathin Al1-xScxN offering high storage density
paired with low voltage operation desired for state of the art ferroelectric
memory devices.

###### Index Terms:

aluminum-scandium-nitride (Al1-xScxN), ultrathin, ferroelectric, gallium
nitride, epitaxial growth

##  I Introduction

The high coercive field ($E_{c}$) and the high, stable remanent polarization
separates the ferroelectric properties recently discovered in materials with
wurtzite-type structure from classical ferroelectrics[1, 2, 3, 4]. This raises
hopes for particularly good scalability of wurtzite-type based ferroelectric
devices. In addition, the CMOS compatibility and the well established
industrial deposition process make Al1-xScxN thin films highly attractive for
building novel neuromorphic computing and memory devices such as ferroelectric
field-effect transistors (FeFET) and ferroelectric tunnel junctions (FTJ)[5,
6, 7, 8, 9]. Furthermore, it is expected that wurtzite-type ferroelectrics
such as Al1-xScxN introduce ferroelectricity into III-N technology, resulting
in a straightforward approach to realize embedded memory. Recently,
ferroelectric all-epitaxial all-wurtzite type Al1-xScxN/GaN heterostructures
were demonstrated[10, 11]. However, a very low film thickness of the
ferroelectric layer is needed for following the general trend of
miniaturization and increasing storage density in all of the aforementioned
devices. In this context, Al1-xScxN offers high scalability due to its high
$E_{c}$, making it possible to tailor the film thickness to the ultrathin
regime to achieve reasonable memory windows and low operating voltages[12].
Furthermore, in terms of device design, ultrathin ferroelectric films are a
prerequisite for building FTJs[7]. However, reducing the thickness to the
ultrathin regime ($<$ 30 nm) is often accompanied with a material-specific
diminution of the remanent polarization or results in a total loss of
ferroelectricity[13]. Up to now, only a small number of studies on the
thickness scaling properties of ferroelectric Al1-xScxN were conducted,
revealing ferroelectricity for $\approx$20 nm film thickness at room
temperature as well as for $\approx$10 nm film thickness at elevated
temperatures (373 K)[14, 15, 16]. Here we present our recent results on 10 nm
thick epitaxial Al0.72Sc0.28N grown on epitaxial Pt/GaN templates using
sputter deposition. Typical butterfly-shaped capacitance-voltage ($C-U$) loops
were recorded at room temperature, demonstrating clearly distinguishable
ferroelectric switching. Thus, the general feasibility of ultrathin Al1-xScxN
for future device integration and low voltage operation is demonstrated.

##  II Device fabrication and characterization methods

All films were grown by sputter deposition in an Oerlikon (now Evatec) MSQ 200
multisource system on top of commercially available GaN/sapphire substrates. A
12 nm thick epitaxial Pt layer (bottom electrode) was DC-sputtered at 500 °C,
600 W, and 50 sccm Kr flow. XRD measurements were performed using
monochromatic K$\alpha$ radiation in a Seifert XRD 3000 PTS system
($\theta$-2$\theta$ scan) as well as in a Rigaku SmartLab diffractometer (Phi-
scan). The Al0.72Sc0.28N thin films were grown by pulsed-DC co-sputtering at
450 °C, details about the process can be found elsewhere[17, 18]. The 10 nm
thick Al0.72Sc0.28N layer was capped in situ with 100 nm thick Pt. Square top-
electrodes were structured via lithography and ion beam etching (IBE, Oxford
Instruments Ionfab 300). Capacitance measurements were performed on a Hewlett
Packard 4284A Precision LCR meter. A cross-section specimen of the in situ
capped 10 nm thick Al0.72Sc0.28N sample has been prepared using a standard
focused ion beam (FIB) procedure and analyzed using high-resolution
transmission electron microscopy (HRTEM, Tecnai F30, operated at 300 kV, field
emission gun) and Fast Fourier Transform (FFT) patterns.

##  III Results and Discussions

The voltage dependence of the relative permittivity ($\varepsilon_{r}$) is
depicted in Fig. 1. The butterfly-shaped curve (black) demonstrates
distinguishable ferroelectric switching of the 10 nm thick Al0.72Sc0.28N film.

Figure 1: Relative permittivity and capacitance as a function of the applied
voltage and electric field for Pt(100 nm)/Al0.72Sc0.28N(10 nm)/Pt(12 nm)
capacitors deposited on a GaN/sapphire template. The non switching half-loops
(red) were recorded by previously switching the capacitor to the respective
polarity. Measured on 10 x 10 µm² pads with a small signal of 1 Vpp at 100¡
kHz swept in 0.1 V steps with a delay of 300 ms per step.

Additionally non switching half-loops (red) were recorded by measuring twice
in the same direction. No hysteretic behaviour is observed for the non-
switching half-loops, thus giving additional evidence that the butterfly-
shaped hysteresis for the switching loop originates from ferroelectric
switching. Previously, the coercive field in Al1-xScxN thin films was observed
to be almost independent of film thickness down to 27 nm[14]. Similarly, the
$E_{c}$ fields determined via $C$-$U$ loops in this study, as depicted in Fig.
2, is not changing considerably when decreasing the film thickness from 100 nm
to 50 nm, while at 20 nm and escpecially at 10 nm film thickness, a moderate
increase of the coercive field by about 20% is observed.

Figure 2: $\varepsilon_{r}$ and $E_{c}$ in dependence of Al0.72Sc0.28N film
thickness. $\varepsilon_{r}$ was measured at 100 kHz, 100 mV at 0 V bias.
$E_{c}$ was determined via measuring $C-U$ loops keeping the time for a full
sweep constant for the various film thicknesses as well as via $J-E$ using a
sinus signal at 80 kHz. The error bars were calculated using an estimated
capacitor side-length error of 0.4 µm, an estimated thickness error of 5 % and
an estimated error for determination of the $E_{c}$. Additionally for $E_{c}$
determination via $C-U$ loops the step width of the respective voltage sweep
was included in the error calculation.

This is in strong contrast to classical ferroelectrics, e.g. the perovskites
(PZT, BTO, etc.), were a pronounced increase of $E_{c}$ with decreasing film
thickness is observed[19, 13]. Furthermore, the often reported degradation of
ferroelectric properties due to a non-switching interfacial layer (dead
layer[20]) in ultrathin films is typically accompanied by a decrease of
$\varepsilon_{r}$[19, 21]. Thus, ferroelectricity in 10 nm thick Al1-xScxN
appears to be not yet influenced considerably by such interfacial effects. The
non-scaling of $E_{c}$ and $\varepsilon_{r}$ with thickness leads to our
conclusion that a critical thickness in Al1-xScxN is not yet approached at 10
nm and stable ferroelectricity can be expected also for films $<$ 10 nm
thickness.

The $\theta$-2$\theta$ scan depicted in Fig. 3 (left) reveals the expected
out-of-plane 111 orientation for Pt and 0002 orientation for GaN. Laue
oscillations, which are sensitive to crystalline disorder, are appearing for
the Pt 111 reflection[22].

Figure 3: a) $\theta$-2$\theta$ scan and b) smoothed phi-scan for the non-
specular Pt 131 reflection (right) of 12 nm Pt on GaN/sapphire . c)
$\theta$-2$\theta$ of the film shown in a) after deposition of 10 nm
Al0.72Sc0.28N followed by 5 nm in-situ capping with Pt.

The thus expected high crystallinity is confirmed by rocking curve (RC)
measurements (not shown here) of the 111 Pt reflection; the full width at half
maximum (FWHM) of 0.1 ° is indicative of a strain and defect-poor crystalline
phase with high out-of-plane orientation. Furthermore, the distinct Laue
oscillations allow to determine the thickness ($t$) of the ordered crystal
volume[22].

$t=\frac{3\lambda}{2\sin{\theta_{1}}-\sin{\theta_{-1}}}$ (1)

Using equation 1, with $\theta_{1}$ and $\theta_{-1}$ corresponding to the
maximum of the peaks right and left from the center Pt 111 reflection, an
average film thickness of 11 nm is determined.

Figure 4: HRTEM image (right) and corresponding averaged thickness profile
(left) of the Pt(100 nm)/Al0.72Sc0.28N(10 nm)/Pt(12 nm)/GaN/sapphire
heterostructure. In the inset, the FFT pattern for Al0.72Sc0.28N (i) and Pt
(ii) is depicted.

The in-plane ordering was resolved by recording a phi-scan of the non-specular
131 Pt reflection, depicted in Fig. 3 (right). A 6-fold symmetry is visible,
evidencing epitaxial (in-plane ordered) growth.

The intensity profile across the HRTEM micrograph is shown in Fig. 4 and film
thicknesses of 10 nm for Al0.72Sc0.28N and 12 nm for Pt are determined on
local average. The latter is in excellent agreement with the Pt thickness
calculated using equation 1. The FFT pattern clearly reveals an epitaxial
growth for both layers. The epitaxial growth in combination with the in situ
capping of Al1-xScxN to avoid oxidation is expected to have a significant
impact on the ferroelectric switching behaviour, especially for ultrathin
films. The blurred interfaces present in the HRTEM image are attributed to the
finite roughness of the heterostructure.

##  IV Conclusion

Distinguishable ferroelectric switching of 10 nm thick Al0.72Sc0.28N was
demonstrated at room temperature. Largely constant coercive fields and
constant $\varepsilon_{r}$ for films ranging from 100 nm down to 10 nm
thickness were confirmed. Thus, no indications of instabilities of the
ferroelectric phase or considerable contributions from dead layers are
present. Therefore, it is expected that further potential remains to scale
ferroelectric Al1-xScxN to even lower thicknesses. These results highlight the
general feasibility of using Al1-xScxN as an active layer for integrated
electronic devices with low operating voltage and high storage density.

## Acknowledgment

This work was supported by the project “ForMikro-SALSA” (grant no. 16ES1053)
from the Federal Ministry of Education and Research (BMBF) and the Deutsche
Forschungsgemeinschaft (DFG) under the scheme of the collaborative research
centers (CRC) 1261 and 1461.

## References

  * [1] S. Fichtner, N. Wolff, F. Lofink, L. Kienle, and B. Wagner, “AlScN: A III-v semiconductor based ferroelectric,” _Journal of Applied Physics_ , vol. 125, no. 11, p. 114103, mar 2019.
  * [2] K. Ferri, S. Bachu, W. Zhu, M. Imperatore, J. Hayden, N. Alem, N. Giebink, S. Trolier-McKinstry, and J.-P. Maria, “Ferroelectrics everywhere: Ferroelectricity in magnesium substituted zinc oxide thin films,” _Journal of Applied Physics_ , vol. 130, no. 4, p. 044101, jul 2021.
  * [3] J. Hayden, M. D. Hossain, Y. Xiong, K. Ferri, W. Zhu, M. V. Imperatore, N. Giebink, S. Trolier-McKinstry, I. Dabo, and J.-P. Maria, “Ferroelectricity in boron-substituted aluminum nitride thin films,” _Physical Review Materials_ , vol. 5, no. 4, p. 044412, apr 2021.
  * [4] D. Wang, P. Wang, B. Wang, and Z. Mi, “Fully epitaxial ferroelectric ScGaN grown on GaN by molecular beam epitaxy,” _Applied Physics Letters_ , vol. 119, no. 11, p. 111902, sep 2021.
  * [5] H. Mulaosmanovic, T. Mikolajick, and S. Slesazeck, “FeFETs for neuromorphic systems,” in _Topics in Applied Physics_. Springer Singapore, 2020, pp. 399–411.
  * [6] T. Schenk, M. Pešić, S. Slesazeck, U. Schroeder, and T. Mikolajick, “Memory technology—a primer for material scientists,” _Reports on Progress in Physics_ , vol. 83, no. 8, p. 086501, jul 2020.
  * [7] E. Y. Tsymbal and H. Kohlstedt, “Tunneling across a ferroelectric,” _Science_ , vol. 313, no. 5784, pp. 181–183, jul 2006.
  * [8] H. Kohlstedt, N. A. Pertsev, J. R. Contreras, and R. Waser, “Theoretical current-voltage characteristics of ferroelectric tunnel junctions,” _Physical Review B_ , vol. 72, no. 12, p. 125341, sep 2005.
  * [9] X. Liu, D. Wang, J. Zheng, P. Musavigharavi, J. Miao, E. A. Stach, R. H. O. III, and D. Jariwala, “Post-cmos compatible aluminum scandium nitride/2d channel ferroelectric field-effect-transistor,” 2020.
  * [10] P. Wang, D. Wang, N. M. Vu, T. Chiang, J. T. Heron, and Z. Mi, “Fully epitaxial ferroelectric ScAlN grown by molecular beam epitaxy,” _Applied Physics Letters_ , vol. 118, no. 22, p. 223504, may 2021.
  * [11] G. Schönweger, A. Petraru, M. R. Islam, N. Wolff, B. Haas, A. Hammud, C. Koch, L. Kienle, H. Kohlstedt, and S. Fichtner, “From fully strained to relaxed: Epitaxial ferroelectric al 1- x sc x n for III-n technology,” _Advanced Functional Materials_ , p. 2109632, feb 2022.
  * [12] H. Mulaosmanovic, E. T. Breyer, T. Mikolajick, and S. Slesazeck, “Ferroelectric FETs with 20-nm-thick HfO2layer for large memory window and high performance,” _IEEE Transactions on Electron Devices_ , vol. 66, no. 9, pp. 3828–3833, sep 2019.
  * [13] H. Qiao, C. Wang, W. S. Choi, M. H. Park, and Y. Kim, “Ultra-thin ferroelectrics,” _Materials Science and Engineering: R: Reports_ , vol. 145, p. 100622, jul 2021.
  * [14] S. Fichtner, G. Schönweger, T.-N. Kreutzer, A. Petraru, H. Kohlstedt, F. Lofink, and B. Wagner, “Ferroelectricity in alscn: Switching, imprint and sub-150 nm films,” in _Proc. 2020 IEEE ISAF_ , 2020.
  * [15] D. Wang, J. Zheng, P. Musavigharavi, W. Zhu, A. C. Foucher, S. E. Trolier-McKinstry, E. A. Stach, and R. H. Olsson, “Ferroelectric switching in sub-20 nm aluminum scandium nitride thin films,” _IEEE Electron Device Letters_ , vol. 41, no. 12, pp. 1774–1777, dec 2020.
  * [16] R. Mizutani, S. Yasuoka, T. Shiraishi, T. Shimizu, M. Uehara, H. Yamada, M. Akiyama, O. Sakata, and H. Funakubo, “Thickness scaling of al0.8sc0.2)n films with remanent polarization beyond 100 $\mu$c cm-2 around 10 nm in thickness,” _Applied Physics Express_ , vol. 14, no. 10, p. 105501, sep 2021.
  * [17] S. Fichtner, T. Reimer, S. Chemnitz, F. Lofink, and B. Wagner, “Stress controlled pulsed direct current co-sputtered al1-xScxN as piezoelectric phase for micromechanical sensor applications,” _APL Materials_ , vol. 3, no. 11, p. 116102, nov 2015.
  * [18] S. Fichtner, N. Wolff, G. Krishnamurthy, A. Petraru, S. Bohse, F. Lofink, S. Chemnitz, H. Kohlstedt, L. Kienle, and B. Wagner, “Identifying and overcoming the interface originating c-axis instability in highly sc enhanced AlN for piezoelectric micro-electromechanical systems,” _Journal of Applied Physics_ , vol. 122, no. 3, p. 035301, jul 2017.
  * [19] N. A. Pertsev, J. R. Contreras, V. G. Kukhar, B. Hermanns, H. Kohlstedt, and R. Waser, “Coercive field of ultrathin pb(zr0.52ti0.48)o3 epitaxial films,” _Applied Physics Letters_ , vol. 83, no. 16, pp. 3356–3358, oct 2003.
  * [20] C. Zhou and D. M. Newns, “Intrinsic dead layer effect and the performance of ferroelectric thin film capacitors,” _Journal of Applied Physics_ , vol. 82, no. 6, pp. 3081–3088, sep 1997.
  * [21] M. Stengel and N. A. Spaldin, “Origin of the dielectric dead layer in nanoscale capacitors,” _Nature_ , vol. 443, no. 7112, pp. 679–682, oct 2006\.
  * [22] X. Castel, M. Guilloux-Viry, A. Perrin, J. Lesueur, and F. Lalu, “High crystalline quality CeO2 buffer layers epitaxied on sapphire for YBa2cu3o7 thin films,” _Journal of Crystal Growth_ , vol. 187, no. 2, pp. 211–220, may 1998.