11institutetext: Rosseland Centre for Solar Physics, University of Oslo, PO
Box 1029 Blindern, 0315 Oslo, Norway
11email<EMAIL_ADDRESS>22institutetext: Institute of Theoretical
Astrophysics, University of Oslo, PO Box 1029, Blindern 0315, Oslo, Norway

# Formation of H $\varepsilon$ in the solar atmosphere

K. Krikova 1122 T. M. D. Pereira 1122 L. H. M. Rouppe van der Voort 1122

###### Abstract

Context. In the solar spectrum, the Balmer series line H $\varepsilon$ is a
weak blend on the wing of Ca ii H. Recent high-resolution H $\varepsilon$
spectroheliograms reveal a reversed granulation pattern and in some cases,
even unique structures. It is apparent that H $\varepsilon$ could potentially
be a useful diagnostic tool for the lower solar atmosphere.

Aims. Our aim is to understand how H $\varepsilon$ is formed in the quiet Sun.
In particular, we consider the particular physical mechanism that sets its
source function and extinction, how it is formed in different solar
structures, and why it is sometimes observed in emission.

Methods. We used a 3D radiative magnetohydrodynamic (MHD) simulation that
accounts for non-equilibrium hydrogen ionization, run with the _Bifrost_ code.
To synthesize H $\varepsilon$ and Ca ii H spectra, we made use of the RH code,
which was modified to take into account the non-equilibrium hydrogen
ionization. To determine the dominant terms in the H $\varepsilon$ source
function, we adopted a multi-level description of the source function. Making
use of the synthetic spectra and simulation, we studied the contribution
function to the relative line absorption or emission and compared it with
atmospheric quantities at different locations.

Results. Our multi-level source function description suggests that H
$\varepsilon$ source function is dominated by interlocking, with the dominant
interlocking transition being through the ground level, populating the upper
level of H $\varepsilon$ via the Lyman series. This makes the H $\varepsilon$
source function partly sensitive to temperature. The H $\varepsilon$
extinction is set by Lyman- $\alpha$. In some cases, this temperature
dependence gives rise to H $\varepsilon$ emission, indicating heating. The
typical absorption profiles show reversed granulation and the H $\varepsilon$
line core reflects mostly the Ca ii H background radiation.

Conclusions. Synthetic H $\varepsilon$ spectra can reproduce quiet Sun
observations quite well. High-resolution observations reveal that H
$\varepsilon$ is not just a weak absorption line. Regions with H $\varepsilon$
in emission are especially interesting to detect small-scale heating events in
the lower solar atmosphere, such as Ellerman bombs. Thus, H $\varepsilon$ can
be an important new diagnostic tool for studies of heating in the solar
atmosphere, augmenting the diagnostic potential of Ca ii H when observed
simultaneously.

###### Key Words.:

Radiative transfer – Line: formation – Sun: photosphere – Sun: chromosphere

##  1 Introduction

The H $\varepsilon$ spectral line is one of the lesser known members of the
Balmer hydrogen series associated with the transition from energy level $n=7$
to $n=2$, at $\lambda_{0}=397.1202$ nm in vacuum. In the solar spectrum, H
$\varepsilon$ is a weak feature blended on the red wing of the Ca ii H
resonance line, one of the strongest lines in the solar spectrum.
Understanding the formation of H $\varepsilon$ holds an intrinsic value, but
it is even more relevant in the context of simultaneous observations with Ca
ii H, a widely used line to study the solar chromosphere. Combining both lines
has the potential to offer additional information about the dynamics of the
solar atmosphere.

Overall, H $\varepsilon$ is rarely studied in the solar context, but it has
received more attention in other stars, where it is used as an indicator of
chromospheric activity (Montes et al., 1995, 1996). Only a handful of studies
have focused on modeling (Ayres & Linsky, 1975; Bayazitov & Sakhibullin, 1991)
and observing H $\varepsilon$ (Evershed, 1930; Turova, 1994) in the Sun.
Observations of H $\varepsilon$ are either connected to studies of solar
flares (Rolli & Magun, 1995; Rolli et al., 1998a, b) or prominences (Nikaidou,
1982; Anan et al., 2017; Zapiór et al., 2022). What all of these studies have
in common is that they highlight cases where H $\varepsilon$ is seen in
emission relative to the Ca ii H background. The most archetypical case of H
$\varepsilon$ emission can be found in the K giant Arcturus (Wilson, 1938;
Wellmann, 1940; Popper, 1956; Ayres & Linsky, 1975) or late-type stars such as
K and M giants (Wilson, 1957).

The work of Ayres & Linsky (1975) is one of the first studies of the formation
of H $\varepsilon$ in the solar atmosphere, but it was limited to 1D plane-
parallel atmospheres and disk-averaged solar observations. One of the main
conclusions drawn from this study is that the H $\varepsilon$ source function
is dominated by the Balmer continuum radiation field, both in the Sun and
Arcturus. In the case of Arcturus, the authors argue that H $\varepsilon$
appears in emission due to a different column mass structure and a stronger
sensitivity to the chromospheric temperature increase. As we show later in the
present work, high spatial resolution solar observations taken at the H
$\varepsilon$ line core show a reversed granulation pattern. If H
$\varepsilon$ is formed under optically-thick conditions, this is difficult to
reconcile with its source function being dominated by the Balmer continuum,
since, in that case, we would expect to observe photospheric granulation in
the H $\varepsilon$ line core. This puzzling aspect has prompted us to revisit
the formation of H $\varepsilon$.

Additional motivation for this work came from a recent study by Joshi et al.
(2020), who show a much higher occurrence of quiet Sun Ellerman bombs (QSEBs)
in H $\beta$ than in H $\alpha$ observations. Joshi et al. (2020) raises the
question of whether H $\varepsilon$ could be an even better diagnostic tool
for detecting QSEBs. Because of its shorter wavelength, H $\varepsilon$ has
the potential to probe photospheric reconnection at even higher spatial
resolution and enhanced intensity contrast, and it is less affected by
extinction from the overlying chromospheric canopy fibrils.

In this study, we want to understand the formation of H $\varepsilon$ in the
dynamic solar atmosphere, using a state-of-the-art 3D radiative magneto-
hydrodynamic (MHD) simulations in combination with non-LTE radiative transfer.
The outline of the paper is as follows. In Sect. 2 we describe recent
observations of the Ca ii H and H $\varepsilon$ lines. Section 3 details our
forward-modeling methods aimed at understanding the formation of H
$\varepsilon$ in a dynamic atmosphere. Our results are presented in Section 4,
where we detail the formation of H $\varepsilon$ for different structures seen
in observations and investigate the influence of Balmer continuum line
blanketing, hydrogen non-equilibrium ionization, and 3D radiative transfer on
the formation of H $\varepsilon$. We present our discussion in Sect. 5,
followed by our concluding remarks in Sect. 6.

##  2 Observations

To compare our synthetic spectra with solar observations, we made use of data
taken with the CHROMIS instrument installed at the Swedish 1-m Solar Telescope
(SST, Scharmer et al., 2003) located at La Palma. CHROMIS is a Fabry-Pérot
filtergraph based on the design by Scharmer (2006) that is capable of fast
wavelength sampling. We observed the Ca ii H line together with H
$\varepsilon$ at 47 wavelength positions with a cadence of 13 s.

CHROMIS has a narrowband spectral transmission profile width of $0.012$ nm. We
covered a spectral range from $396.63$ nm to $397.09$ nm. Spectral sampling
was performed with a critical sampling of $0.006$ nm steps within $\pm 0.072$
nm around the Ca ii H line core and between an offset of $+0.132$ to $+0.168$
nm relative to the Ca ii H line core, covering the H $\varepsilon$ line. The H
$\varepsilon$ line is located at $+0.156$ nm redward of the Ca ii H line core.
The sampling positions outside of the ranges mentioned above were taken with
$0.012$ nm steps. An additional wavelength position at $400$ nm was used to
sample the continuum. CHROMIS includes a wide band (WB) channel that is used
for image restoration. The WB imaging filter has a center wavelength of 395
nm, a full width at half maximum (FWHM) of 1.32 nm, and displays a
photospheric scene.

CHROMIS data were processed using the SSTRED data processing pipeline,
(Löfdahl et al., 2021), which includes multi-object multi-frame blind
deconvolution (MOMFBD, van Noort et al., 2005) image restoration. Here, we
analyzed spectral scans that were recorded under excellent seeing conditions
and the spatial resolution is close to the diffraction limit of the telescope
($\lambda/D=0\aas@@fstack{\prime\prime}08$ at the wavelength of Ca ii H for
the $D=0.97$ m clear aperture of the SST). The SST adaptive optics wavefront
sensor measured the seeing quality (see Scharmer et al., 2019). The Fried’s
parameter $r_{0}$ for the ground-layer seeing peaked above 50 cm during the
best seeing moments. CHROMIS has a pixel size of
$0\aas@@fstack{\prime\prime}0375$.

The targets of our CHROMIS observations were a quiet Sun and a pore region.
The quiet Sun observation was taken on June 22, 2021 at 15:28 UT near the disk
center at (x, y) = ($12\arcsec$, $14\arcsec$). The pore observation was taken
on June 28, 2022 starting at 08:30 UT close to AR 13040 at (x, y) =
($232\arcsec$, $-220\arcsec$) and is part of a longer time series. Two
snapshots of this series are shown in this work. The first snapshot was taken
around 08:49:45 UT and the second around 09:17:52 UT. For the purposes of our
work, we only show a small part of the total field of view with a spatial
extent of $\approx$24 Mm and $\approx$17 Mm (for the 09:17:52 UT snapshot).

##  3 Methods

###  3.1 Model atmospheres

We used two different model atmospheres in our study. First, the semi-
empirical model C from Fontenla et al. (1993) which we will refer to as the
FALC model. The 1D time-independent FALC atmosphere was constructed to match
temporal and spatial averaged spectra in the UV. This “simple” atmosphere
gives a good starting point to understand the formation of H $\varepsilon$ and
compare synthetic profiles to temporal and spatial averaged observations of
the Sun, such as observations from the Fourier transform spectrometer (FTS)
atlas (Kurucz et al., 1984).

The second model atmosphere we use is the publicly available Bifrost model
from Carlsson et al. (2016). Bifrost is a 3D radiation magneto-hydrodynamic
simulation code that solves the resistive MHD equation on a Cartesian grid
including various important physical effects in the solar atmosphere (for more
details, see Gudiksen et al., 2011). One important physical ingredient is
hydrogen non-equilibrium ionization (HNEI), influencing the hydrogen level
populations and, therefore, the formation of H $\varepsilon$. The Cartesian
grid of the numerical model consists of $504\times 504\times 496$ grid points,
spanning a computational box of $24\times 24\times 14.4$ Mm. The simulation
consists of multiple snapshots in time. We used snapshot 385 at simulation
time $t=3850$ s, which is the first published snapshot with non-equilibrium
hydrogen ionization.

###  3.2 Model atoms

The H $\varepsilon$ line is a weak line blended in the strong Ca ii H wing.
Therefore, we have to properly model the Ca ii H line with partial frequency
redistribution (PRD) over the line profile as the inner wings of Ca ii H are
affected by PRD at the H $\varepsilon$ wavelength. To model Ca ii H, we used
the five-level plus continuum Ca model atom from Carlsson & Leenaarts (2012).
We modeled the infrared triplet of Ca II (at $849.8$, $854.2$, and $866.2$ nm)
with the assumption of complete frequency redistribution (CRD), as well as the
Ca II K line. The resonance doublet and infrared triplet of Ca II are affected
by the cross-redistribution (XRD) of photons between these two spectral lines,
generally known as Raman scattering. Therefore, XRD will affect the inner
wings of the Ca ii H and K line mapping heights where the Ca triplet lines are
formed. These give a different escape route of Ca ii H and K photons and will
decrease the inner wing intensity of the lines. The influence of XRD effects
in Ca ii was studied in some detail by Bjørgen et al. (2018). We tested the
influence of XRD on the Ca ii H wing intensity at the wavelength where H
$\varepsilon$ is located and found minor intensity changes if we treat the
doublet and triplet lines in XRD. Consequently, we only model Ca ii H in PRD
and treat all the other lines in CRD to decrease the computational time. We
used the approximation of angle-dependent PRD of Leenaarts et al. (2012b).

To construct an efficient hydrogen model atom to synthesize H $\varepsilon$,
we started with a nineteen-level plus continuum model atom and collapsed
collectively levels starting from the uppermost hydrogen energy level down to
the upper energy level of H $\varepsilon$ ($n=7$). The goal was to find the
smallest hydrogen model atom, which is computationally the most efficient and
models the H $\varepsilon$ line best, as compared with the nineteen-level
hydrogen atom. We found that an eight-level plus continuum hydrogen atom,
where the eighth-level consists of the collapsed levels from $n=19$ to $n=8$,
provides a good approximation of the H $\varepsilon$ line when compared with
the nineteen-level atom. Other collapsed model atoms seem to model the H
$\varepsilon$ line less accurately (although, generally, the difference is
small), by representing the recombination cascades of hydrogen less correctly.
We tested if the assumption of PRD for other hydrogen transitions has an
influence on the emergent H $\varepsilon$ line profile. Tests show that the
influence of PRD in H $\varepsilon$ is negligible and we therefore treated all
hydrogen transitions in CRD. To summarize, we chose to model H $\varepsilon$
with a collapsed eight-level plus continuum atom where H $\varepsilon$ is
calculated in CRD as well as all other hydrogen transitions.

###  3.3 Total source function

The total source function, $S^{\mathrm{t}}_{\nu}$, at the wavelength of H
$\varepsilon$ can be written as an expression of the line source function,
$S^{\mathrm{l}}_{\nu}$, and background source function,
$S^{\mathrm{b}}_{\nu}$,

$S^{\mathrm{t}}_{\nu}=\delta_{\nu}\,S^{\mathrm{l}}_{\nu}+\left(1-\delta_{\nu}\right)\,S^{\mathrm{b}}_{\nu},$
(1)

where

$\delta_{\nu}=\frac{\chi^{\mathrm{l}}_{\nu}}{\left(\chi^{\mathrm{l}}_{\nu}+\chi^{\mathrm{b}}_{\nu}\right)},$
(2)

which is the relative line extinction. The relative line extinction gives the
probability per total extinction that a photon is destroyed by a line process.
The background source function, $S^{\mathrm{b}}_{\nu}$, contains all
background emissivity and extinction contributions (Ca ii H, continuum and
others) at the H $\varepsilon$ wavelength not related to the H $\varepsilon$
transition.

###  3.4 Radiative transfer

We synthesized the Ca II H line with overlapping H $\varepsilon$ wavelength
points using the RH code (Uitenbroek, 2001) and the 1.5D version (Pereira &
Uitenbroek, 2015). This code solves the non-LTE radiative transfer problem
with overlapping wavelength points for bound-bound and bound-free transitions.
Forward modeling with the RH 1.5D code will treat every vertical column of a
3D model atmosphere as an independent plane-parallel 1D atmosphere (1.5D
approximation) with no horizontal transfer of radiation. Because of this
shortcoming, we went on to use the 3D branch of RH to synthesize Ca ii H and H
$\varepsilon$ to check for 3D effects in a few regions.

To study the formation of the H $\varepsilon$ spectral line, we need to
distinguish between the region of formation of the emergent intensity; in our
case, this is the Ca ii H wing intensity and the region where the H
$\varepsilon$ line depression is formed. As we are dealing with a “weak” line
blend in a strong background continuum, the Ca ii H line, the use of
contribution function (CF) to intensity can be misleading (Magain, 1986) and
can give a misleading impression of where the line depression is formed. In
the case of faint or weak spectral lines, one has to distinguish between the
origin of the continuum and the origin of the region where the line depression
is formed. Otherwise, this could lead to a wrong interpretation where the
emergent intensity of a spectral line is formed (de Jager, 1952; Gurtovenko et
al., 1974). Magain (1986) gives a derivation of the relative contribution
function to intensity (line depression or line emission):

$R_{\nu}=\frac{\left(I^{\mathrm{b}}_{\nu}-I^{\mathrm{l}}_{\nu}\right)}{I^{\mathrm{b}}_{\nu}},$
(3)

and expresses the transfer equation for $R_{\nu}$, identifying the integrand
of the formal solution as the CF to relative line depression (for a detailed
derivation, see, Magain, 1986). Also, $I^{\mathrm{b}}_{\nu}$ is the background
(continuum) intensity and $I^{\mathrm{l}}_{\nu}$ is the intensity of the
spectral line. Thus, the formal solution and relative contribution function
for emergent line depression can be expressed as:

$\displaystyle
R_{\nu}=\int_{z_{0}}^{z_{1}}S^{\mathrm{R}}_{\nu}\,e^{\tau^{\mathrm{R}}_{\nu}}\,\chi^{\mathrm{R}}_{\nu}\,\textrm{d}z,$
(4) $\displaystyle
C_{\mathrm{R}}(\nu,z)=\frac{\textrm{d}I_{\mathrm{R}}}{dz}=S^{\mathrm{R}}_{\nu}\,e^{-\tau^{\mathrm{R}}_{\nu}}\,\chi^{\mathrm{R}}_{\nu},$
(5)

with $S^{\mathrm{R}}_{\nu}$ and $\chi^{\mathrm{R}}_{\nu}$ given as:

$\displaystyle
S^{\mathrm{R}}_{\nu}=\left(1-\frac{S^{\mathrm{l}}_{\nu}}{I^{\mathrm{b}}_{\nu}}\right)/\left(1+\frac{\chi^{\mathrm{b}}_{\nu}}{\chi^{\mathrm{l}}_{\nu}}\frac{S^{\mathrm{b}}_{\nu}}{I^{\mathrm{b}}_{\nu}}\right),$
(6)
$\displaystyle\chi^{\mathrm{R}}_{\nu}=\chi^{\mathrm{l}}_{\nu}+\chi^{\mathrm{b}}_{\nu}\,\frac{S^{\mathrm{b}}_{\nu}}{I^{\mathrm{b}}_{\nu}}.$
(7)

The relative contribution function can be simplified to

$C_{\mathrm{R}}(\nu,z)=\frac{\textrm{d}I_{\mathrm{R}}}{dz}=\chi^{\mathrm{l}}_{\nu}\,\left(1-\frac{S^{\mathrm{l}}_{\nu}}{I^{\mathrm{b}}_{\nu}}\right)\,e^{-\tau^{\mathrm{R}}_{\nu}},$
(8)

which gives a more intuitive way to understand relative spectral line
formation. We see a contribution to line depression or emission if
$\chi^{\mathrm{l}}_{\nu}\neq 0$ or $S^{\mathrm{l}}_{\nu}\neq
I^{\mathrm{b}}_{\nu}$. The former says that some absorber has to be present
for a spectral line to form and the latter tells us that the photons created
by the line processes should not be equal to the photons created by the
background processes at the region where the line is formed. Further, the
expression in the brackets of Eq. (8) defines whether there is a contribution
to relative line emission $S^{\mathrm{l}}_{\nu}>I^{\mathrm{b}}_{\nu}$ or
depression $S^{\mathrm{l}}_{\nu}<I^{\mathrm{b}}_{\nu}$. In our case,
$I^{\mathrm{b}}_{\nu}$ represents the height-dependent Ca ii H background wing
intensity.

###  3.5 Statistical equilibrium and approximation to HNEI

In a dynamical atmosphere, such as the solar chromosphere, the statistically
time-dependent state of the plasma is given by:

$\sum_{j\neq i}n_{j}\,P_{ji}-n_{i}\sum_{j\neq i}P_{ij}=\frac{\partial
n_{i}}{\partial t},$ (9)

where $n_{i}$, $n_{j}$ are the number density of particles in state $i$ or
$j$. The left-hand side gives the collision term of the kinetic equilibrium
equation and is known as the statistical equilibrium equation if time-
dependent changes are negligible. Then, $P_{ij}$ is the total rate of
transition (collisions and radiative) from state $i$ to $j$. Radiative
transfer codes such as RH solve only the statistical equilibrium equation,
neglecting time-dependent changes which are important for species with large
second-level excitation energy, such as hydrogen and helium (Carlsson & Stein,
2002; Golding et al., 2014). These time-dependent changes are available from
state-of-the-art radiative MHD simulations. Currently, it is not
computationally feasible to solve the exact radiative transfer problem at the
same time with the MHD equations. Therefore, we can use an approximation for
HNEI given by

$\sum_{j\neq i,c}n_{j}\,P_{ji}-n_{i}\sum_{j\neq i}P_{ij}=n_{c}\,P_{ci},$ (10)

solving only for the bound levels and fixing the ionized hydrogen number
density $n_{c}$ given by the Bifrost simulation which treats hydrogen in non-
equilibrium ionization. Here, $P_{ci}$ describes the transition rate out of
the continuum. We implemented this approximation to HNEI into the RH code when
solving the statistical equilibrium equation.

The inherent degeneracy of solutions to the statistical equilibrium equation
allows us to formulate solutions in a variety of ways (White, 1961). However,
not all of them are convenient for the interpretation of line source
functions. A particularly interesting solution is given by Jefferies (1960,
1968), written in terms of indirect transition probabilities between levels.
This formulation helps to find the strongest indirect transition contributing
to the interlocking term in the multi-level source function description (see
Sect. 3.6). A general formal solution in terms of level ratios was already
given by Rosseland (1926) and in terms of indirect transition probabilities
can be written as

$\frac{n_{u}}{n_{l}}=\frac{\sum_{k}P_{lk}\,q_{ku,l}}{\sum_{k}P_{uk}\,q_{kl,u}},$
(11)

with $P_{lk}$ the transition rate out of the lower level $l$ into individual
levels $k$ and $q_{ku,l}$ the indirect transition probability that a
transition starting from the level $k$ arrives at the upper level $u$ before
transitioning back to the lower level $l$. The indirect transition
probabilities can be calculated by solving a set of linear equations:

$q_{ij,k}=\sum_{l\neq i}p_{il}\,q_{lj,k},$ (12)

with the the transition probability $p_{il}$ given as:

$p_{il}=\frac{P_{il}}{\sum_{k}P_{ik}}.$ (13)

###  3.6 Multi-level source function

The expression for the solution of the statistical equilibrium equation given
by Eq. (11) can be used to express the multi-level line source,
$S^{\mathrm{l}}_{\nu}$, for the complete redistribution in a more intuitive
way:

$S^{\mathrm{l}}_{\nu_{0}}=\sigma\,\overline{J}_{\nu_{0}}+\epsilon\,B_{\nu_{0}}(T_{\mathrm{e}})+\eta\,B_{\nu_{0}}(T^{\star}),$
(14)

where the term $\sigma\overline{J}_{\nu_{0}}$ describes the contribution to
line photons by scattering from the profile-averaged mean radiation field
$\overline{J}_{\nu_{0}}$, while $\epsilon B_{\nu}(T_{\mathrm{e}})$ is the
fraction of thermally created line photon following the Planck function and
the last term, $\eta B_{\nu}(T^{\star}),$ is the contribution to line photons
from interlocking transitions. Interlocking transitions populate the upper or
lower level of a given transition indirectly via an intermediate level. Here,
we express the multi-level line source function $S^{\mathrm{l}}_{\nu_{0}}$ in
terms of total extinction, while in the literature (e.g., Jefferies, 1968) the
multi-level line source function is often expressed relative to scattering.
The coefficients $\sigma,\epsilon$, and $\eta$ relative to the total
extinction are given by:

$\displaystyle\kappa=A_{ul}+C_{ul}\left(1-\mathrm{exp}\left(-\frac{h\nu_{0}}{k_{\mathrm{B}}T_{\mathrm{e}}}\right)\right)+\left(g_{u}\sum_{u}-g_{l}\sum_{l}\right),$
(15) $\displaystyle\sigma=\frac{A_{ul}}{\kappa},$ (16)
$\displaystyle\epsilon=\frac{C_{ul}\left(1-\mathrm{exp}\left(-\frac{h\nu_{0}}{k_{\mathrm{B}}T_{\mathrm{e}}}\right)\right)}{\kappa},$
(17)
$\displaystyle\eta=\frac{\left(g_{u}\sum_{u}-g_{l}\sum_{l}\right)}{\kappa},$
(18)

with $C_{ul}$ and $A_{ul}$ having the usual meaning of the collisional and
radiative deexcitation rates. Also, $T_{\mathrm{e}}$ is the electron
temperature. The terms $\sum_{u},\sum_{l}$, and $B_{\nu}(T^{\star})$, which
determine interlocking, are given as:

$\displaystyle\sum_{u}=\sum_{j\neq l}P_{uj}\,q_{jl,u},$ (19)
$\displaystyle\sum_{l}=\sum_{j\neq u}P_{lj}\,q_{ju,l},$ (20) $\displaystyle
B_{\nu}(T^{\star})=\frac{2h\nu^{3}}{c^{2}}\left(\frac{g_{u}\sum_{u}}{g_{l}\sum_{l}}-1\right)^{-1}.$
(21)

The quantities $\sum_{u},\sum_{l}$ are the sum of all indirect transitions
that start from the upper level $u$ (lower level $l$) and populate indirectly
the lower level $l$ (upper level $u$) via some intermediate transitions. These
two quantities help determine which intermediate levels play a significant
role in the contribution to interlocking, whereas $B_{\nu}(T^{\star})$
represents a Planck function with a characteristic temperature $T^{\star}$
determined by the dominant intermediate interlocking processes.

To determine how strong each intermediate level contributes to interlocking
$B_{\nu}(T^{\star})$ in Eq. (21) we can further break down the expression into
a sum of Planck functions each with a characteristic temperature, $T^{\star}$,
as:

$B_{\nu}(T^{\star})=\sum_{i}B^{i}_{\nu}(T^{\star}_{i})\gamma_{i}=\sum_{i}B^{i}_{\nu}(T^{\star}_{i})\frac{\beta_{i}}{\beta_{\mathrm{tot}}},$
(22)

where $i$ is the sum over the intermediate levels and
$B^{i}_{\nu}(T^{\star}_{i})$, $\beta_{i}$, and $\beta_{\mathrm{tot}}$ are
given by:

$\displaystyle
B^{i}_{\nu}(T^{\star}_{i})=\frac{2h\nu^{3}}{c^{2}}\left(\frac{g_{u}\,P_{ui}\,q_{il,u}}{g_{l}\,P_{li}\,q_{iu,l}}-1\right)^{-1},$
(23) $\displaystyle\beta_{i}=g_{u}\,P_{ui}\,q_{il,u}-g_{l}\,P_{li}\,q_{iu,l},$
(24)
$\displaystyle\beta_{\mathrm{tot}}=\sum_{i}g_{u}\,P_{ui}\,q_{il,u}-g_{l}\,P_{li}\,q_{iu,l},$
(25)

where $B^{i}_{\nu}(T^{\star}_{i})$ describe Planck functions following
characteristic temperatures $T^{\star}$ determined by the radiative and
collisional processes coupling the upper and lower level through a loop with
the intermediate level $i$. Here, $\beta_{i}$ and $\beta_{\mathrm{tot}}$ are
necessary to split the ratio of $g_{u}\sum_{u}/g_{l}\sum_{l}$ into several
different Planck functions. Physically, $\beta_{i}/\beta_{\mathrm{tot}}$
describes the fractional contribution to interlocking extinction for each
intermediate loop via level $i$. $\beta_{i}$ the contribution to interlocking
extinction corrected for indirect transitions into the upper level for each
individual indirect transition loop.

Finding the physical process and transitions that set the line source function
in NLTE radiative transfer can be quite challenging if indirect transitions
play a significant role in populating atomic levels. It is difficult to
quantify how strong intermediate transitions affect the line source function,
especially for model atoms with a large number of transitions. The given
description of the multi-level source function, including indirect transition
probabilities for the interlocking term, gives a straightforward way to
quantify which transition(s) dominate the line source function.

Figure 1: H $\varepsilon$ line core observations of a quiet Sun and pore
region. The images are histogram equalized to the synthetic H $\varepsilon$
line core image in Fig. 6. _Left panel:_ Quiet Sun observation with a wide
band insert at the top right corner. The WB insert covers the area outlined by
the light blue square and includes a few magnetic bright points. _Right
panel:_ More active scene close to a pore region (taken at 08:49:45 UT) with a
zoomed-in insert at the top right corner (the insert is more saturated to
highlight dark fibrils). Crosses indicate positions of spectral profiles shown
in Fig. 2 with the same color coding.

##  4 Results

###  4.1 H $\varepsilon$ observations

In this section, we discuss how the morphology of H $\varepsilon$ line core
images can guide us to understand the formation of the spectral line. To
explain the apparent structures seen in H $\varepsilon$ line core images, we
need to determine where we are looking in the solar atmosphere, in terms of
extinction and what sets the total source function.

Figure 1 displays such H $\varepsilon$ line core images. We choose these two
particular regions for comparison with our synthetic spectra. Figure 2
compares synthetic degraded spectra with selected observed spectral profiles.

Figure 2: Ca ii H plus H $\varepsilon$ spectral profiles for different solar
structures indicated with crosses in Fig. 1. _Solid lines:_ Observations taken
with the CHROMIS instrument. _Dashed lines:_ Degraded synthetic spectra. The H
$\varepsilon$ emission profiles are shown in orange and light blue. Orange
synthetic profile (1.5D HNEI) is taken from Fig. 6 at $x\approx 13.6$ Mm and
$y\approx 17.7$ Mm (red contours). Light blue synthetic profile (1.5D HNEI)
shows the same as the one in the bottom right panel of Fig. 12. The red
spectral profile is taken from a magnetic bright point. The synthetic profile
(3D SE) is shown in the bottom left panel of Fig. 10. The turquoise spectral
profile is taken from a dark fibrilar structure. The synthetic profile (1.5D
HNEI) is shown in the bottom right panel of Fig. 10.

We start with a discussion of where we are looking with regard to extinction
in the quiet Sun case. The H $\varepsilon$ line core images reveal a reversed
granulation pattern. Reversed granulation originates from the mid-photosphere
due to a reversal in the temperature structure, which leads to a reversal in
the intensity images compared to the photospheric granulation pattern
(Nordlund, 1984; Cheung et al., 2007). The fact that the H $\varepsilon$ line
core image shows reversed granulation, as well as the Ca ii H blue wing (see
Fig. 3) suggests that the layers above reversed granulation are optically thin
to H $\varepsilon$ radiation. As we are looking at atmospheric layers close to
the temperature minimum, the H $\varepsilon$ extinction should be
significantly reduced due to the temperature sensitivity of the hydrogen $n=2$
population. A forward synthesis of H $\varepsilon$ from the Bifrost simulation
would help clarify which observed structures (if any) could be optically thick
to H $\varepsilon$ radiation. Furthermore, as no chromospheric canopy
structure (the fibrilar chromosphere that is characteristic for a
chromospheric line such as H $\alpha$) is apparent in H $\varepsilon$ quiet
Sun images, the chromosphere is mostly transparent, suggesting a mid-
photospheric origin of H $\varepsilon$.

Figure 3: H $\varepsilon$ and Ca ii H pore observation covering a region with
stronger magnetic fields (taken at 09:17:52 UT). The wavelength positions are
indicated in the top left corners on an average spectral profile with light
blue crosses. _Left panel:_ Ca ii H blue ($-0.15$ nm) wing image, showing no
chromospheric structures. _Middle panel:_ H $\varepsilon$ line core image,
showing both reversed granulation and chromospheric structures. _Right panel:_
Ca ii H line core image, dominated by chromospheric structures.

The pore region shows a slightly different picture. Here, we see some imprints
of dark fibrilar structures. This hints at a chromospheric origin for such
fibrilar structures observed in H $\varepsilon$ (see top-right insert in the
right panel of Fig. 1). To illustrate that H $\varepsilon$ can show some
chromospheric structures, we show the Ca ii H blue wing, H $\varepsilon$, and
Ca ii H line core images side-by-side in Fig. 3. The blue wing Ca ii H blue
image shows only reversed granulation, whereas the H $\varepsilon$ image shows
fine imprints of dark fibrils. These dark fibrils are connected to
chromospheric fibrilar structures observed in the Ca ii H line core image,
illustrating a chromospheric contribution in the H $\varepsilon$ images. Not
all dark fibrilar structures seem to be optically thick. Through some fibrilar
structures, we observe the background with reversed granulation. In more
active solar regions, the H $\varepsilon$ extinction appears to increase at
higher layers, making the chromosphere opaque. We note the images shown were
taken during excellent seeing conditions, making the thin dark fibrilar
structures clearly visible.

Next, we address the H $\varepsilon$ source function. Ayres & Linsky (1975)
suggested that the H $\varepsilon$ source function is dominated by the Balmer
continuum radiation field. This would result in apparent images resembling
granulation if H $\varepsilon$ is formed under optically-thick conditions, as
the Balmer continuum is highly scattering and the thermalization depth lies in
the photosphere. Nevertheless, we see reversed granulation in H $\varepsilon$
images, which suggests most radiation is coming from the Ca ii H background
intensity and not H $\varepsilon$ itself. Detailed radiative transfer
calculations can reveal which hydrogen atomic transition dominates the H
$\varepsilon$ source function.

The conclusion we can draw from the H $\varepsilon$ observations helps to
narrow down the question, how the H $\varepsilon$ line is formed in the solar
atmosphere and gives some restriction on the total source function (Eq. 1).
The total source function $S_{\nu}^{\mathrm{t}}$ has to be sensitive to
temperature but it is a linear combination of the line source function,
$S^{\mathrm{l}}_{\nu_{0}}$, and the background source function,
$S^{\mathrm{b}}_{\nu}$, weighted by their relative extinctions. This leaves us
with two options. Either the Ca ii H extinction dominates and
$S_{\nu}^{\mathrm{t}}$ is set by $S^{\mathrm{b}}_{\nu}$, which follows the
Planck function. Or, instead, the H $\varepsilon$ extinction dominates, and
$S_{\nu}^{\mathrm{t}}$ follows $S^{\mathrm{l}}_{\nu_{0}}$, which is coupled to
temperature through one of the terms in Eq. (14).

We address the formation of H $\varepsilon$ with the forward synthesis of H
$\varepsilon$ from state-of-the-art radiative MHD simulations. We specify
where we are looking in the solar atmosphere in terms of extinction and
determine the dominant processes in the total source function. We evaluate the
effects of Balmer line blanketing, HNEI, and 3D effects on the formation of H
$\varepsilon$. In addition, we discuss the appearance of bright points and
elongated dark fibrils apparent in the H $\varepsilon$ images with example
spectra shown in Fig. 2 and why and when we see H $\varepsilon$ in emission
against the Ca ii H background. We mainly focus on quiet Sun regions but will
try to give some idea about the formation of H $\varepsilon$ in more active
regions.

###  4.2 Balmer continuum

The Balmer continuum radiation is the main driver of the hydrogen rate system
and the dominant source of hydrogen ionization in the solar chromosphere
(Carlsson & Stein, 2002). Electrons recombine into high energy levels and
cascade downwards, strongly affecting the transition rates and level
populations of the hydrogen atom and, therefore, the H $\varepsilon$ line. It
is thus essential to properly model the Balmer continuum radiation field
including line blends. RH 1.5D has the option to add line blends in the form
of a Kurucz line list 111For details see:
http://kurucz.harvard.edu/linelists.html. The spectral lines from the Kurucz
line list are calculated assuming LTE, including scattering. To speed up the
computation we created a list of lines that affect the Balmer continuum
radiation significantly.

Figure 4: Synthetic Ca ii H and H $\varepsilon$ spectral profiles compared
against a FTS atlas observation. _Top panel:_ FTS atlas observations (orange)
compared to the synthesised Ca ii H plus H $\varepsilon$ profiles from the
FALC model. Blue profile for the case with no Balmer continuum blanketing,
light blue for Balmer blanketing with all lines included in the Kurucz line
list, and dark blue for the optimized Kurucz line list including only lines
affecting the Balmer continuum radiation significantly. _Bottom panel:_
Averaged Bifrost Ca ii H plus H $\varepsilon$ profiles for the SE (light blue)
and HNEI (dark blue) case against FTS atlas observations (orange).

Figure 4 compares the synthesized Ca ii H plus H $\varepsilon$ spectral region
from the FALC model atmosphere against the FTS atlas normalized to the
continuum. The inclusion of line blanketing in the Balmer continuum results in
stronger H $\varepsilon$ absorption closely matching the observations. The
bottom panel compares the FTS atlas against the average line profile summed
over the Bifrost simulation calculated with SE and HNEI. There is a strong
effect of the Balmer radiation field strength on the H $\varepsilon$ line
profile, but this does not necessarily mean that the H $\varepsilon$ source
function is set by the Balmer radiation field. It only highlights that the
hydrogen transitions (as H $\varepsilon$) are sensitive to changes in the
Balmer radiation field. To show that the Balmer radiation field sets the H
$\varepsilon$ source function, it is necessary to show that the divergence
between the upper and lower level of H $\varepsilon$ with height is set by the
Balmer radiation field. Further, we have to evaluate what sets the total
source function in Eq. (1) and which term dominates the multi-level source
function in Eq. (14) for H $\varepsilon$.

Figure 5: Total source function of H $\varepsilon$, effect of Balmer continuum
line blanketing, and H $\varepsilon$ multi-level source function for the FALC
model. _Upper panel:_ Line ($S_{H\epsilon}$), background ($S_{b}$), and total
source function ($S_{t}$), for the FALC model including Balmer continuum line
blanketing. The dotted grey line shows the ratio $\delta$ between line
extinction and total extinction (Eq. (2)). The three vertical ticks mark the
total optical depth heights $\tau_{\lambda}=3,1,0.3$. The $32$% label shows
that at the $\tau_{\lambda}=1$ height $32$% of the total source function
consists of the H $\varepsilon$ source function. _Middle panel:_ Fraction of
the total, line, and background source function including Balmer blanketing
and no Balmer blanketing. The change in the total source function with Balmer
blanketing results from a change in the H $\varepsilon$ source function, only
making up $32$% of the total source function. _Bottom panel:_ H $\varepsilon$
source function split up in the multi-level source function terms (Eq. (14);
scattering $\sigma\overline{J}$, thermal $\epsilon$B(Te), and interlocking
$\eta$B(Te)). The H $\varepsilon$ source function is dominated by interlocking
and not by the Balmer edge (dotted green), with the hydrogen ground level as
the dominant intermediate level ($\eta(\gamma$B1(T${}^{\star}_{1}$))).

Figure 5 presents what sets the total source function (top panel), the
influence of the Balmer line blanketing on the total source function (middle
panel), and what sets the H $\varepsilon$ source function for the FALC solar-
atmosphere model. The total source function is dominated by the Ca ii H
background source function throughout the lower solar atmosphere. At the
formation height of H $\varepsilon$, only around $32\%$ of the photons are
actually from the H $\varepsilon$ transition. This stems from the fact that Ca
ii H is a resonance line, whereas H $\varepsilon$ lower level is an excited
level with an energy jump of $10.2$ eV. The Ca ii H extinction will dominate
the total extinction in the cooler parts of the atmosphere, with the dip in
$\delta$ coinciding with the photospheric temperature minimum. This is an
indication that the solar layers above reversed granulation are optically thin
to H $\varepsilon$ radiation.

Figure 6: Synthetic H $\varepsilon$ line core images at $397.1202$ nm. _Left
panel:_ SE case. _Right panel:_ HNEI case. Contours mark regions of H
$\varepsilon$ emission at different column mass regions. _Yellow:_ column mass
$\geq$ $2\cdot 10^{-1}$ g cm-2. _Red:_ $2\cdot 10^{-1}$ g cm${}^{-1}>$ column
mass $>2\cdot 10^{-2}$ g cm-2. _Cyan:_ column mass $\leq$ $2\cdot 10^{-2}$ g
cm-2. _Dashed lines:_ position of the vertical slice shown in Fig. 7.

The strong effect on the emergent H $\varepsilon$ absorption stems from the
Balmer line blanketing changing the total source function and not the
extinction. The Balmer line blanketing has minor effects on the lower level
populations of the H $\varepsilon$ and Ca ii H transitions that make up the
extinction. As a result, the atmospheric layers we are looking at in the FALC
model do not change, but the number of photons available as described by the
source function does change. The middle panel illustrates that the change in
the total source function, resulting from the Balmer line blanketing is mainly
due to a change in the H $\varepsilon$ source function – and not the
background source function. Even when the total source function is dominated
by the Ca ii H source function. An increase in the strength of the Balmer
radiation field (no Balmer blanketing) will lead to an increase in the upper-
level population of the H $\varepsilon$ transition, therefore, increasing the
source function but keeping the H $\varepsilon$ extinction the same. However,
this does not show that the H $\varepsilon$ source function is set by the
Balmer radiation field, only that the H $\varepsilon$ source function is
somehow affected by the Balmer radiation field. To illustrate which process
dominates the H $\varepsilon$ source function we have to evaluate the multi-
level source function.

The bottom panel of Fig. 5 divides the H $\varepsilon$ source function into
the multi-level source function terms. The contributions from the scattering
$\sigma\overline{J}$ and thermal $\epsilon B(T_{\mathrm{e}}$) part are orders
of magnitude lower than the one from interlocking $\eta B(T_{\mathrm{e}}$).
Interlocking dominates the H $\varepsilon$ source function throughout the
whole FALC model atmosphere. To highlight that the H $\varepsilon$ source
function is actually not following the Balmer radiation field we plotted the
Balmer mean radiation field close to the Balmer edge in dotted green. The
dominant interlocking term in Eq. (22) is $\gamma_{1}B_{1}(T^{\star}_{1})$,
with the ground level as an intermediate level. The Planck function
$B_{1}(T^{\star}_{1})$ is set by the ratio $P_{71}q_{12,7}/P_{21}q_{17,2}$,
where the probability $q_{17,2}$ characterizes the H $\varepsilon$ source
function. The probability $q_{17,2}$ represents the sum of the transition
probabilities, $p_{1i}$, Eq. (13) from the ground level into other hydrogen
levels multiplied by the indirect transition probability to end up in the
upper level of H $\varepsilon$. The transition probabilities are dominated by
the radiative rates out of the ground level, not the collisional rates.

In summary, in the FALC model, the total source function is dominated by the
Ca ii H source function at the formation height of H $\varepsilon$. Our multi-
level source function description suggests that the H $\varepsilon$ source
function is dominated by interlocking, with the ground level being the
dominant intermediate level radiatively populating the upper level of H
$\varepsilon$ and not by the Balmer continuum radiation field as suggested by
Ayres & Linsky (1975).

###  4.3 Effect of HNEI on H $\varepsilon$

The solar chromosphere is pervaded by upwardly propagating magneto-acoustic
shocks, making the effect of time-dependent hydrogen ionization important
already above the solar photosphere. Non-equilibrium computations of hydrogen
ionization strongly affect the temperature structure and ionization of the
atmosphere due to the slow recombination rate at low-temperature intershock
phases (Leenaarts et al., 2007). The ion population in turn is strongly
coupled to the $n=2$ population of hydrogen, which sets the line extinction of
H $\varepsilon$ and therefore the formation height of Balmer series lines
(Carlsson & Stein, 2002). In some cases, HNEI can have strong effects
depending on where we are looking at in the atmosphere and what we are
observing, thereby directly influencing the H $\varepsilon$ source function.

Figure 6 shows synthesized H $\varepsilon$ line core (at $397.1202$ nm) images
from the Bifrost simulation (left panel for SE, and right panel for HNEI). The
SE case uses the standard version of RH to solve the statistical equilibrium
equation, which loses information about the time-dependent hydrogen
ionization. Our HNEI implementation in RH uses the proton density given by the
Bifrost simulation (see Sect. 3.5).

Figure 7: Hydrogen ionization (_top_) and $n=2$ populations (_bottom_) from a
vertical slice through the simulation (indicated in Fig. 6). _Left column_ :
Results for statistical equilibrium (SE). _Right column_ : Results for non-
equilibrium hydrogen ionization (HNEI). The solid orange and dashed grey lines
show the $\tau=1$ heights for both cases, and the orange crosses and light
blue crosses mark regions where H $\varepsilon$ is in emission. The dark blue
line shows the relative variation of the H $\varepsilon$ line core intensity.

The two images look almost identical. Both are mapping reversed granulation,
which originates in the mid-photosphere. The synthetic images match the
observed H $\varepsilon$ images from Fig. 1 quite well, which suggests a
similar formation of H $\varepsilon$ between simulation and observations.
However, the effect of HNEI seems negligible for the observed reversed
granulation pattern, which in the simulation is mainly made up of weak H
$\varepsilon$ absorption lines. To evaluate the effect of HNEI on the H
$\varepsilon$ source function and extinction, we look at a vertical cut
through the Bifrost simulation, denoted by two dashed lines in Fig. 6.

Figure 8: Departures from statistical equilibrium from a vertical slice
through the Bifrost simulation. _Top panel:_ Departure of the total source
function from SE. _Middle panel:_ Departure of the hydrogen $n=2$ populations
from SE. _Bottom panel:_ Temperature structure of the Bifrost simulation. The
vertical slice is taken at the position marked with white dashed lines in Fig.
6. _Red solid:_ H $\varepsilon$ $\tau=1$ height for SE. _White dashed:_ H
$\varepsilon$ $\tau=1$ height for HNEI.

Figure 7 compares the hydrogen ionization structure and $n=2$ level population
for the cases of SE and HNEI. The HNEI case shows a different and smoother
ionization structure through the atmosphere compared to the SE case. This is
also reflected in the $n=2$ populations – especially in the lower part of the
atmosphere (below $1$ Mm), where the dark voids vanish under HNEI. The reason
for the smoother $n=2$ populations is that under HNEI the $n=2$ populations
follow the hydrogen ionization (Leenaarts et al., 2007), whereas for SE, the
$n=2$ population follows the temperature (compare temperature panel in Fig.
8). This temperature sensitivity below $1$ Mm stems from the fact that
Lyman-$\alpha$ is in radiative balance $n_{1}/n_{2}=R_{21}/R_{12}$ with
$S\approx J\approx B$. Therefore, the $n=2$ hydrogen population follows LTE:
Saha-Boltzmann partitioning for the predominately neutral lower solar
atmosphere, with high Lyman-$\alpha$ extinction (Vernazza et al., 1981;
Rutten, 2016, 2017).

The $n=2$ level population increases strongly for hot pockets in the lower
solar atmosphere in SE. This raises the $\tau=1$ height, represented by the
orange line in Fig. 7. These elevated $\tau=1$ locations are connected to H
$\varepsilon$ emission indicated by crosses. Under HNEI the H $\varepsilon$
extinction does not “feel” this temperature increases as much as in SE. This
not only has a strong influence on the extinction but also on the source
function itself. The influence of HNEI on the $n=2$ populations and total
source function in terms of departure from SE can be seen in Fig. 8. We show
the departures from SE for the lower solar atmosphere in the formation region
of H $\varepsilon$. The $\tau=1$ heights for H $\varepsilon$, where most of
the radiation escapes are similar between SE and HNEI, which illustrates that
the increase in the $n=2$ populations under HNEI is not enough to change where
we are looking at in the simulation. Further, the total source functions at
the $\tau=1$ height do not change much between SE and HNEI. Therefore, time-
dependent hydrogen ionization effects are negligible for the H $\varepsilon$
extinction and the total source function where H $\varepsilon$ absorption is
formed. However, this is not the case for H $\varepsilon$ emission.

H $\varepsilon$ emission is co-located with enhanced temperatures in the lower
solar atmosphere, as illustrated in the lower panel of Fig. 8 (H $\varepsilon$
emission is marked with crosses in Fig. 7). This leads to two orders of
magnitude higher $n=2$ populations at these heating events, shifting up the
$\tau=1$ formation heights. At these particular atmospheric heights, the total
source function under SE is significantly higher (up to a factor of 2)
compared to HNEI, resulting in stronger H $\varepsilon$ emission. Only two out
of five emission locations (crosses in Fig. 7) are still in emission in HNEI.
Even though the $\tau=1$ height is lower in HNEI, for these two locations in
emission, under SE, the H $\varepsilon$ emission profile has a significant
contribution from these higher-lying layers, where the total source function
is increased.

HNEI has a sizeable effect on the $n=2$ populations throughout the whole
atmosphere, from the lower atmosphere up to the transition region. It changes
not only the extinction but also the H $\varepsilon$ source function. This
effect can be seen in Fig. 8 (top panel). Similarly, we find that HNEI also
affects the line source functions of other lines in the Balmer series (not
shown in the Figures). Depending on where the lines are formed, the difference
between the emergent Balmer series line profiles in SE or HNEI can be
significant. Another important side effect of HNEI is time-dependent
ionization, where the hydrogen populations can have a memory of a previous
atmospheric state, decoupling the observed line profiles from the
instantaneous conditions in the atmosphere. This has implications for the
interpretation of H $\alpha$ spectroheliograms, as they contain information
about previous events, especially heating events that increase the H $\alpha$
extinction (Rutten, 2016).

The strong effect of HNEI on H $\varepsilon$ emission can be seen in Fig. 6,
where the over-plotted contours mark regions where H $\varepsilon$ is in
emission. The different colored contours group the column masses at $\tau=1$
height into three different groups related to the column mass panel of Fig. 9.
The yellow contours highlight regions where the column mass at $\tau=1$ is
greater than or equal to $2\cdot 10^{-2}$ g cm-1, red regions greater than
$2\cdot 10^{-1}$ g cm-2 and less than $2\cdot 10^{-2}$ g cm-2, and the cyan
color less than or equal to $2\cdot 10^{-2}$ g cm-2. Figure 9 shows a kernel
density estimation (Rosenblatt, 1956) for the distributions of formation
heights, column masses, temperatures, and electron densities at the $\tau=1$
heights for H $\varepsilon$ emission and absorption lines for the Bifrost
simulation. H $\varepsilon$ in absorption is formed just above the Ca ii H
background, whereas H $\varepsilon$ emission is formed across a large range of
column masses. This range is connected to enhanced temperatures in the lower
solar atmosphere, which shifts the formation height upwards under SE and is
connected to increased ionization in HNEI. Generally, HNEI reduces the amount
of H $\varepsilon$ emission formed higher up in the atmosphere.

In summary, H $\varepsilon$ synthetic images are similar to solar
observations, showing mainly reversed granulation. Why we observe reversed
granulation in the H $\varepsilon$ line core, we will address in the Sect.
4.4. Furthermore, we have to explain the small bright points seen in the H
$\varepsilon$ line core images; HNEI has little effect on the reversed
granulation pattern but influences the H $\varepsilon$ emission. As H
$\varepsilon$ emission can form relatively high up in the lower solar
atmosphere, we have to check for probable 3D radiative transfer effects (Sect.
4.5).

Figure 9: Distributions of atmospheric quantities at the H $\varepsilon$
$\tau=1$ height. Panels show the kernel density estimation as a function of
height, column mass, temperature, and electron density at the $\tau=1$ height
for the H $\varepsilon$ rest wavelength. The distributions are separated by
the type of H $\varepsilon$ profile, as indicated in the legend: absorption in
HNEI, emission in HNEI, emission in SE, and the Ca ii H background continuum
at H $\varepsilon$ rest wavelength in HNEI.

###  4.4 H $\varepsilon$ absorption

In this section, we discuss the formation of H $\varepsilon$ absorption
profiles. We analyze why we observe reversed granulation in weak H
$\varepsilon$ absorption profiles and why strong H $\varepsilon$ absorption
lines are linked to magnetic elements.

Figure 10 shows a synthetic H $\varepsilon$ line core image from a small
region from the Bifrost simulation. The region covers a magnetic element and
parts of the reversed granulation pattern. To check for 3D effects we
synthesized the same region in 3D, but assuming SE and treating Ca ii H in
CRD. We adopted these approximations because it was difficult to get RH to
converge under HNEI and PRD in regions with strong horizontal inhomogeneities,
such as magnetic elements. Representative profiles are shown in Fig. 10
including one particularly interesting structure, which goes from H
$\varepsilon$ emission in 1.5D to absorption in 3D (this structure crosses the
red cross at an approximately 45∘ angle). We address the effects of 3D
radiative transfer on H $\varepsilon$ emission in Sect. 4.5.

Figure 10: Synthetic Ca ii H and H $\varepsilon$ images and spectral profiles
for reversed granulation, magnetic elements, and a dark fibrilar structure.
_Topmost row:_ Synthesized H $\varepsilon$ line core images from a cutout of
the Bifrost simulation calculated in 1D (_left_) and 3D (_right_). _Bottom two
rows:_ Ca ii H plus H $\varepsilon$ spectral profiles from locations marked
with coloured crosses in the H $\varepsilon$ line core images. _Middle-left:_
Granular spectral profile. _Middle-right:_ Intergranular spectral profile.
_Bottom-left:_ Spectral profile located at a magnetic field concentration.
_Bottom-right:_ Location where H $\varepsilon$ emission in 1D turns to
absorption in 3D. Figure 11: H $\varepsilon$ four-panel diagrams for an
intergranular lane (_left_) and a magnetic bright point (_right_). The
locations of each point in the simulation are marked with orange and blue
crosses in Fig. 10. Each diagram is organized as follows. _Top left panels:_
Relative $\chi_{\nu}/\tau_{\nu}$; sensitive to velocity gradients. _Top right
panels:_ $S_{\nu}$: Relative source function. _Solid white:_ Total source
function at rest wavelength. _Dashed yellow:_ H $\varepsilon$ source function.
_Dashed blue:_ Background source function. Source functions are displayed in
brightness temperature units with the scale at the top. _Solid grey:_ Ratio of
H $\varepsilon$ to total extinction (if H $\varepsilon$ extinction dominates
the grey curve is on the left side of the plot; if Ca ii H extinction
dominates on the right). _Dotted red:_ Atmospheric temperature profile.
_Bottom-left panels:_ Relative $\tau_{\nu}e^{(-\tau_{\nu})}$; indicates where
most of the contribution to relative line intensity comes from. _Solid green:_
$\epsilon B_{\nu}(T_{e})$. _Solid violet:_ $\sigma J_{\nu}$. _Solid gold:_
$\eta B_{\nu}(T^{\star})$. The three terms are displayed in brightness
temperature units with the scale at the top. _Bottom-right panels:_ $C_{I}$;
contribution to relative line depression. _Solid white:_ Ca ii H plus H
$\varepsilon$ line profile. _Dashed white:_ Ca ii H background line profile.
_Solid red:_ Atmospheric vertical velocity profile. _Solid light blue:_
$\tau=1$ height.

In general, we observe very weak H $\varepsilon$ absorption lines inside
granules that get stronger the closer we get to the granular edge. 3D effects
are negligible on the H $\varepsilon$ line profiles (see middle panels of Fig.
10). To discuss the formation of weak H $\varepsilon$ absorption lines we
created relative four-panel formation diagrams inspired by Carlsson & Stein
(1997). Figure 11 displays such diagrams, which we henceforth refer to as
“four-panel diagrams.” The four-panel diagrams decompose the contribution to
relative line depression or emission into three components:

$C_{\mathrm{I}}=\frac{\mathrm{d}I_{\mathrm{R}}(\nu,z)}{\mathrm{d}z}=\frac{\chi^{\mathrm{R}}_{\nu}(z)}{\tau^{\mathrm{R}}_{\nu}(z)}\cdot\tau^{\mathrm{R}}_{\nu}(z)\,\mathrm{exp}(-\tau^{\mathrm{R}}_{\nu}(z))\cdot
S^{\mathrm{R}}_{\nu}(z),$ (26)

based on Eq. (5). The first component
$\chi^{\mathrm{R}}_{\nu}(z)/\tau^{\mathrm{R}}_{\nu}(z)$ is the ratio between
relative extinction and relative optical depth. This term has large values
where there are many absorbing or emitting particles (line and background;
$\chi^{\mathrm{l}}_{\nu}+\chi^{\mathrm{b}}_{\nu}$) at low relative optical
depth. This component is sensitive to line-of-sight velocity gradients that
move plasma, so that they effectively increase extinction at varying Doppler
offsets at wavelength regions with low optical depths.

The second component
$\tau^{\mathrm{R}}_{\nu}(z)\,\mathrm{exp}(-\tau^{\mathrm{R}}_{\nu}(z))$ has
the the largest values around the $\tau=1$ height. In this panel, we over-
plotted the three components of the multi-level source function (Eq. 14) for H
$\varepsilon$. This helps to determine the dominant line-formation process.

The third component $S^{\mathrm{R}}_{\nu}(z)$ describes the relative source
function. The relative source function can have positive or negative values.
This is a major difference from standard four-panel formation diagrams, where
the source function panel only has positive values. Here, positive values
signify that the line source function, $S^{\mathrm{l}}_{\nu}$, is lower than
the background intensity, $I^{\mathrm{b}}_{\nu}$, resulting in a contribution
to relative line depression at a particular height in the atmosphere. Negative
values signify a relative contribution to line emission. However, the relative
source function does not specify what sets the total source function at the
formation height of H $\varepsilon$. Based on the upper panel of Fig 5, we
included the different components expressing the total source function at the
H $\varepsilon$ line core wavelength: the total source function, the H
$\varepsilon$ source function, the background source function, and
additionally included the ratio of H $\varepsilon$ extinction to total
extinction. The last panel shows $C_{I}$, which is the contribution to the
relative line depression or emission. We note that $C_{I}$ is positive in
regions that contribute to line absorption and negative in regions that
contribute to line emission.

First, we focus on the left four-panel diagram in Fig. 11, which illustrates
the formation of a relatively weak H $\varepsilon$ absorption line in the
intergranular lane (orange cross in Fig. 10). Most contribution to relative
line depression is formed just above the Ca ii H background intensity. This
highlights that the added H $\varepsilon$ extinction on top of the Ca ii H
extinction is very small, making the atmosphere nearly transparent to H
$\varepsilon$ radiation. Further, we have to evaluate what sets the total
source function at the line formation region specified by the relative
extinction over-plotted in the relative source function panel. The grey curve
points out that more than $50\%$ of line photons are actually coming from the
Ca ii H transition. The total source function is more tightly coupled to the
Ca ii H source function than to the H $\varepsilon$ source function. The Ca ii
H source function is still strongly coupled to LTE, which means that in this
case more than $50\%$ of photons at the H $\varepsilon$ wavelength are formed
under LTE conditions. The optical depth $\tau_{\lambda}$ of the atmosphere
above the height where the Ca ii H background is formed is just below unity
($\tau\approx 0.78$ for the intergranular location), thus giving the
appearance of reversed granulation in the H $\varepsilon$ line core, which
happens for most for quiet Sun locations.

The closer we get to the granular center, the weaker the H $\varepsilon$
absorption gets (see the middle row in Fig. 10). This stems from the fact that
the mid-photosphere is cooler at the center of granules where reversed
granulation is formed. As the H $\varepsilon$ extinction is strongly coupled
to Saha-Boltzmann due to Lyman-$\alpha$ thermalization, the H $\varepsilon$
contribution to the total extinction becomes smaller, essentially removing the
H $\varepsilon$ line depression. The decrease in H $\varepsilon$ extinction
will couple the total source function even closer to the Ca ii H source
function ($\approx$80%).

The change in H $\varepsilon$ absorption strength from the granular center to
the granular edge is therefore dependent on the temperature structure and not
so much on the strength of the Balmer radiation. The Balmer radiation will
indeed increase the H $\varepsilon$ source function in both SE and HNEI. In
SE, a stronger Balmer radiation field will not change the H $\varepsilon$
extinction ($n=2$ population), but increase the $n=7$ populations, thereby
increasing the H $\varepsilon$ source function. In HNEI, the behavior is
different: the H $\varepsilon$ extinction and $n=7$ populations decrease for a
stronger Balmer radiation field but still result in an increased H
$\varepsilon$ source function as the H $\varepsilon$ extinction decreases more
rapidly with height than the $n=7$ populations. However, the H $\varepsilon$
source function for weak H $\varepsilon$ absorption is dominated by
interlocking illustrated by the gold solid line in the
$\tau^{\mathrm{R}}_{\nu}(z)\,\mathrm{exp}(-\tau^{\mathrm{R}}_{\nu}(z))$ panel.
Similarly, as in Sect. 4.2, our multi-level source function description
suggests the interlocking term is dominated by the Lyman series with indirect
transitions into the upper level H $\varepsilon$ through the ground level –
and not by the Balmer radiation field. The Balmer radiation field will affect
the mean radiation fields in the Lyman series, but not the H $\varepsilon$
source function directly.

We now move to the formation of strong H $\varepsilon$ absorption lines. As a
representative example, we consider the bottom left line profile in Fig. 10
(blue cross), and show its corresponding four-panel diagram in Fig. 11 (right
panels). This strong H $\varepsilon$ absorption line is connected to a
magnetic field concentration. Due to the low gas density inside magnetic
elements, H $\varepsilon$ forms at lower heights, therefore, it is associated
with higher temperatures. As the Ca ii H wing intensity is well approximated
by LTE and follows the temperature, the intensity above the magnetic elements
will be significantly increased compared to that in reversed granulation. The
reason is similar to why magnetic elements appear bright in the wings of H
$\alpha$ (Leenaarts et al., 2006). However, H $\varepsilon$ is formed on top
of the Ca ii H wing. The added H $\varepsilon$ extinction shifts the formation
height upwards making magnetic elements optically thick to H $\varepsilon$
radiation. The total source function strongly decreases with height, being
mainly set by the H $\varepsilon$ source function due to the increased H
$\varepsilon$ extinction above magnetic elements. This mapping of higher
layers with a lower source function compared to the Ca ii H background leads
to strong H $\varepsilon$ absorption lines. Therefore, H $\varepsilon$ can be
used to identify and track weak magnetic elements.

Our multi-level source function description also suggests that the H
$\varepsilon$ source function is dominated by interlocking. The dominant
interlocking process is an indirect transition loop through the ground level
for strong H $\varepsilon$ absorption profiles connected to magnetic field
concentrations, the same process as found for weak H $\varepsilon$ absorption
profiles.

The fourth profile, marked with a red cross in Fig. 10 is a good example that
some H $\varepsilon$ profiles can suffer significant 3D effects. In this
particular case, the bright elongated structure seen in 1D will turn into a
faint fibrilar structure if 3D radiative transfer effects are taken into
account. In the next section, we address 3D effects on H $\varepsilon$
emission regions, where the formation height is increased compared to the
usual reversed granulation.

###  4.5 3D effects on H $\varepsilon$ emission

The Balmer series lines, for example H $\alpha$, are strongly affected by
scattering. The characteristic fibrilar chromosphere seen in H $\alpha$ can
only be reproduced if the radiation field is treated in 3D geometry (Leenaarts
et al., 2012a). However, as shown in Sect. 4.2, the H $\varepsilon$ source
function is not dominated by two-level scattering but instead by interlocking
through the Lyman series. The Lyman series lines are the strongest scattering
lines in the solar spectrum, having the lowest collisional destruction
probability per line photon extinction. Hence, we evaluated 3D Lyman
scattering effects on the emergent H $\varepsilon$ line profiles formed higher
up in the atmosphere (at lower column masses than $2\cdot 10^{-2}$ g cm-2).

Figure 12: Synthetic Ca ii H and H $\varepsilon$ images and spectral profiles
for H $\varepsilon$ emission formed at low mass densities. _Top two rows:_
Synthesized H $\varepsilon$ line core images from a cutout of the Bifrost
simulation calculated in 1D (_top row_) and 3D (_middle row_). Contours mark
column mass regions, same as in Fig. 6. _Bottom row:_ Ca ii H plus H
$\varepsilon$ spectral profiles from locations marked with red crosses in the
H $\varepsilon$ line core images. Figure 13: H $\varepsilon$ four-panel
diagrams for the locations marked with crosses in Fig. 12. _Top row:_ H
$\varepsilon$ synthesized in 1.5D. _Bottom row:_ Same location synthesized in
3D. _Top left panels:_ Relative $\chi_{\nu}/\tau_{\nu}$; sensitive to velocity
gradients. _Top right panels:_ $S_{\nu}$; relative source function. Yellow
colormap for negative source function values. Light blue colormap for positive
source function values. _Solid white:_ Total source function at rest
wavelength. _Dashed yellow:_ H $\varepsilon$ source function. _Dashed blue:_
Background source function. Source functions are displayed in brightness
temperature units with the scale at the top. _Solid grey:_ Ratio of H
$\varepsilon$ to total extinction (if H $\varepsilon$ extinction dominates the
grey curve is on the left side of the plot; if Ca ii H extinction dominates on
the right). _Dotted red:_ Atmospheric temperature profile. _Bottom left
panels:_ Relative $\tau_{\nu}e^{(-\tau_{\nu})}$; indicates where most of the
contribution to relative line intensity comes from. _Solid green:_ $\epsilon
B_{\nu}(T_{e})$. _Solid violet:_ $\sigma J_{\nu}$. _Solid gold:_ $\eta
B_{\nu}(T^{\star})$. The three terms are displayed in brightness temperature
units with the scale at the top. _Bottom right panels:_ $C_{I}$; contribution
to relative line depression/emission. Yellow colormap contribution to line
depression. Light blue colormap contribution to line emission. _Solid white:_
Ca ii H plus H $\varepsilon$ line profile. _Dashed white:_ Ca ii H background
line profile. _Solid red:_ Atmospheric upward velocity profile. _Solid light
blue:_ $\tau=1$ height.

To investigate 3D effects on the emergent H $\varepsilon$ profiles, we chose
small regions of interest from the Bifrost simulation. From these regions, we
synthesized H $\varepsilon$ with the assumption of HNEI in 1D and 3D geometry.
We carried out spectral synthesis in 3D only in these selected regions to
avoid the large computational costs of a full 3D synthesis, and because there
are only a few regions where H $\varepsilon$ is formed higher up in the
atmosphere. Furthermore, we synthesized Ca ii H in CRD instead of PRD in 3D
geometry as the running with PRD led to numerical instabilities and the
problem often failed to converge. Because we assume CRD in 3D, the intensity
at the Ca ii H inner wings will deviate between 3D and 1.5D.

First, we discuss the case of low column mass, where 3D effects become
important at certain locations. Figure 12 shows two regions of interest from
the Bifrost simulation, where H $\varepsilon$ is in emission at low column
mass. The left column shows a region with strong 3D effects, while the right
column shows a region with negligible 3D effects. We show the corresponding
four-panel diagrams in Fig. 13.

The left column of Fig. 13 illustrates the formation of H $\varepsilon$ with
strong 3D effects. In the 3D case, $C_{\mathrm{I}}$ outlines significantly
more contribution to line depression than for the 1.5D case on top of the
relative line emission region. The answer to the question of why we have more
contribution to line depression in 3D is twofold. First the H $\varepsilon$
extinction $\chi^{\mathrm{l}}_{\nu}$ plays an important role (see. Eq. (8)).
The $n=2$ populations are in radiative equilibrium with Lyman-$\alpha$ as
$n_{2}=n_{1}R_{12}/R_{21}$ is set by the Lyman-$\alpha$ mean radiation field
and the exponential density decrease of the atmosphere. The Lyman-$\alpha$
mean radiation field in 3D is slightly increased compared to 1.5D. Therefore,
the H $\varepsilon$ extinction is mapping more of the positive relative source
function giving rise to relative line depression. However, the increase in H
$\varepsilon$ extinction is not the main reason for the enhanced contribution
to the relative line depression in 3D highlighted by similar
$\tau_{\nu}\mathrm{exp}(-\tau_{\nu})$ panels.

The main reason why we see enhanced relative line depression in 3D is that the
H $\varepsilon$ line source function is significantly lower than the 1.5D case
above a certain atmospheric height. This decrease in $S^{\mathrm{l}}_{\nu}$ in
3D causes a negative relative line source function at lower heights compared
to 1.5D (illustrated in the top right $S_{\nu}$ panels in Fig. 13). The
relative source function becomes negative when the H $\varepsilon$ source
function $S^{\mathrm{l}}_{\nu}$ is smaller then the background Ca ii H wing
intensity, specified by $1-S^{\mathrm{l}}_{\nu}/I^{\mathrm{b}}_{\nu}$. The
background intensity $I^{\mathrm{b}}_{\nu}$ flattens out with height (around
$\approx 0.5$ Mm), having lower values in 3D. Our multi-level source function
description suggests that the H $\varepsilon$ source function is dominated by
interlocking, set by the Lyman series (same process as for H $\varepsilon$
absorption), or more precisely, a linear combination of the transition
probability, $p_{1i}$, multiplied by an indirect transition probability,
$q_{i7,2}$, due to a Lyman series transition. The transition probabilities are
dominated by the radiative rates. The net rate out of the ground level is set
by the Lyman-$\alpha$ radiative rate which is orders of magnitude higher than
other Lyman series radiative rates. Therefore, the transition probabilities
are given by $p_{1i}=R_{1i}/R_{12}$ and reflect the ratio between the Lyman
series and the Lyman-$\alpha$ mean radiation field. The main difference
between 1.5D and 3D is that this ratio decreases more strongly with height in
3D than in 1.5D. This puts the H $\varepsilon$ source function,
$S^{\mathrm{l}}_{\nu}$, below the Ca ii H background intensity,
$I^{\mathrm{b}}_{\nu}$, at lower heights giving rise to a negative source
function creating the strong relative absorption contribution in 3D.

The negligible 3D effects on H $\varepsilon$ for the right column of Fig. 13
stem from the fact that the temperature increase is situated lower in the
atmosphere than for the left column. This temperature rise at lower heights
leads to a coupling of the Lyman-$\alpha$ mean radiation field with the Planck
function to heights where the temperature is decreasing again, above the
temperature peak. Whereas for the case of strong 3D effects, the
Lyman-$\alpha$ mean radiation decouples from the Planck function at the
temperature peak in the atmosphere. In both cases, the Lyman-$\alpha$ mean
radiation decouples approximately at the same height from the Planck function.
As the $n=2$ level population follows the Lyman-$\alpha$ mean radiation via
$n_{2}=n_{1}R_{12}/R_{21}$, the H $\varepsilon$ line extinction
$\chi^{\mathrm{l}}_{\nu}$ decreases after the temperature rise. This strong
decrease in line extinction, as illustrated in the
$\tau_{\nu}\mathrm{exp}(-\tau_{\nu})$ panel, is the main reason why we do not
see a strong contribution to relative line depression above the emission
region. The H $\varepsilon$ emission formed at greater column masses than
$2\cdot 10^{-2}$ g cm-2 (yellow and red contours in Fig. 6) show negligible 3D
effects.

In summary: if a temperature increase in the lower solar atmosphere extends
too high up in the atmosphere and is not confined in the lower solar
atmosphere (for the Bifrost simulation this would be approximately a height of
$0.75$ Mm), the H $\varepsilon$ line will suffer from 3D effects, showing a H
$\varepsilon$ absorption line instead of a H $\varepsilon$ emission line in
1.5D, given the density stratification of the Bifrost atmosphere.

This explains the strong change in the appearance of the bright fibrilar
structure in 1.5D in Fig. 10, compared with the dark fibrilar structure seen
in 3D. Figure 8 show a temperature cut through this fibrilar structure (white
dashed lines in Fig. 6). The strong temperature enhancement throughout the
lower atmosphere connected to this structure leads to an increased H
$\varepsilon$ line extinction via $n_{2}=n_{1}R_{12}/R_{21}$. The 3D effects
lead to a lower H $\varepsilon$ line source function than the Ca ii H
background intensity $I^{\mathrm{b}}_{\nu}$ and give a strong contribution to
relative line depression above the emission region compared to the 1.5D case.
As a result, the 1.5D H $\varepsilon$ emission lines at the fibrilar structure
become absorption lines in 3D.

##  5 Discussion

###  5.1 The H $\varepsilon$ source function

Numerical experiments such as the Bifrost simulation represent a solar-like
model atmosphere with properties close to the real solar atmosphere. We focus
our discussion on the similarities between synthetic and observed H
$\varepsilon$ profiles. The synthetic H $\varepsilon$ profiles combined with
the Bifrost simulation can help us understand the basic formation processes in
the solar atmosphere that lead to features observed in the quiet Sun, such as
reversed granulation, bright points, H $\varepsilon$ emission, and dark
fibrilar structures.

Several decades ago, Thomas (1957) suggested dividing the resonance lines and
strong lines such as H $\alpha$ into collisionally or photoelectrically
dominated. For H $\alpha$, this division meant that the source function is set
by collisions, the second term in Eq. (14), or by the third term
(photoelectric) in Eq. (14). The third term describes the interlocking
contribution to the line source function and, following Thomas (1957), it
states that the indirect path that dominates the H $\alpha$ source function is
photoionization from $n=2$ and a recombination cascade into $n=3$ (thus, the
description of it as photoelectric). However, Rutten & Uitenbroek (2012)
demonstrated that the H $\alpha$ source function is dominated by scattering,
described by the first term in Eq. (14), with strong contributions from
chromospheric backscattering. However, our main question concerns which
process dominates the source function of higher order Balmer series lines,
such as H $\varepsilon$.

Ayres & Linsky (1975) first modeled the formation of H $\varepsilon$ based on
1D numerical radiative transfer calculations for two model atmospheres, one
representing the Sun (H $\varepsilon$ in absorption) and one Arcturus (H
$\varepsilon$ in emission). They concluded that the Balmer continuum radiation
dominates the H $\varepsilon$ source function via photoionization, the same
mechanism proposed by Thomas (1957) for the H $\alpha$ source function.
However, our H $\varepsilon$ forward modeling from the FALC and Bifrost model
atmosphere suggests otherwise. Our multi-level source function description
suggests that the H $\varepsilon$ source function is dominated by interlocking
via the ground level of hydrogen due to the Lyman series. Not the Balmer
continuum radiation field.

Our synthetic H $\varepsilon$ line core images (Fig. 6) reproduce the observed
H $\varepsilon$ reversed granulation pattern (Fig. 1) quite well. Both show
the characteristic dark granules and bright intergranular lanes suggesting
that the Bifrost simulation has the basic physical properties representing a
quiet Sun lower atmosphere. The reversed granulation pattern consists of weak
H $\varepsilon$ absorption lines. The weakest H $\varepsilon$ absorption
profiles can be found inside granules and the closer we get to the
intergranular lanes, the stronger the H $\varepsilon$ absorption becomes. We
illustrated that the total source function for weak H $\varepsilon$ absorption
lines is dominated by the background Ca ii H source function and not the H
$\varepsilon$ source function. The temperature structure of the reversed
granulation sets the extinction and total source function and therefore the
behavior of weak H $\varepsilon$ absorption lines.

We tested the influence of HNEI on the emergent H $\varepsilon$ line profiles
against SE. The influence is negligible at the reversed granulation pattern,
which is represented by weak H $\varepsilon$ absorption lines, formed deep in
the solar atmosphere where time-dependent ionization is not important.
However, for hydrogen lines formed higher up in the solar atmosphere, HNEI
will have severe effects on the line extinctions and source functions. The
main effect of HNEI is that the proton density sets the level populations and
not the temperature structure. Therefore, structures not reflecting reversed
granulation should indicate high ionization regions (leading to increased
extinction), showing higher atmospheric layers than reversed granulation. The
most prominent features seen sticking out of the reversed granulation
background in Fig. 1 and Fig. 3 are bright points and dark fibrilar
structures.

###  5.2 Bright points

The bright points in our synthetic and observational H $\varepsilon$ images
are connected to small magnetic elements. The wide band insert in Fig. 1
illustrates that the magnetic structures are connected to H $\varepsilon$
bright points. We explained the formation of H $\varepsilon$ above magnetic
elements in Sect. 4.4 as the result of density depletion and a different
optical depth scale compared to their surroundings. Therefore, we get large Ca
ii H wing intensities (formed under LTE) with strong H $\varepsilon$
absorption lines on top. The strong H $\varepsilon$ absorption forms higher in
the atmosphere, with a strong decrease of the total source function versus
height that is dominated by the H $\varepsilon$ source function. If we compare
our synthetic H $\varepsilon$ profile in Fig. 2 with observations of a
magnetic element (red profiles), they look remarkably similar. Therefore, we
conclude that bright points in H $\varepsilon$ observations represent magnetic
elements where the Ca ii H wing intensities reflect the temperature at the Ca
ii H wing formation height, and the strong H $\varepsilon$ absorption lines
reflect an enhanced H $\varepsilon$ extinction above magnetic elements where
the H $\varepsilon$ source function drops rapidly with increasing height.

###  5.3 Dark fibrilar structures

We could only reproduce a dark fibrilar structure with the Bifrost simulation
if we include 3D radiative transfer effects; still, these structures are rare
in the synthetic spectra. One example is shown in Fig. 10, where the bright
fibrilar structure in 1D turns into a (very faint) dark fibrilar structure if
3D effects are taken into account. This feature goes diagonally through the
red cross in Fig. 10, extending nearly to the top and right borders of the
image. The 3D radiative transfer effects lead to stronger relative absorption
contribution above emission (see Fig. 13) compared to 1D. The reason is the
temperature structure of the atmosphere. The temperature structure (or proton
densities under HNEI) sets the $n=2$ populations due to the Lyman-$\alpha$
mean radiation field. If we have high-enough temperatures in the lower
atmosphere, up to heights where the Lyman-$\alpha$ decouples from the Planck
function, it will increase the H $\varepsilon$ line extinction. This increases
the contribution from higher layers, where the H $\varepsilon$ source function
is lower than the Ca ii H background intensity, thereby giving rise to strong
contributions to absorption. Such dark fibrilar structures could in principle
be visible in solar observations and not only in synthetic H $\varepsilon$
images from simulation. They would mark regions with strong heating: high
temperatures that extend from the lower atmosphere up to the lower
chromosphere and are not confined to the lower atmosphere.

It is important to avoid confusing the Bifrost dark fibrilar structure shown
in Fig. 10 with the dark fibrilar structures observed in more active regions
(Fig. 3). These fibrilar structures become partly opaque due to increased mass
density similar to the H $\alpha$ fibrilar chromosphere (Leenaarts et al.,
2012a). These fibrilar structures outline the chromospheric canopy. There is a
clear correlation between the dark fibrilar structures seen in the Ca ii H and
H $\varepsilon$ line core images. On the other hand, the dark fibrilar
structure from Bifrost is a result of heating. Dark fibrilar structures
similar to those seen in the simulation should be observable in the quiet Sun,
but they may be easily confused with other dark structures connected to the
reversed granulation background. Dark fibrilar structures should stand out
clearly from the Ca ii H reversed granulation and not be aligned with the
reversed granulation pattern. One example of this is the dark fibrilar
structure seen in Fig. 10 connecting two different granules. It is conceivable
that a strong episode of flux emergence may be an event that shows clear
signatures of dark fibrilar structures.

###  5.4 H $\varepsilon$ in emission

Lastly, we address locations of H $\varepsilon$ emission that mark regions of
confined heating in the lower atmosphere (otherwise 3D effects would become
important and we would observe a H $\varepsilon$ absorption line). Locations
with H $\varepsilon$ in emission are difficult to identify only by looking at
H $\varepsilon$ line core images. We have to look at the individual line
profiles to identify locations with emission. In Fig. 6, we outlined H
$\varepsilon$ emission locations with contours. All locations covered by
contours in the simulation are related to temperature enhancements in the
lower atmosphere, at heights between $0.3$ Mm and $0.75$ Mm in the Bifrost
simulation (see an example in Fig. 13). We found a significant amount of H
$\varepsilon$ emission in the Bifrost simulation, but a question arises about
whether such emission can also be found in observations covering quiet Sun
regions, as H $\varepsilon$ is generally thought to be a weak absorption
feature in the quiet Sun. To illustrate that the quiet Sun does show regions
with emission, we show two example profiles in Fig. 2 and compare them to some
of our synthetic profiles.

Observing H $\varepsilon$ in emission requires high resolution. We interpret
the regions with H $\varepsilon$ emission as locations marking confined
heating events: weak temperature enhancements in the lower solar atmosphere,
as illustrated by the formation of H $\varepsilon$ emission from the Bifrost
simulation. We are not aware of any spectral line that shows such a clear
signature of weak small-scale heating ($\approx 6000$ K) in the lower solar
atmosphere. The Mg ii triplet lines are sensitive to such weak heating events
formed at even lower column masses (Pereira et al., 2015) and can serve as a
fitting complement to H $\varepsilon$.

##  6 Conclusions

We revisited the formation of H $\varepsilon$, inspired by the work of Ayres &
Linsky (1975), and we propose that the H $\varepsilon$ source function is not
set by the Balmer radiation field but instead by interlocking via the ground
level. We present high-resolution H $\varepsilon$ line core images that appear
to be dominated by reversed granulation, mostly coming from the Ca ii H
background radiation as the contribution of H $\varepsilon$ is very small
under typical quiet Sun conditions. In several cases, we find locations with a
different appearance: bright points or dark fibrilar structures that are
related to enhanced temperatures and reflect locations where a majority of
photons come from the H $\varepsilon$ transition. These structures are
optically thick to H $\varepsilon$ radiation. Therefore, H $\varepsilon$ line
core images represent a stacked view of different atmospheric layers where the
dominant source of photons varies with hydrogen ionization.

Our main findings are as follows:

  * •

In most quiet Sun locations, the H $\varepsilon$ extinction is minor and the
atmospheric layers above the formation height of Ca ii H are transparent to H
$\varepsilon$ radiation.

  * •

In magnetic elements and small-scale heating events, the H $\varepsilon$
extinction dominates and the atmospheric layers are optically thick to H
$\varepsilon$ radiation.

  * •

Our multi-level source function description suggests that the H $\varepsilon$
source function is dominated by interlocking via the term
$P_{71}q_{12,7}/P_{21}q_{17,2}$, with the hydrogen ground level as
intermediate interlocking level – and not by the Balmer continuum radiation.

  * •

The H $\varepsilon$ source function is dominated by the Lyman series via the
term $q_{17,2}$, which represents a linear combination of the transition
probabilities out of the ground level.

  * •

The H $\varepsilon$ source function and extinction are both coupled to
temperature in the lower solar atmosphere. The extinction via Lyman-$\alpha$,
and the source function via the whole Lyman series.

  * •

HNEI couples the hydrogen populations to the proton densities, whereas SE
couples the population to temperature in the lower solar atmosphere. HNEI has
a strong effect on all hydrogen lines, affecting both source functions and
extinctions.

  * •

We can synthesize H $\varepsilon$ in 1.5D with the assumption of HNEI as long
as H $\varepsilon$ is not formed too high up in the atmosphere, where the 3D
effect can become important.

  * •

HNEI has the strongest effect at locations where H $\varepsilon$ is formed
higher up in the solar atmosphere, such as magnetic elements, dark fibrilar
structures, and H $\varepsilon$ emission locations.

  * •

H $\varepsilon$ emission marks locations of confined heating events in the
lower solar atmosphere.

The H $\varepsilon$ transition is special compared to H $\alpha$ and H $\beta$
because the extinction and source function are coupled to temperature. For H
$\alpha$ and H $\beta$, only the extinctions are coupled to temperature, while
their source functions are dominated by scattering, except by strong heating
events in the lower atmosphere such as Ellerman bombs (Ellerman, 1917; Rutten
et al., 2013). In H $\alpha$ and H $\beta$, the line core extinction is high
enough that we observe the chromospheric magnetic canopy. For H $\varepsilon$,
this is not the case. This gives the unique opportunity for a Balmer series
line to look through the magnetic canopy and observe heating events in the
lower solar atmosphere. Another interesting fact is that the source function
of higher-order Balmer series lines (higher than H $\alpha$ and H $\beta$)
seem to become dominated by interlocking via the Lyman series, and therefore
temperature-sensitive. This extinction “gap” and temperature sensitivity of
the H $\varepsilon$ source function opens a new window onto the possibility of
observing small-scale weak heating events in the lower atmosphere marked by H
$\varepsilon$ emission regions, for instance, QSEB reconnection events (Rouppe
van der Voort et al., 2016; Joshi et al., 2020), shocks, and others in high
spatial and temporal cadence. Another positive aspect is that it is relatively
easy to observe H $\varepsilon$ together with Ca ii H. For example, the
CHROMIS instrument at the SST can observe both lines without sacrificing too
much temporal cadence making the combination of spectral lines a strong
diagnostic tool. The H $\varepsilon$ and Ca ii H combination is also well
covered by the ViSP instrument (de Wijn et al., 2022) at DKIST (Rimmele et
al., 2020).

The H $\varepsilon$ line is not only of solar interest but also for the wider
astrophysics community. Generally, H $\varepsilon$ is used to study stellar
flare activity (Pavlenko et al., 2019), chromospheric activity in M and K
dwarfs (Houdebine et al., 2009; Maldonado et al., 2017), active galactic
nuclei (AGNs) (Ilić et al., 2012), young active stars (Biazzo et al., 2009),
Cepheids (Kovtyukh et al., 2015), and chromospheric active binaries (Montes et
al., 1995, 1996). Montes et al. (1996) showed a clear surface flux correlation
between Ca ii H and H $\varepsilon$ excess emission for chromospherically
active binaries. We can explain this correlation qualitatively. The Ca ii H
source function for these types of stars is collisionally dominated, formed
under LTE conditions, and therefore the intensity is proportional to
temperature. Meanwhile, the H $\varepsilon$ source function is set by
interlocking with a partial sensitivity to temperature and therefore the
intensity is coupled to the temperature structure of the atmosphere. The
sensitivity of both lines to temperature results in a strong correlation
between the two line excess fluxes and explains the steep temperature increase
needed to produce H $\varepsilon$ emission against the Ca ii H background
intensity.

We demonstrate that the H $\varepsilon$ can be used to detect weak small-scale
heating in the lower solar atmosphere at an unprecedented resolution. This
could be of particular interest in the studies of QSEBs and it is worthwhile
exploring how H $\varepsilon$ could help to shed light on these new phenomena
in greater detail. Further questions that ought to be considered in future
studies include: what higher order Balmer series lines can tell us about QSEBs
and whether stellar H $\varepsilon$ emission forms the same way as in the
solar atmosphere as we have predicted, as well as when (and why) the Balmer
series lines actually become dominated by interlocking.

###### Acknowledgements.

This work has been supported by the Research Council of Norway through its
Centers of Excellence scheme, project number 262622, and through project
number 325491. Computational resources have been provided by Sigma2 – the
National Infrastructure for High-Performance Computing and Data Storage in
Norway. We acknowledge funding support by the European Research Council under
ERC Synergy grant agreement No. 810218 (Whole Sun). K. K. acknowledges the
insight Mats Carlsson gave during his Stellar Atmospheres II course. The
Swedish 1-m Solar Telescope is operated on the island of La Palma by the
Institute for Solar Physics of Stockholm University in the Spanish
Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de
Canarias. The Institute for Solar Physics is supported by a grant for research
infrastructures of national importance from the Swedish Research Council
(registration number 2021-00169). The Scientific color maps (Crameri, 2021)
are used in this study to prevent visual distortion of the data and exclusion
of readers with colourvision deficiencies (Crameri et al., 2020).

## References

  * Anan et al. (2017) Anan, T., Ichimoto, K., & Hillier, A. 2017, A&A, 601, A103
  * Ayres & Linsky (1975) Ayres, T. R. & Linsky, J. L. 1975, ApJ, 201, 212
  * Bayazitov & Sakhibullin (1991) Bayazitov, U. S. & Sakhibullin, N. A. 1991, Soviet Astronomy Letters, 17, 160
  * Biazzo et al. (2009) Biazzo, K., Frasca, A., Marilli, E., et al. 2009, A&A, 499, 579
  * Bjørgen et al. (2018) Bjørgen, J. P., Sukhorukov, A. V., Leenaarts, J., et al. 2018, A&A, 611, A62
  * Carlsson et al. (2016) Carlsson, M., Hansteen, V. H., Gudiksen, B. V., Leenaarts, J., & De Pontieu, B. 2016, A&A, 585, A4
  * Carlsson & Leenaarts (2012) Carlsson, M. & Leenaarts, J. 2012, A&A, 539, A39
  * Carlsson & Stein (1997) Carlsson, M. & Stein, R. F. 1997, ApJ, 481, 500
  * Carlsson & Stein (2002) Carlsson, M. & Stein, R. F. 2002, ApJ, 572, 626
  * Cheung et al. (2007) Cheung, M. C. M., Schüssler, M., & Moreno-Insertis, F. 2007, A&A, 461, 1163
  * Crameri (2021) Crameri, F. 2021, Scientific colour maps, Zenodo
  * Crameri et al. (2020) Crameri, F., Shephard, G. E., & Heron, P. J. 2020, Nature Communications, 11, 5444
  * de Jager (1952) de Jager, C. 1952, Recherches Astronomiques de l’Observatoire d’Utrecht, 1
  * de Wijn et al. (2022) de Wijn, A. G., Casini, R., Carlile, A., et al. 2022, Sol. Phys., 297, 22
  * Ellerman (1917) Ellerman, F. 1917, ApJ, 46, 298
  * Evershed (1930) Evershed, J. 1930, MNRAS, 90, 762
  * Fontenla et al. (1993) Fontenla, J. M., Avrett, E. H., & Loeser, R. 1993, ApJ, 406, 319
  * Golding et al. (2014) Golding, T. P., Carlsson, M., & Leenaarts, J. 2014, ApJ, 784, 30
  * Gudiksen et al. (2011) Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154
  * Gurtovenko et al. (1974) Gurtovenko, E., Ratnikova, V., & de Jager, C. 1974, Sol. Phys., 37, 43
  * Houdebine et al. (2009) Houdebine, E. R., Junghans, K., Heanue, M. C., & Andrews, A. D. 2009, A&A, 503, 929
  * Ilić et al. (2012) Ilić, D., Popović, L. Č., La Mura, G., Ciroi, S., & Rafanelli, P. 2012, A&A, 543, A142
  * Jefferies (1960) Jefferies, J. T. 1960, ApJ, 132, 775
  * Jefferies (1968) Jefferies, J. T. 1968, Spectral line formation
  * Joshi et al. (2020) Joshi, J., Rouppe van der Voort, L. H. M., & de la Cruz Rodríguez, J. 2020, A&A, 641, L5
  * Kovtyukh et al. (2015) Kovtyukh, V., Szabados, L., Chekhonadskikh, F., Lemasle, B., & Belik, S. 2015, MNRAS, 448, 3567
  * Kurucz et al. (1984) Kurucz, R. L., Furenlid, I., Brault, J., & Testerman, L. 1984, Solar flux atlas from 296 to 1300 nm
  * Leenaarts et al. (2007) Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625
  * Leenaarts et al. (2012a) Leenaarts, J., Carlsson, M., & Rouppe van der Voort, L. 2012a, ApJ, 749, 136
  * Leenaarts et al. (2012b) Leenaarts, J., Pereira, T. M. D., & Uitenbroek, H. 2012b, A&A, 543, A109
  * Leenaarts et al. (2006) Leenaarts, J., Rutten, R. J., Sütterlin, P., Carlsson, M., & Uitenbroek, H. 2006, A&A, 449, 1209
  * Löfdahl et al. (2021) Löfdahl, M. G., Hillberg, T., de la Cruz Rodríguez, J., et al. 2021, A&A, 653, A68
  * Magain (1986) Magain, P. 1986, A&A, 163, 135
  * Maldonado et al. (2017) Maldonado, J., Scandariato, G., Stelzer, B., et al. 2017, A&A, 598, A27
  * Montes et al. (1995) Montes, D., de Castro, E., Fernandez-Figueroa, M. J., & Cornide, M. 1995, A&AS, 114, 287
  * Montes et al. (1996) Montes, D., Fernandez-Figueroa, M. J., Cornide, M., & de Castro, E. 1996, A&A, 312, 221
  * Nikaidou (1982) Nikaidou, Y. 1982, Sol. Phys., 80, 259
  * Nordlund (1984) Nordlund, A. 1984, in Small-Scale Dynamical Processes in Quiet Stellar Atmospheres, ed. S. L. Keil, 181
  * Pavlenko et al. (2019) Pavlenko, Y. V., Suárez Mascareño, A., Zapatero Osorio, M. R., et al. 2019, A&A, 626, A111
  * Pereira et al. (2015) Pereira, T. M. D., Carlsson, M., De Pontieu, B., & Hansteen, V. 2015, ApJ, 806, 14
  * Pereira & Uitenbroek (2015) Pereira, T. M. D. & Uitenbroek, H. 2015, A&A, 574, A3
  * Popper (1956) Popper, D. M. 1956, ApJ, 123, 377
  * Rimmele et al. (2020) Rimmele, T. R., Warner, M., Keil, S. L., et al. 2020, Sol. Phys., 295, 172
  * Rolli & Magun (1995) Rolli, E. & Magun, A. 1995, Sol. Phys., 160, 29
  * Rolli et al. (1998a) Rolli, E., Wülser, J. P., & Magun, A. 1998a, Sol. Phys., 180, 361
  * Rolli et al. (1998b) Rolli, E., Wülser, J. P., & Magun, A. 1998b, Sol. Phys., 180, 343
  * Rosenblatt (1956) Rosenblatt, M. 1956, The Annals of Mathematical Statistics, 27, 832
  * Rosseland (1926) Rosseland, S. 1926, ApJ, 63, 218
  * Rouppe van der Voort et al. (2016) Rouppe van der Voort, L. H. M., Rutten, R. J., & Vissers, G. J. M. 2016, A&A, 592, A100
  * Rutten (2016) Rutten, R. J. 2016, A&A, 590, A124
  * Rutten (2017) Rutten, R. J. 2017, A&A, 598, A89
  * Rutten & Uitenbroek (2012) Rutten, R. J. & Uitenbroek, H. 2012, A&A, 540, A86
  * Rutten et al. (2013) Rutten, R. J., Vissers, G. J. M., Rouppe van der Voort, L. H. M., Sütterlin, P., & Vitas, N. 2013, in Journal of Physics Conference Series, Vol. 440, Journal of Physics Conference Series, 012007
  * Scharmer (2006) Scharmer, G. B. 2006, A&A, 447, 1111
  * Scharmer et al. (2003) Scharmer, G. B., Bjelksjo, K., Korhonen, T. K., Lindberg, B., & Petterson, B. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4853, Innovative Telescopes and Instrumentation for Solar Astrophysics, ed. S. L. Keil & S. V. Avakyan, 341–350
  * Scharmer et al. (2019) Scharmer, G. B., Löfdahl, M. G., Sliepen, G., & de la Cruz Rodríguez, J. 2019, A&A, 626, A55
  * Thomas (1957) Thomas, R. N. 1957, ApJ, 125, 260
  * Turova (1994) Turova, I. P. 1994, Sol. Phys., 150, 71
  * Uitenbroek (2001) Uitenbroek, H. 2001, ApJ, 557, 389
  * van Noort et al. (2005) van Noort, M., Rouppe van der Voort, L., & Löfdahl, M. G. 2005, Sol. Phys., 228, 191
  * Vernazza et al. (1981) Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635
  * Wellmann (1940) Wellmann, P. 1940, ZAp, 19, 236
  * White (1961) White, O. R. 1961, ApJ, 134, 85
  * Wilson (1938) Wilson, O. C. 1938, PASP, 50, 245
  * Wilson (1957) Wilson, O. C. 1957, ApJ, 126, 46
  * Zapiór et al. (2022) Zapiór, M., Heinzel, P., & Khomenko, E. 2022, ApJ, 934, 16