Bandit Pareto Set Identification: the Fixed Budget Setting

Cyrille Kone Emilie Kaufmann Laura Richert

Institution 1 Institution 2 Institution 3  ]

We study a multi-objective pure exploration problem in a multi-armed bandit model. Each arm is associated to an unknown multi-variate distribution and the goal is to identify the distributions whose mean is not uniformly worse than that of another distribution: the Pareto optimal set. We propose and analyze the first algorithms for the fixed budget Pareto Set Identification task. We propose Empirical Gap Elimination, a family of algorithms combining a careful estimation of the “hardness to classify” each arm in or out of the Pareto set with a generic elimination scheme. We prove that two particular instances, EGE-SR and EGE-SH, have a probability of error that decays exponentially fast with the budget, with an exponent supported by an information theoretic lower-bound. We complement these findings with an empirical study using real-world and synthetic datasets, which showcase the good performance of our algorithms.

§ INTRODUCTION

The multi-armed bandit problem has been extensively studied in the literature, predominantly as a single-objective stochastic optimization problem. In this framework, an agent sequentially collects samples from a set of $K$ (unknown) probability distributions, called arms. Two common goals are either to maximize the sum of observations (regret minimization, <cit.>) or to identify the arm with the largest expected value (best arm identification, <cit.>). In many real-world problems, however, there are multiple (possibly conflicting) objectives to optimize simultaneously and there might not exist a unique arm that maximizes all these objectives. This leads to a scenario in which the observations are multi-dimensional, and there may be several Pareto optimal arms. An arm belongs to the Pareto optimal set $\cS^\star$ if it is not dominated by any other arm (in the sense that the mean values of all the objectives would be larger for the latter). While some line of work aims at minimizing the Pareto regret by strategically selecting arms from the Pareto set along the agent's trajectory <cit.>, we are interested in the pure exploration counterpart of this problem, first studied by <cit.>, in which the goal is to identify (relaxations of) the set $\cS^\star$, as quickly and accurately as possible. This multi-objective exploration problem has a wide range of applications including adaptive clinical trials (with multiple indicators of efficacy), software and hardware design (with multiple criteria such as speed and energy consumption), or A/B/n testing for recommender systems (with different metrics of users engagement).

Notably, prior work <cit.> has studied the Pareto Set Identification (PSI) problem under the fixed-confidence setting: given a risk $\delta \in (0, 1)$, the agent stops at a random time $\tau$ and recommends a set $\widehat S_\tau$ such that $\bP(\widehat S_\tau \neq \cS^\star) \leq \delta$, while minimizing its expected stopping time $\bE[\tau]$. However, to the best of our knowledge, no algorithms have been proposed for the dual fixed-budget setting in which the agent should output a guess for the Pareto set after a given budget $T$ while minimizing its probability of error $\bP(\widehat S_T \neq \cS^\star)$. For some applications, the fixed-budget setting can be viewed as more practical: e.g., in clinical trials the maximal number of patients is typically defined in advance (due to some financial constraints), and the stopping time used in the fixed-confidence formulation could be prohibitively large.

We fill this gap by proposing the first algorithm(s) for Pareto Set Identification in the fixed-budget setting. Our approach may be viewed as an extension of Successive Rejects (SR) <cit.> and Sequential Halving (SH) <cit.>, two state-of-the-art algorithms for fixed-budget best arm identification (i.e. PSI in the uni-dimensional setting), that are based on arm eliminations. We propose a novel arm elimination criterion based on the computation of an Empirical Gap. This quantity can be viewed as a measure of the hardness to classify an arm as Pareto optimal or sub-optimal and is inspired by the gaps featured in the analysis of fixed-confidence algorithms <cit.>. While SR and SH successively eliminate arms with smallest empirical means, in Empirical Gap Elimination (EGE), the arms with the largest empirical gaps are eliminated and appropriately classified as optimal or sub-optimal.

We prove an upper bound on the error probability of a generic version of EGE using any arm allocations satisfying some budget constraints (<ref>). We explicit our bound for two instances using the same arm allocations as SR and SH respectively, called EGE-SR and EGE-SH, showing that their error probability decrease exponentially fast with the budget (<ref>). We derive an information-theoretic lower bound showing that these algorithms are essentially optimal in the worse case (<ref>). On the practical side, we report in Section <ref> some experimental results showing that EGE-SR/SH outperform uniform sampling and are much more robust than , an adaptation of an existing fixed-confidence algorithm that requires some prior knowledge of the problem complexity as input (discussed in Appendix <ref>). Finally, we propose in Section <ref>, , a variant of EGE-SR that tackles a relaxation of PSI first considered by <cit.> in which one should identify at most $k$ Pareto optimal arms, with a provably reduced probability of error.

Related work
The Best Arm Identification problem (BAI) is one of the most studied pure exploration problems in the bandit literature. The fixed-confidence formulation of BAI is by now well understood as there are efficient algorithms with guaranteed bounds on their sample complexity, $\bE[\tau]$ (e.g., <cit.>) and some algorithms that attain the minimal sample complexity on any (parametric) problem instance in the asymptotic regime $\delta \rightarrow 0$ (e.g., <cit.>). In contrast, the theoretical understanding of fixed-budget BAI is more elusive. While the error probability of any reasonable strategy (e.g. uniform sampling) decays exponentially fast with the budget $T$, <cit.> shows that there is no optimal exponent that can be attained by a single algorithm on all bandit instances.
Still, the Successive Reject and Sequential Halving algorithms mentioned above have their error probability scaling in $\exp(-T/\log(K)H(\nu))$ for some complexity quantity $H(\nu)$ depending on the bandit instance $\nu$ and <cit.> proved a (worse-case) lower bound showing that there exists instances $\nu$ in which this rate is un-improvable (up to constant factors in the exponent) when the budget $T$ is large.

In recent years, there has been a growing interest in pure exploration problems in a multi-dimensional setting. Adaptive algorithms for Pareto Set Identification were first proposed in the literature on black-box multi-objective optimization. Besides several evolutionary heuristics (e.g., <cit.>), the works of <cit.> have studied PSI under a Gaussian Process modeling assumption, obtaining some error bounds that decay polynomially in $T$ and depend on some notion of information gain that is specific to the GP setting. In the bandit literature, the first fixed-confidence PSI algorithm proposed by <cit.> is based on uniform sampling and an accept/reject mechanism. <cit.> then proposed a more adaptive LUCB-like algorithm for PSI and some relaxations introduced therein. <cit.> studied an extension of PSI to any $\bR^D$ partial order defined by an ordering cone, still in the fixed confidence setting.

Alongside PSI, other multi-dimensional pure exploration problems have been investigated such as feasible arm identification <cit.> in which the goal is to identify the arms which belong to a given polyhedron $\cP\subset \bR^D$ or top feasible arm identification <cit.> in which one seeks a feasible arm that further maximizes some linear combinations of the different objectives. While the former work considers the fixed-budget setting and also proposes some counterpart of Successive Rejects, the latter considers the fixed-confidence setting. In the fixed-budget setting, a special case was recently considered by <cit.> in a  2-dimensional setting in which the goal is to find the arm maximizing one attribute under the constraint that the second attribute is larger than a given threshold. They propose a SR like algorithm that is similar in spirit to our Empirical Gap Elimination, but with a very different notion of empirical gap.

§ SETTING

We formalize the fixed budget Pareto Set Identification problem and introduce some notation.

We are given $K$ distributions $\nu_1, \dots, \nu_K$ over $\bR^D$ with means (resp.) $\vmu_1, \dots, \vmu_K$, called arms. Let $\nu:=(\nu_1, \dots, \nu_K)$ be the bandit instance and $\bm\Theta := (\vmu_1 \dots \vmu_K)^\T$. We use Boldfaced symbols for $\bR^D$ elements. As in prior work prior on multi-objective pure exploration, we assume that the marginal distributions of each arm are all $\sigma$-subgaussian.
We recall that a random variable $X$ is $\sigma$-subgaussian if $\forall\; u\in \bR$, $\bE[\exp\left( u (X -\bE[X]) \right)] \leq \exp\left( {u^2\sigma^2}/2\right)$.

Given two arms $i,j \in [K]$, $i$ is weakly (Pareto) dominated by $j$ (denoted by $\vmu_i \leq  \vmu_j$) if for any $d\in \{1, \dots, D\}$, $\mu_{i}^d \leq \mu_j^d$. The arm $i$ is (Pareto) dominated by $j$ ($\vmu_i \preceq \vmu_j$ or $i\prec j$) if $i$ is weakly dominated by $j$ and there exists $d \in \{1, \dots, D\}$ such that $\mu_i^d < \mu_j^d$. The arm $i$ is strictly (Pareto) dominated by $j$ ($\vmu_i\prec \vmu_j$ or $i \prec j$) if for any $d \in \{1, \dots, D\}$, $\mu_i^d < \mu_j^d$.

The Pareto set $\cS^\star(\nu)$ is
\begin{equation*}
\cS^\star(\nu) := \left\{ i \in [K] \tm{ s.t } \nexists j \in [K]: \vmu_i  \prec \vmu_j \right\},
\end{equation*}

and will be denoted by $\cS^\star$ when $\nu$ is clear from the context. Any arm $a\in \cS^\star$ will be called (Pareto) optimal and an arm $a\notin \cS^\star$ is called sub-optimal.

At each time $t=1, 2, \dots,T$, the agent chooses an arm $a_t$ and observes an independent draw $\bX_t \sim \nu_{a_t}$ with $\bE(\bX_{a_t}) = \vmu_{a_t}$. We denote by $\bP_{\nu}$ the law of the stochastic process $\{\bX_t\}_{t\geq 1}$ and by $\bE_{\nu}$, the expectation under $\bP_{\nu}$. Let $\cF_t:= \sigma(a_1, \bX_1, \dots, a_t, \bX_t)$ be the $\sigma$-algebra representing the history of the process up to time $t$ and $\cF:= \{\cF_t\}_{t\geq 1}$, the filtration of the process. The agent's strategy is adaptive in the sense that $a_t$ is $\cF_{t-1}$ measurable. At time $T$, the agent outputs a guess $\widehat{S}_T$ (which is $\cF_T$ measurable) for $\cS^\star$ which should minimize the error probability $e_T(\nu):= \bP_{\nu}(\widehat S_T\neq \cS^\star(\nu))$.

To characterize the difficulty of the problem, we introduce
below some quantities that allow to measure how much an arm is dominated. For any arms $i,j$, we let

\begin{eqnarray*}
\m(i,j) &:=& \min_d \left[\mu_j^d - \mu_i^d\right], \\
\M(i,j) &:=& \max_d \left[ \mu_i^d - \mu_j^d\right].
\end{eqnarray*}

If $\vmu_i\nprec \vmu_j$, $\M(i,j)$ is the smallest uniform increase of $j$ that makes it dominate $i$. If $\vmu_i\prec \vmu_j$ then $\m(i,j)$ is the smallest increase of any component of $i$ which makes it non-dominated by $j$.
<cit.> proposed an algorithm in the fixed-confidence setting whose sample complexity is characterized by some sub-optimality gaps $\Delta_i$'s defined as follows. For a sub-optimal arm $i\notin \cS^\star$,

\begin{equation}
\label{eq:def-gap-sub}
\Delta_i := \Delta_i^\star:= \max_{j \in \cS^\star} \m(i,j),
\end{equation}

which is the smallest quantity that should be added component-wise to $\vmu_i$ to make $i$ appear Pareto optimal w.r.t $\{\vmu_k: k\in [K]\setminus \{i\} \}$.
For an optimal arm $i \in \cS^\star$,

\begin{equation}
\label{eq:def-gap-opt}
\Delta_i :=
\min (\delta_i^+, \delta_i^-)
\end{equation}

\begin{eqnarray*}
\delta_i^+ &:=& \min_{j\in \cS^\star \setminus \{i\}} \min(\M(i,j), \M(j,i))\;, \\
\delta_i^- &:=& \min_{j\in [K]\setminus \cS^\star}\{(\M(j,i))^+ + \Delta_{j}\}\;,
\end{eqnarray*}
with the convention $\min_\emptyset = +\infty$.
$\delta_i^+$ accounts for how much $i$ is close to dominate (or to be dominated by) another optimal arm while $\delta_i^-$ translates in a sense the smallest “margin" from an optimal arm $i$ to the sub-optimal arms.  <cit.> show that the difficulty (i.e. near-optimal sample complexity) of fixed-confidence PSI is characterized by the complexity term

$$ H(\nu) :=  \sum_{a=1}^K  \frac{1}{\Delta_{a}^{2}}.$$

In this work, we show that it is also a relevant complexity measure for the fixed-budget setting.

§ ALGORITHMS

In this section, we introduce the family of Empirical Gap Elimination (EGE) algorithms.

§.§ Empirical Gap Elimination

An Empirical Gap Elimination algorithm uses a round-based structure. The algorithm is parameterized by a number of rounds $R>0$, an arm schedule vector $\bm\lambda = (\lambda_1,\dots,\lambda_R,\lambda_{R+1})\in [K]^{R+1}$ satisfying $\lambda_{r} > \lambda_{r+1}$ where $\lambda_r$ indicates how many arms are active in round $r$, and an allocation vector $\bm t= (t_1,\dots,t_R) \in [T]^{R}$, where $t_r$ is the number of samples gathered from each active arm in round $r$. Given the maximal budget $T$, these vectors should further satisfy the following:

\begin{align}
&\lambda_1 = K \quad \text{and}\quad \lambda_{R+1} \in \{0, 1\}
\label{eq:cond1} \\
&\text{and} \quad \sum_{r=1}^R \lambda_r t_r \leq T,	\label{eq:cond2}
\end{align}

ideally with an equality.

In each round $r$, the algorithm maintains a set of active arm, denoted by $A_r$, that is of size $\lambda_r$, and collects $t_r$ new samples from each arm in $A_r$. The total number of samples from each arm in $A_r$ gathered in the first $r$ rounds[We emphasize that EGE algorithms do not discard samples between rounds] is therefore

$$n_r = \sum_{s= 1}^r t_s.$$

For each $a\in A_r$, the empirical estimate of $\bm\theta_a$ at the end of round $r$ is denoted by
$\vmuh_a(r) := \frac{1}{n_r}\sum_{s=1}^{n_r} \bX_{a, s}$
where $\bX_{a,s}$ denote the $s$-th observation drawn from distribution $\nu_a$. These estimates are used to carefully decide which arm to explore in the next round, based on an appropriate notion of empirical gap.

These gaps are empirical variants of the (fixed-confidence) gaps introduced in the previous section. We first propose a rewriting of these quantities to remove the explicit dependency on $\cS^\star(\nu)$.

For any arm $i \in [K]$,

\[\Delta_i = \left\{\begin{array}{ll}
\Delta_i^\star = \max_{j\in [K]}\m(i, j) & \text{ if } i\notin \cS^\star \\
\delta_i^\star & \text{ if } i\in \cS^\star    \\
\end{array}\right.\;,
\]

where $\delta_i^\star := \min_{j\neq i} [\M(i, j)\land (\M(j, i)^+ +(\Delta_j^\star)^+)].$

We introduce the empirical quantities

\begin{eqnarray*}
\mh(i,j; r) &:=& \min_{d}\left[\muh_j^d(r) -  \muh_i^d(r)\right], \\
\Mh(i,j; r) &:=& \max_d \left[(\muh_i^d(r) - \muh_j^d(r))\right].
\end{eqnarray*}

and the empirical Pareto set

\begin{eqnarray*}S_r &:=& \{i \in A_r : \nexists j \in A_r : \vmuh_i(r) \prec \vmuh_j(r) \},\\
& = &\{i \in A_r: \forall j  \in A_r \backslash \{i\}, \Mh(i,j;r) > 0 \}\;.
\end{eqnarray*}

Finally, we define for any arm $i \in A_r$

\begin{eqnarray*}
\widehat \Delta_{i, r}^\star &:=& \max_{j\in A_r\backslash\{ i\}}\mh(i,j;r), \\
\widehat \delta_{i,r}^\star &:=& \!\!\min_{j \in A_r \backslash\{i\}} [\Mh(i,j;r) \;\land\; (\Mh(j,i;r)^+ + (\widehat \Delta_{i, r}^\star)^+)]
\end{eqnarray*}

and for any $i\in A_r$

\begin{eqnarray}
\label{eq:def-emp-gap}
\widehat \Delta_{i,r} := \begin{cases}
\widehat \Delta_{i, r}^\star & \text{ if } i \in A_r \backslash{S_r},\\
\widehat\delta_{i,r}^\star & \text{ else. }
\end{cases}
\end{eqnarray}

Empirical Gap Elimination (EGE)

Pareto set
parameters $R, \bm t, \bm \lambda$
Set $A_1 = \{1, \dots, K\}, B_1 := \emptyset$, $D_1 = \emptyset$
$r=1, 2, \dots, R$

Collect $t_r$ samples from each arm $a \in A_r$

Compute $S_r$ the empirical Pareto set

Let $A_{r+1}$ be the set of $\lambda_{r+1}$ arms in $A_r$ with the smallest empirical gaps $\widehat\Delta_{a,r}$

ties broken in favor of arms in $S_r$

$ B_{r+1} \gets B_r \cup \{ S_r \cap (A_r \backslash A_{r+1})\}$

$ D_{r+1} \gets D_r \cup \{  (A_r \backslash A_{r+1}) \backslash S_r \}$

$\widehat S_T = B_{R+1}  \cup A_{R+1}$

At the end of round $r$, the algorithm sorts the arms by increasing order of their empirical gaps: $\widehat{\Delta}_{(1),r} \leq \dots \leq  \widehat{\Delta}_{(\lambda_r),r}$ and we define $A_{r+1} = \left\{(1),\dots,(\lambda_{r+1})\right\}$. In case of ties, i.e. if $\widehat{\Delta}_{(\lambda_{r+1}),r} = \widehat{\Delta}_{(\lambda_{r+1} + m),r}$ for some $m$, we first add arms in $S_r$ to $A_{r+1}$. We emphasize that this tie-breaking rule is crucial in our analysis.
Arms in $A_{r+1}\backslash A_r$ are further classified as optimal (and added to the set $B_{r+1}$) or sub-optimal (and added to the set $D_{r+1}$) based on whether or not they belong to $S_r$. The output of the algorithm is the set $B_{R+1}\cup A_{R+1}$.

§.§ Particular Instances

Arm elimination algorithms have been proposed for different (mostly uni-dimensional) fixed-budget identification tasks <cit.> with different elimination rules, that could also be rewritten featuring some (simpler) gaps. In these works, two different types of arm schedule and sampling allocations have been mostly investigated:

* In Successive Rejects (SR), one arm is de-activated in each round, ie. $R=K-1$ and $\lambda^{\text{SR}}_{r}=K-r+1$ for all $r \leq K$. The sampling allocation proposed by  <cit.> satisfies $t^{\text{SR}}_r = n^{\text{SR}}_{r} - n^{\text{SR}}_{r-1}$ with $n_r = \left\lceil  \frac{1}{\overline{\log (K)}} \frac{T-K}{K+1-r}\right\rceil$ where $\overline{\log (K)} := 2^{-1} +\sum_{i=2}^{K}i^{-1}$ and $n_0=0$.

* In Sequential Halving (SH), one half of the active arm is de-activated in each round, that is $R=\lceil\log_2(K)\rceil$ and for all $r\in \{1,\dots, \lceil\log_2(K)\rceil\}$, $ \lambda^{\text{SH}}_{r+1} := \lceil \lambda_r^{\text{SH}}/2\rceil$ (we easily very that $\lambda^{\text{SH}}_{R +1} = 1$). The sampling allocation proposed by <cit.> is uniformly spread across rounds, that is $t_r^{\text{SH}} := \left\lfloor \frac{T}{\lvert A_r \lvert \lceil\log_2(K)\rceil } \right\rfloor$.

We refer to EGE-SR (resp. EGE-SH) as the instances of EGE using the same allocation as SR (resp. SH). In Appendix <ref> we propose a third instantiation using the geometric allocation of <cit.>.

For $D=1$, the PSI problem coincides with BAI and EGE-SR (resp. EGE-SH) coincides with SR (resp. SH[Besides the fact that the original SH algorithm for BAI discards samples collected in previous rounds to compute the empirical means.]). Indeed, in that case $S_r = \{\widehat a_r\}$ reduces to the empirical best arm, $\widehat \Delta_{i, r} = \muh_{\widehat a_r, r} - \muh_{i, r}$ for $i\neq \hat{a}_r$, $\widehat \Delta_{\widehat a_r, r} = \min_{i\in A_r\backslash\{\widehat a_r\}} \widehat \Delta_{i, r}.$
Then, our tie-breaking rule ensures that no arm is accepted before the last round and at each round $r$, $A_{r+1}$ is defined as the $\lambda_{r+1}$ arms in $A_r$ with the largest empirical means and the final survival arm is recommended as optimal.

§.§ Alternative Approach

Another idea to tackle fixed-budget PSI is to adapt an existing fixed-confidence algorithm to that setting. The APE algorithm of <cit.> takes  $\delta \in (0, 1)$ as input and uses a stopping time $\tau_\delta$ such that
\[ \bP\left(\cS^\star =  S(\tau_\delta) \text{ and } \tau_\delta \leq C H(\nu)f(\delta)\right) \geq 1 - \delta, \]

where $S(\tau_\delta)$ is the empirical Pareto set at time $\tau_\delta$ and $f$ is a function of $\delta$ and $C$ is a constant. The idea is then to choose $\delta$ so as $C H(\nu)f(\delta) \leq T$, that is tune $\delta$ w.r.t $H(\nu)$ and $T$. This is roughly the approach used
by UGapEb <cit.> in BAI. In Appendix <ref>, we analyze , a fixed-budget version of APE that takes as input a parameter $a\geq 0$ and we prove an upper-bound on its probability of error when $a\leq \frac{25}{36}\frac{T-K}{H(\nu)}$. However, such an “oracle” tuning is not very satisfying, as assuming $H(\nu)$ to be known is quite unrealistic in practice. This is why in the sequel we focus on presenting the analysis of which does not require any prior knowledge about $\nu$.

§ THEORETICAL GUARANTEES

We first propose an analysis of for a generic number of rounds $R$, arm schedule $\bm\lambda$ and sampling allocation $\bm t$ satisfying (<ref>) and (<ref>). It features the quantity

\[\widetilde T^{R,\bm t, \bm \lambda}(\nu)  := \min_{r\in [R]} \left(\sum_{s=1}^{r}t_s\right)\Delta_{(\lambda_{r+1} +1)}^2,\]

in which the dependency in $\nu$ is captured in the gaps.

Let $T\geq K$ and $\nu$ be a bandit with marginally $\sigma$-subgaussian arms. Then
Empirical Gap Elimination with parameters $R, \bm\lambda$ and $\bm t$ satisfies

$$ e^{\text{EGE}}_T(\nu) \leq 2(K-1) \lvert \cS^\star\lvert R D  \exp\left( - \frac{\widetilde T^{R,\bm t, \bm\lambda}(\nu)}{144\sigma^2} \right).$$

This result shows that the probability of failure of decreases exponentially fast with $\widetilde T^{R,\bm t, \bm \lambda}(\nu)$.
In Appendix <ref>, we further show that for both EGE-SR and EGE-SH, $\widetilde T^{R,\bm t, \bm \lambda}(\nu)$ is of order $T/H_2(\nu)$ with

\[   H_2(\nu) := \max_{i \in [K]} i \Delta_{(i)}^{-2},\]

where $(\cdot)$ is a permutation such that $\Delta_{(1)}\leq \dots\leq \Delta_{(K)}$. More precisely, we obtain the following.

Let  $T\geq K$ and $\nu$ be a bandit with $\sigma$-subgaussian marginals. Then EGE-SR satisfies

$$ e_T^\text{SR}(\nu) \leq  2(K-1)^2 \lvert \cS^\star \lvert D\exp\left( - \frac{T-K}{144\sigma^2H_2(\nu)\overline\log(K)}\right),$$

and for EGE-SH, $e_T^\text{SH}(\nu)$ is upper-bounded by

$$2(K-1)\lceil\log_2(K) \rceil \lvert \cS^\star \lvert  D\exp\left( - \frac{T}{288\sigma^2 H_2(\nu) \lceil\log_2(K)\rceil}\right).$$

The complexity measure $H_2(\nu)$ featured in our error exponent satisfies $H_2(\nu) \leq H(\nu) \leq H_2(\nu)\log(2K)$ as proved by <cit.>. For BAI ($D=1$), we essentially recover the existing guarantees for SR and SH, whose error bounds also feature $H_2(\nu)$, up to constant factors inside the exponential and an extra multiplicative $K$ factor for Sequential Halving. Still, to our knowledge this is the first analysis of the variant of SH that does not discard samples between rounds, which often performs (much) better in practice.

In the general case ($D\geq 1$) we remark that the bounds obtained for EGE-SR and EGE-SH are hard to compare: the latter has an improved polynomial dependence ($K\log_2(K)$ instead of $K^2$) but a worse constant inside the exponential. As we shall see in the experiments, both algorithms have actually pretty close performance (like in BAI, see <cit.>). Moreover, they both outperform a simple baseline using Uniform Allocation (UA) and recommending the Pareto set of the empirical means. We remark that Theorem <ref> yields an upper bound on the error probability of this strategy (by choosing $R=1$, $t_1 =T/K$ and $\lambda_2=0$, for which $\widetilde{T}^{1,\bm t,\bm\lambda}(\nu) = n_1\Delta_{(\lambda_2 +1)}^2 = T/(K\Delta_{(1)}^{-2})$), which we add to our summary in Table <ref>. This bound can be much worse than that for EGE-SR/SH when the gaps are distinct.

Algorithm     Error probability

EGE-SR    $K^2\lvert\cS^\star\lvert D\!\exp\!\left(\!-T/(H_2(\!\nu\!)\log K)\right)$
EGE-SH    $ K\!\log(\!K\!)\lvert\cS^{\!\star}\lvert D\!\exp\!\left(\!-T/(2H_2(\!\nu\!)\log\!K)\!\right)$
APE-FB$^\star$     $K\log(T)D\!\exp\!\left(\!-T/H(\!\nu\!)\right)$

UA     $K\lvert\cS^\star\lvert D\!\exp\!\left(\!-T/(K \Delta_{(1)}^{-2})\right)$

Upper bounds on $e_T(\nu)$ for different algorithms (up to constants). $^\star$tuned with $H(\nu)$.

§.§ Lower Bound

§.§ Sketch of proof of Theorem <ref>

We define for any arms $i, j$ and round $r$ the events

\begin{eqnarray*}
\xi_{i,j, r} \!\!\!\!&:=&\!\!\!\! \left\{\left \|(\vmuh_{i, n_r} \!- \vmuh_{j, n_r}) - (\vmu_i \!-\vmu_j)  \right\|_{\infty}\!\!\! \leq c \Delta_{(\lambda_{r+1} +1)}\right\} \\
\cE_c^1 \!&:=&\! \bigcap_{r\in [R]}\bigcap_{i \in \cS^\star}\bigcap_{j \in [K]}\;  \xi_{i, j, r}, \text{ for any } c>0.
\end{eqnarray*}

We shall prove that there exists some $c>0$ such that does not make any error on the event $\cE_c^1$. That is, no sub-optimal arm is added to $B_r$ and no  optimal arm is added to $D_r$, in any round $r$, and the possibly remaining arm in $A_{R+1}$ is an optimal arm.

To do so, an important step is to justify that any sub-optimal arm should be de-activated before the optimal arm that dominates it the most. More formally, for any sub-optimal arm $i$, we let $i^\star \in \argmax_{j\in \cS^\star} \m(i, j)$, which by definition is such that $\Delta_i = \m(i, i^\star)$. For a sub-optimal arm $i$, we know that $i^\star \in \cS^\star$ always exists. More importantly, $i^\star$ could be the only arm dominating $i$.
Therefore it is crucial to ensure that $i$ is no longer active before discarding $i^\star$, otherwise $i$ could appear as optimal w.r.t the remaining active arms.

Defining $\cP_r := \{ \forall i \notin \cS_\star, i \in A_r \Rightarrow i^\star \in A_r\}$, we first prove the following concentration result.

Assume that $\cE_c^1$ holds. Let $r\in [R]$ such that $\cP_r$ holds. Then, for any sub-optimal arm $i\in A_r$,

$$\lvert \widehat \Delta_{i, r}^\star - \Delta_i^\star \lvert \leq 2c\Delta_{(\lambda_{r+1}+1)} $$

and for any optimal arm $i \in A_r$,
$$ \widehat \delta_{i,r}^\star \geq \Delta_i -2c \Delta_{(\lambda_{r+1} +1)}.$$

This result then permits to prove by induction that $\cP_r$ holds in any round $r$, when $c$ is small enough.

lemmalemref Let $c < 1/6$.
On the event $\cE_c^1$, for any $r\in [R+1]$, $\cP_r$ holds.
In particular, for any sub-optimal arm $i$, $i^\star$ cannot be deactivated before $i$.

A first consequence is that if $A_{R+1}$ contains one arm (i.e. $\lambda_{R+1}=1$) and $\cE_c^1$ holds, then it is an optimal arm.

Note that in BAI ($D=1$), for any sub-optimal arm $i$, $i^\star$ is the unique best arm and <ref> is sufficient to ensure the correctness of EGE. But in the general case we need to ensure that no optimal arm is rejected and no sub-optimal arm is accepted. This is done by using <ref> and <ref>.
lemmanewlemma Let $c < 1/6$.
On the event $\cE_c^1$, $\widehat{S}_T  = \cS^\star$.

<ref> then follows by upper-bounding the probability of $\bar{\cE_c^1}$ using Hoeffding's inequality.

§ RELAXING PSI

In this section we explain how EGE-SR can be slighlty adapted to tackle the “at most $k$ optimal arms” relaxation (or PSI-$k$) first introduced by <cit.> in the fixed confidence setting.

In this problem the goal is to return a subset $\widehat{S}_T$ of $\cS^\star$ of size $k$ , or $\cS^\star$ itself if its size is smaller than $k$. In the fixed budget setting, we define the $k$-relaxed expected loss as $e_{T, k}(\nu) := \bE_\nu[\cL(\widehat S_T, k)]$ where

\begin{eqnarray*}
\cL(\widehat S, k) := \begin{cases}
\ind\{\widehat S \subset \cS^\star\} & \text{if} \quad  \lvert \widehat S \lvert = k, \\
\ind\{\widehat S = \cS^\star\} & \text{else}.
\end{cases}
\end{eqnarray*}

To minimize $e_{T, k}(\nu)$ for any parameter $k$ and budget $T$, we propose , a variant of EGE-SR which may stop at some round $r < K-1$ if $\lvert B_{r+1}\lvert = k$ and recommend $B_{r+1}$ (see ).

In Appendix <ref>, we prove an upper-bound on the expected loss $e_{T, k}(\nu)$. To introduce it, we define $\omega_{(k)}$ to be the $k$-th largest gap among the optimal arms :
$ \omega_{(k)}:= \max_{i\in \cS^\star}^{k} \Delta_i$ with $\omega_{(k)} = 0$ if $|\cS^\star| < k$. Our bound features the complexity measure $H_2^{(k)}(\nu) :=\max_{i\in [K]} i(\Delta^{(k)}_{(i)})^{-2}$ with the $k$-relaxed gaps

\begin{equation}
\Delta_i^{(k)} :=
\begin{cases}
\max(\Delta_i, \omega_{(k)}) & \text{if} \quad i \in \cS^\star \\
\Delta_i &\text{else.}
\end{cases}
\end{equation}

At most $k$ Pareto optimal arms
$k\in [K]$, Set $A_1= [K]$, $B_1 = \emptyset$, $D_1 = \emptyset$

$r=1, 2, \dots, K-1$

Collect $t_r$ samples from each arm $a \in A_r$

Compute $S_r$ the empirical Pareto set

Choose $i_r \in \argmax_{i\in A_r} \widehat \Delta_{i, r}$  ties broken in favor of arms in $S_r$

$A_{r+1} \gets A_r \backslash\{ i_r\}$

$B_{r+1} \gets B_r \cup \{ S_r \cap \{i_r\}\}$

$D_{r+1} \gets D_r \cup \{\{i_r\} \cap (A_r\backslash S_r) \}$

$\lvert B_{r+1} \lvert=k$
break and return  $\widehat S_T = B_{r+1}$

$\widehat S_T = B_{K}  \cup A_K$

Let $k\in [K]$. satisfies
$$e_{T, k}(\nu) \!\leq 2(K-1)^{2}\lvert \cS^\star \!\lvert D\exp\!\left(\!- \frac{T-K}{144\sigma^2H_2^{(k)}(\nu)\overline\log(K)}\!\right)\!.$$

This result is particularly insightful when there are many optimal arms and some of them are easy to identify as such (large gaps).
Indeed when $\lvert \cS^\star \lvert \approx K$ and $k\ll \lvert \cS^\star\lvert$, $H_2^{(k)}(\nu)$ can be an order of magnitude smaller than $H_2(\nu)$.
We also note that when  $k>\lvert \cS^\star \lvert$ (and PSI-$k$ reduces to PSI), $H_2^{(k)}(\nu) = H_2(\nu)$ and we recover the result of <ref>.

The stopping time of is
$$ \tau := \inf\left\{r: \lvert B_{r+1}\lvert  = k\right\} \land (K-1).$$

Letting $N_{\tau}$ denote the total number of samples used at termination, we upper-bound $\bE_\nu[\tau]$ and $\bE_\nu[N_\tau]$ in Appendix <ref>, showing that when the budget is large, $\bE_\nu[\tau]\leq q$ and $\bE_\nu[N_\tau] \leq N_q$, with $q:=K-\lvert \cS^\star \lvert + k$. Intuitively it suggests that in the worst case the $(K-\lvert \cS^\star \lvert)$ sub-optimal arms are discarded before $k$ optimal arms are accepted as they might be needed to dominate some sub-optimal arms. Moreover, we remark that $q$ can be way smaller than $(K-1)$. For instance when $[K]= \cS^\star$, we have $q=k$.

§ EXPERIMENTAL STUDY

We evaluate our algorithms on synthetic and real-world tasks.
We compare to Uniform Allocation (UA) and for three parameters
of the form $a_c = c \frac{25}{36}\frac{T -K}{H(\nu)}$, $c\in \{1/10, 1, 10\}$. Our guarantees on are only valid for $c\leq 1$ and are optimal for $c=1$. We consider a heuristic version of which adaptively estimates the hardness $H(\nu)$. We refer to this algorithm as APE-FB-ADAPT. We run the experiments 4000 times with different seeds and we report the $\log_{10}$ of the average mis-identification rate.

Application 1: COV-BOOST trial

Application 2: Sorting Networks dataset.

Arms on a convex Pareto set.

Each sub-optimal $i$ is only dominated by $i^\star$.

$K=200$ arms on the unit circle.

High dimension ($D=10$) with 2 group of arms.

§.§ Real-world Datasets

COV-BOOST <cit.> is a phase II vaccine clinical trial conducted on 2883 participants to measure the immunogenicity of different Covid-19 vaccines as a booster (third dose). Combining the first two doses received and the third dose investigated in the trial, there were $K=20$ arms and the authors reported the sample mean response and confidence intervals (based on log-normal assumption of the data) of each arm to a bunch of immunogenicity indicators. <cit.> further extracted and processed the responses to 3 indicators (cellular response, anti-spike IgG and NT$_{50}$) to generate a multivariate normal bandit with $K=20, D=3$ and a diagonal covariance matrix. We use their dataset in our simulations with a total budget of $T=2500$.

Hardware Design
We use the SNW (Sorting Network Width) dataset <cit.> to generate a bandit model. The dataset is made of 206 different sorting networks (devices used to sort items). The authors reported  area and throughput of each network when synthetized on a FPGA (Field-Programmable Gate Array). The area is the number of FPGA slices (resources units) used during execution and the throughput is the number of samples treated per second. The goal is then to optimize both objectives simultaneously (reduced area and large throughput). The simulation process is costly and the measures may slightly vary for the same network due to randomness in the circuits. We simulate a Gaussian bandit with $K=206, D=2$ with the SNW dataset. We use a total budget of $T=5000$.

§.§ Synthetic Benchmark

We run the algorithms on different synthetic instances. For each of them, we compute the complexity $H(\nu)$ and following <cit.> we set the total budget to $T= H(\nu)$.

Experiment 1: Arms on a convex Pareto set.
$K=60, D=2$. We choose $x_1,\dots, x_{20}$ equally spaced in $[0.55, 0.95]$ and $i=1,\dots, 10$, $\vmu_i:=(x_i^2, 1/(4x_i^2))^\T$. $\vmu_{11}, \dots, \vmu_{60}$ are chosen from $\left\{(x, y)\in [0.1, 0.8]^2: xy\leq \frac15\right\}$.

Experiment 2: Group of sub-optimal arms where each sub-optimal arm $i$ has a unique $i^\star$. $K=10, D=2, \lvert \cS^\star\lvert = 2$. For each sub-optimal arm $i$, there is a unique $j$ such that $\vmu_i \prec \vmu_j$. We choose $\vmu_1 := (0.4, 0.75)^\T, \vmu_2 := (0.75, 0.4)^\T$  and for $i=1, \dots, 4$, $\vmu_{2i+1} := (0.45 + 0.2^i, 0.35 - 0.2^i)^\T$, $\vmu_{2i+2}:=(0.10 + 0.20^i, 0.70 - 0.20^i)^\T$.

Experiment 3: Many arms on the unit circle. $K=200, D=2$ and generate n isotropic multivariate normal instance with $\sigma= 1/4$.  We choose $\beta_1, \dots, \beta_{20}$ equally in $[\pi/12, \pi/2 - \pi/12]$ and $\beta_{21},\dots, \beta_{200}$ equally spread in $[\pi/2+\pi/6, 2\pi - \pi/6]$.  For each arm $i, \vmu_i:= (\cos(\beta_i), \sin(\beta_i))^\T$.

Experiment 4: $D=10$ We set $K=50$ and generate $(\vmu_1, \dots, \vmu_{30}) \sim \mathcal{U}\left([0.2, 0.45]^{10}\right)^{\otimes 30}$ and $(\vmu_{31}, \dots, \vmu_{50}) \sim \mathcal{U}\left([0.55, 0.75]^{10}\right)^{\otimes 20}$

§.§ Results

In all experiments, the uniform allocation (UA) baseline is largely outperformed by both EGE-SR and EGE-SH (with no clear ordering between the two algorithms). This is particularly the case when there are many arms and for complex Pareto sets. By estimating the hardness to classify each arm, eventually allocates more samples to arms that are difficult to classify, leading to a smaller error probability.
is a good competitor to our algorithms however it is not robust to the hyper-parameter $a$, which requires the knowledge of $H(\nu)$ to be properly tuned. Our proposed heuristic that estimates the complexity online fails dramatically.

In Appendix <ref> we discuss the computational complexity of and we detail the implementation setup. In Appendix <ref> we report additional experimental results on synthetic datasets. We illustrate the PSI-$k$ relaxation on the SNW dataset, showing that can efficiently identify a subset of the Pareto set.

§ CONCLUSION

We proposed the first algorithms for Pareto set identification in the fixed-budget setting.

Our generic algorithm can be coupled with different allocation and arm elimination schemes. We proved that for two instantiations, EGE-SR and EGE-SH, the probability of error decays exponentially fast with the budget, with an exponent that is un-improvable for some bandit instances. We conducted experiments showing that these algorithms consistently outperform a uniform sampling baseline and are competitive with an oracle algorithm that knows the complexity of the underlying instance.

Extensions of this work could consider variance-aware estimates of the gaps or a Bayesian setting where prior information on the means could be used to improve performance <cit.>. Finally, just like Successive Rejects or Sequential Halving, EGE algorithms heavily relies on the knowledge of $T$ for their tuning. In future work, we are interested in proposing anytime algorithms, that can have a small error probability for any budget $T$. We may take inspiration from some anytime algorithms recently proposed for simpler pure exploration tasks <cit.>.

§.§.§ Acknowledgements

Cyrille Kone is funded by an Inria/Inserm PhD grant. Emilie Kaufmann acknowledges the support of the French National Research Agency under the projects BOLD (ANR-19-CE23-0026-04) and FATE (ANR22-CE23-0016-01).


# Contribution of spicules to solar coronal emission

Shanwlee Sow Mondal Astronomy and Astrophysics Division, Physical Research
Laboratory, Ahmedabad 380009, India  Indian Institute of Technology,
Gandhinagar, Gujarat 382355, India  James A. Klimchuk Heliophysics Science
Division, NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD
20771, USA  Aveek Sarkar Astronomy and Astrophysics Division, Physical
Research Laboratory, Ahmedabad 380009, India

###### Abstract

Recent high-resolution imaging and spectroscopic observations have generated
renewed interest in spicules’ role in explaining the hot corona. Some studies
suggest that some spicules, often classified as type II, may provide
significant mass and energy to the corona. Here we use numerical simulations
to investigate whether such spicules can produce the observed coronal emission
without any additional coronal heating agent. Model spicules consisting of a
cold body and hot tip are injected into the base of a warm ($0.5$ MK)
equilibrium loop with different tip temperatures and injection velocities.
Both piston- and pressure-driven shocks are produced. We find that the hot tip
cools rapidly and disappears from coronal emission lines such as Fe XII $195$
and Fe XIV $274$. Prolonged hot emission is produced by pre-existing loop
material heated by the shock and by thermal conduction from the shock.
However, the shapes and Doppler shifts of synthetic line profiles show
significant discrepancies with observations. Furthermore, spatially and
temporally averaged intensities are extremely low, suggesting that if the
observed intensities from the quiet Sun and active regions were solely due to
type II spicules, one to several orders of magnitude more spicules would be
required than have been reported in the literature. This conclusion applies
strictly to the ejected spicular material. We make no claims about emissions
connected with waves or coronal currents that may be generated during the
ejection process and heat the surrounding area.

methods: numerical, Sun: corona, Sun: chromosphere, Sun: atmosphere, Sun:
magnetic fields, Sun: UV radiation

##  1 Introduction

Defying decades of continued efforts, many aspects of coronal heating remain
unanswered (Klimchuk, 2006, 2015; Viall et al., 2021). Even the basic
mechanism is a matter of debate. Despite the fact that all the coronal mass is
sourced at the chromosphere, agreement on how the chromospheric mass is heated
and transported up to the corona has not been reached. An early observation of
the solar chromosphere reported the existence of several small jet-like
features (Secchi, 1877). They were later named spicules (Roberts, 1945). With
improved observations, these spicules were seen to propagate upwards (Pneuman
& Kopp, 1977, 1978) with speed $20$ \- $50$ km s-1. They were also seen to
survive for about $5$ to $10$ minutes and carry almost $100$ times the mass
needed to balance the mass loss in the solar corona due to the solar wind.
Further studies of the spicules (Athay & Holzer, 1982) suggested a pivotal
role in transferring energy from the inner layers of the solar atmosphere to
the lower solar corona. However, the proposal was not pursued further because
these traditional spicules lack emission in the Transition Region (TR) and
coronal lines (Withbroe, 1983).

About a decade ago, using high-resolution imaging and spectroscopic
observations from the Hinode and Solar Dynamic Observatory missions, De
Pontieu et al. (2007, 2011) further discovered jet like features traveling
from the chromosphere to the corona. These features appear all over the Sun
with a lifetime between $10-150$ s and a velocity of $50-150$ km s-1. De
Pontieu et al. (2007) termed them type II spicules and suggested that they are
capable of connecting the relatively cooler solar chromosphere with the hot
corona.

Since their discovery, multiple observations have identified type II spicules
and reported on their characteristics. However, nothing conclusive has yet
been established about their origin. Only recently, Samanta et al. (2019) have
identified the near-simultaneous origin of spicules and the emergence of
photospheric magnetic bipoles. The tips of the originated spicules eventually
appear in the coronal passband, suggesting that the plasma is heated to
coronal temperatures. A 2.5D radiative MHD simulation of type II spicules
(Martínez-Sykora et al., 2017) has reproduced many of their observed features.
This simulation also suggests that ambipolar diffusion in the partially
ionized chromosphere may play a crucial role in the origin of type II
spicules. On the other hand, a recent work (Dey et al., 2022) based on
radiative MHD simulation and laboratory experiment suggests that quasi-
periodic photospheric driving in the presence of vertical magnetic fields can
readily generate spicules in the solar atmosphere. Their work, devoid of any
chromospheric physics, can still account for the abundance of wide varieties
of spicules, as seen in the observations.

The evolution of spicules during their propagation is understood through
multi-wavelength studies (e.g., De Pontieu et al., 2011; Skogsrud et al.,
2015). Observations of De Pontieu et al. (2011) suggest that spicule plasma
emanating from the chromosphere gets heated to transition region (TR)
temperatures and even up to coronal temperatures. Such heating may happen for
two reasons:

  1. (a)

Spicule propagation can produce shocks, compressing the material lying ahead
of it. In such a scenario, it is not the ejected spicule material but the pre-
existing coronal material in front of it that gets compressed by the shock to
contribute to the hot emission (Klimchuk, 2012; Petralia et al., 2014);

  2. (b)

The tip of the spicule may get heated during the ejection process, on-site,
through impulsive heating and produce emissions in the coronal lines. In the
latter scenario, the emission indeed comes from the ejected spicule material
(De Pontieu et al., 2007).

The radiative MHD simulations of Martínez-Sykora et al. (2018) suggest that
spicules and surrounding areas get heated by ohmic dissipation of newly
created currents and by waves. Note, however, that the currents in the
simulations are relatively large-scale volume currents and would not be
dissipated efficiently at the many orders of magnitude smaller resistivity of
the real corona. Heating in the real corona involves magnetic reconnection at
thin current sheets, of which there are at least $100,000$ in a single active
region (Klimchuk, 2015). It is not known whether the ohmic heating in the
simulations is a good proxy for the actual reconnection-based heating.

Klimchuk (2012) considered a simple analytical model for the evolution of
spicules with a hot tip. He argued that if a majority of observed coronal
emission were from such hot tips, it would be inconsistent with several
observational features (see also Tripathi & Klimchuk (2013); Patsourakos et
al. (2014)). The result was further supported by hydrodynamic simulations
(Klimchuk & Bradshaw, 2014). Using these simulations, they studied the
response of a static loop to impulsive heating in the upper chromosphere,
which produces localized hot material that rapidly expands upward and might
represent the hot tip of a spicule. Noticing the inability of a single hot
spicule tip to explain the observations, Bradshaw & Klimchuk (2015) further
explored the role of frequently recurring chromospheric nanoflares. The study
was motivated by the suggestion that rapidly repeating type II spicules might
accumulate enough hot plasma to explain the coronal observations (De Pontieu
et al., 2011). However, the simulations were still inconsistent with
observations.

In both the analytical model and the simulations, the dynamics of the hot
material is due entirely to an explosive expansion from the locally enhanced
pressure. There is no additional imposed force to bodily eject the material.
The consequences of such a force were investigated by Petralia et al. (2014).
Their study involves injecting cold and dense chromospheric material into the
corona with an initial velocity. The result indicates that the production of a
shock can give rise to coronal emission. However, the emission is from the
preexisting coronal material rather than the spicule itself. The injected
material has no hot component.

The studies mentioned above have investigated the dynamics of either the hot
tip of a spicule without any initial velocity or a spicule with a cold tip and
finite injection velocity. Our work combines these two effects. The spicule is
now injected in a stratified flux tube with high velocity and consists of both
a hot tip and a cold body ($T=2\times 10^{4}$ K). We further investigate the
possibility that most of the observed hot emission from the corona can be
explained by such spicules. Through forward modelling, we quantitatively
compare the simulations with observations to answer this question.

The rest of this paper is organized as follows. The numerical setup is
described in Section 2. We report on the simulation results in Section 3.
Finally we summarize and discuss our results in Section 4.

##  2 Numerical Setup

Spicules are seen to follow magnetic field lines. To simulate their dynamics,
we consider a straight 2D magnetic flux tube consisting of uniform $10$ G
magnetic field. We impose a gravity corresponding to a semi-circular loop such
that the vertical component of the gravitational force is maximum at both ends
and smoothly becomes zero in the middle of the tube. Two ends of the tube are
embedded in the chromosphere. The loop is symmetric about the center, which
corresponds to the apex. We use Cartesian coordinates, and therefore the loop
actually corresponds to an infinite slab. This is a reasonable approximation
because we are interested in how the plasma evolves within an effectively
rigid magnetic field appropriate to the low $\beta$ corona. The slab dimension
corresponding to the loop length is $100$ Mm. The other dimension is $0.42$
Mm, but this is not relevant. Rigid wall boundary conditions are imposed at
the sides, and the evolution is essentially equivalent to 1D hydrodynamics, as
discussed below. The first 2 Mm of both ends of the loop are resolved with a
fine uniform grid with 10 km cells, while the coronal part is resolved with a
stretched grid containing 1500 cells on each side. The fine grid close to both
the footpoints allows us to resolve the steep transition region more
accurately.

The spicule simulation begins with an initial static equilibrium atmosphere
obtained with the double relaxation method described in Appendix A. We choose
a relatively low temperature and low density loop because we wish to test the
hypothesis that the observed coronal emission comes primarily from spicules.
The apex temperature of the loop is $0.5$ MK. Figure 1 shows the background
loop profile that is used in most of our simulations. The chromosphere is
$470$ km deep - approximately half a gravitational scale height. It merely
acts as a mass reservoir. Detailed chromospheric physics like partial
ionization and optically thick radiation are not implemented in the code as we
are solely interested in coronal emission. We use a modified radiation loss
function to maintain a chromospheric temperature near $2\times 10^{4}$ K, as
described in Appendix A.

The propagation of a spicule in the loop is emulated through an injection of
dense material from the left footpoint. The injected material follows
specified density, velocity and temperature profiles in time which are
described below. At this injection boundary, all plasma parameters, except the
density and pressure, are set to their initial values once the injection is
over. The density is set to have the prescribed value at the end of the
injection phase, and pressure is determined from the ideal gas equation of
state. On the other hand, at the right footpoint, all the plasma parameters
maintain the initially prescribed values throughout the entire simulation.

We solve the compressible MHD equations inside our simulation domain using the
PLUTO code (Mignone et al., 2007) with ideal gas environment. Plasma inside
the domain is cooled by radiation and field aligned thermal conduction. The
CHIANTI (Landi et al., 2013) radiative loss function for coronal abundance is
used to model the radiative cooling. For anisotropic conduction, the thermal
conductivity, $\kappa_{\parallel}$ = $5.6\times 10^{-7}T^{5/2}$ erg s-1 K-1
cm-1 is considered along the magnetic field lines, whereas $\kappa_{\perp}$ is
taken to be zero. Also, for the saturated conductive flux used in PLUTO,
$F_{sat}$ = 5$\phi\rho C_{s}^{3}$, where we have considered the value of the
free parameter $\phi$ to be 0.9, which represents effective thermal conduction
in the system, and $C_{s}$ is the isothermal sound speed. The MHD equations
are solved in Cartesian coordinates.

The photospheric magnetic flux is found to be localized and clumpy, whereas,
in the corona, it fills out space uniformly. Such nature of the magnetic flux
at different layers of the solar atmosphere dictates that the flux tubes
expand laterally at the junction of the chromosphere and corona, where the
plasma $\beta$ changes from being greater than one to less than one. This type
of expansion of flux tubes is realized in 2D MHD simulations of coronal loops
(e.g., Guarrasi et al., 2014). Through an area expansion factor, this has also
been incorporated in 1D or 0D models (Mikić et al., 2013; Cargill et al.,
2022). We do not include expansion in our model because we are interested in
the spicule dynamics in the corona, and the simplification should not affect
our results significantly. We note that the plasma $\beta$ is less than unity
throughout the evolution, so no expansion from the spicule injection would be
expected. Additionally, the initial atmosphere and injection profile are
uniform along the horizontal (cross-field) axis. Hence the plasma remains
nearly uniform in the lateral direction, effectively making our simulations
similar to 1D hydrodynamic simulations. Nevertheless, we ran all our
computations using the 2D MHD set up because of our familiarity with the
powerful PLUTO code. The limited number of grid points in the cross-field
direction keeps the computational demands relatively low.

Two main components of our simulations are: (a) a background loop in
hydrostatic and energy equilibrium representing a tenuous coronal atmosphere,
and (b) the propagation of injected material resembling spicule propagation
along the loop. Our experimental spicule consists of a hot dense tip followed
by cold dense material injected from the base of the model. Here we
investigate how changing the temperature of the hot tip and injection speed
can alter the intensities and profiles of the Fe XII (195 Å) and Fe XIV(274 Å)
coronal spectral lines.

We have performed six sets of simulations where the spicule tip temperatures
are considered to be at $2$, $1$, and $0.02$ MK, followed by cold material
with a temperature of $0.02$ MK. All the runs are performed with two injection
velocities: $50$ and $150$ km s-1 (see Table 1). Since we assume that spicules
might have been generated deep inside the chromosphere, we inject a high-
density material in the loop to emulate the spicule. The density scale height
of the spicule is chosen to be six times the gravitational scale height at the
base of the equilibrium loop. To impose such conditions on the ejected
spicule, its density follows a time profile given by,

$\rho(t)=\begin{cases}\rho_{0}\exp\Big{[}\frac{v(t)t}{6H}\Big{]},&\ 0<t\leq
t_{5}\\\ \rho(t_{5}),&\ t_{5}<t\\\ \end{cases}~{},$ (1)

where $\rho(t)$ and $\rho_{0}$ are the injected density at time $t$ and the
base density of the equilibrium loop, respectively. The time profile of the
injection velocity is given by,

$v(t)=\begin{cases}v_{inj}\times\Big{(}\frac{t}{t_{1}}\Big{)},&\ 0<t\leq
t_{1}\\\ v_{inj},&\ t_{1}<t\leq t_{4}\\\
v_{inj}\times\Big{(}\frac{t_{5}-t}{t_{5}-t_{4}}\Big{)},&\ t_{4}<t\leq t_{5}\\\
0,&\ t_{5}<t\\\ \end{cases}~{},$ (2)

where $v_{inj}$ corresponds to $50$ or $150$ km s-1 (depending on the
simulation). $H$ represents the gravitational scale height given by

$H=\frac{k_{B}T_{base}}{\mu m_{H}g_{\odot}}~{},$ (3)

where $T_{base}=0.02$ MK is the base temperature of the loop, $k_{B}$ is the
Boltzmann constant, while $m_{H}$ and $g_{\odot}$ represent mass of the
hydrogen atom and solar surface gravity, respectively, and $\mu=0.67$ denotes
the mean molecular weight of the plasma. The temperature of the ejected
spicule also follows a time profile given by

$T(t)=\begin{cases}T_{base}+(T_{tip}-T_{base})\times\Big{(}\frac{t}{t_{1}}\Big{)},&\
0<t\leq t_{1}\\\ T_{tip},&\ t_{1}<t\leq t_{2}\\\
T_{tip}+(T_{base}-T_{tip})\times\Big{(}\frac{t-t_{2}}{t_{3}-t_{2}}\Big{)},&\
t_{2}<t\leq t_{3}\\\ T_{base},&\ t_{3}<t\\\ \end{cases}~{},$ (4)

where $T_{base}$ is the temperature of the cold material (bottom part) of the
spicule ($=0.02$ MK) and $T_{tip}$ is the spicule tip temperature which can
take values $2$, $1$, or $0.02$ MK depending on the run being performed. In
all the above equations, times $t_{1}$, $t_{2}$, $t_{3}$, $t_{4}$ and $t_{5}$
are chosen to be $2$, $10$, $12$, $90$ and $100$ s, respectively. Times are
chosen so that the top $10\%$ of the spicule’s body emits in coronal lines as
is generally observed (De Pontieu et al., 2011). The total injection duration
is also motivated by the observed lifetime of type II spicules (De Pontieu et
al., 2011). The ramping up of velocity, density, and temperature ensures a
smooth entry of the spicules into the simulation domain. Similarly, the
ramping down at the end of the injection avoids any spurious effects. Figure 2
shows one such example of velocity, density, and temperature profiles when the
spicule is ejected with velocity $150$ km s-1, and its hot tip is at $2$ MK.
Likewise, different injection time profiles have been used for other injection
velocities and temperatures. The initial equilibrium loop remains the same in
all cases, unless specified.

Figure 1: Density and temperature profiles of the initial static equilibrium
loop. Figure 2: Injection velocity, density and temperature profiles as a
function of time at the bottom boundary from Run1.

##  3 Results

The large velocity of the spicule and high pressure compared to the ambient
medium give rise to a shock, which propagates along the loop and heats the
material ahead of it. Depending on the sound speed of the ambient medium
(i.e., the preexisting loop plasma) and the temperature of the injected
spicule material, the generated shock turns out to be of different kinds: (a)
Piston driven shock – in which case the shock speed is nearly equal to the
injection speed (e.g., simulation with $T_{tip}=0.02$ MK), and (b) Pressure
driven shock – in which case the shock speed exceeds the injection speed
(e.g., when $T_{tip}=2$ and $1$ MK). Emission from the shock heated plasma
differs depending on the nature of the shock.

We compare different simulations to understand the coronal response to
spicules with different injection parameters. Our discussion starts with the
results from Run1 where the hot tip of the injected spicule has a temperature
$T_{tip}=2$ MK and injection velocity $v=150$ km s-1. The injection profiles
are those already shown in Figure 2.

###  3.1 Dynamics and heating

Insertion of dense, high temperature plasma ($T_{tip}=2$ MK) with high
velocity ($v=150$ km s-1, Run1) into the warm loop produces a shock. Figure 3
shows the temperature, density and plasma velocity along the loop at $t=70$ s.
The dashed lines mark the location of the shock front. It is evident from the
figure that the high compression ratio exceeds the ratio of an adiabatic
shock. The compression ratio of an adiabatic shock should always be $\leq 4$.
To understand the nature of the shock, we perform a shock test with Rankine-
Hugoniot (RH) conditions, which read

$\frac{\rho_{2}}{\rho_{1}}=\frac{\gamma+1}{\frac{2}{M^{2}}+(\gamma-1)}=\frac{v_{1}}{v_{2}}~{}.$
(5)

Here $\rho_{1}$ and $\rho_{2}$ are the pre- and post-shock plasma mass
densities respectively, and $v_{1}$ and $v_{2}$ are likewise the pre- and
post-shock plasma velocities in the shock rest frame. Furthermore, $\gamma$ is
the ratio of the specific heats, $c_{s}=\sqrt{\frac{\gamma P_{1}}{\rho_{1}}}$
is the upstream sound speed, where $P_{1}$ is the upstream pressure, and
finally $M=v_{1}/c_{s}$ is the upstream Mach number in the shock reference
frame. Injection of high temperature plasma accelerates the shock with a speed
much larger than the injection speed of the spicule material giving rise to a
pressure driven shock front.

Figure 3: Plasma variables along the loop at $t=70$ s in Run1. The colored
dashed line marks the location of shock front. The observed density and
velocity jumps match well with $\gamma=1.015$ indicating the shock is locally
isothermal, consistent with the temperature profile.

Figure 3 demonstrates an abrupt change in plasma variables at the shock. The
shock speed at this instant is $562$ km s-1. It also shows that at the
discontinuity location ($s=36.9$ Mm, $s$ being the coordinate along the loop),
the density and velocity ratios are $10.7$ and $0.094$, respectively, in the
shock rest frame. The inverse relationship of these ratios indicates a
constant mass flux across the shock front, in accordance with equation (5).
The Mach number in the shock frame at the same location is $3.37$. With this
Mach number, the RH condition gives density and velocity ratios in accordance
with those in the simulation ($10.5$ and $0.095$) when $\gamma=1.015$. In
other words, consistency is achieved with this value of $\gamma$. Being close
to unity, it implies a nearly isothermal shock. Efficient thermal conduction
carries a large heat flux from the shock front to its surroundings, giving
rise to the locally smooth, near-isothermal temperature profile in Figure 3.
It is worth mentioning here that RH-jump conditions do not consider any heat
loss/gain, such as thermal conduction or radiative loss. However, our system
includes these sink terms in the energy equation. It is because of these loss
functions the shock-jump is larger. Limited thermal conduction would bring the
jump condition closer to the adiabatic approximation but would also affect the
thermal profile ahead and behind the shock. Our result is consistent with that
of Petralia et al. (2014), where the signature of the shocks in front of the
spicule has been reported. As we show later, the initially hot material in the
spicule tip cools dramatically. Only ambient material heated by the shock is
hot enough to produce significant coronal emission.

Interestingly, the high compression ratio at the shock front depends more on
the temperature difference and corresponding pressure difference between the
injected and ambient plasma material than on the velocity with which it is
injected. Table 1 shows a study of how the compression ratio (or the shock
strength) varies when the tip of the spicules are at different temperatures
and are injected with different velocities. As mentioned earlier, the
injection conditions give rise to two different types of shocks. When the
injected plasma temperature is high (e.g., spicule tips with temperatures $2$
and $1$ MK), the excess pressure gives rise to a pressure-driven shock. On the
other hand, injection of a cold material (tip temperature equal to that at the
loop footpoint, i.e., $0.02$ MK) produces a piston-driven shock. Our test runs
identify both kinds of shocks. For example, when we inject spicules with a
fixed injection velocity of $150$ km s-1, but with different tip temperatures
(viz. $2$, $1$ and $0.02$ MK), the average shock speed is $520$, $400$ and
$210$ km s-1, respectively (see Figure 11). The first two shocks are pressure-
driven as the average shock speeds exceed the injection speed by a wide
margin. The third shock maintains a speed close to the injection speed and can
be categorized as a piston-driven shock. The shock speed depends not only on
the injected tip temperature, but also on the properties of the ambient
material in which it is propagating, which vary along the loop. This is
discussed further in Appendix B.

Table 1: Dependence of compression ratio on the injected hot tip temperature and speed. Run |  $T_{tip}$ |  $v$ |  Compression
---|---|---|---
|  (MK) |  (km s-1) |  ratio
1 |  2 |  150 |  11.2
2 |  2 |  50 |  8.9
3 |  1 |  150 |  8.7
4 |  1 |  50 |  6.2
5 |  0.02 |  150 |  3.6
6 |  0.02 |  50 |  1.7

###  3.2 Loop emission

Thermally conducted energy from the shock front heats the material lying ahead
of it. Therefore, a magnetic flux tube subjected to spicule activity could
produce hot emission from newly ejected material at the spicule’s hot tip and
from pre-existing coronal material in both the pre and post-shock regions. We
now examine the contributions from these three different sources. We identify
the leading edge of the hot spicule tip by finding the location in the loop
where the column mass integrated from the right footpoint equals the initial
column mass of the loop. Recall that the spicule is injected from the left
footpoint. The spicule compresses the material in the loop, but does not
change its column mass. We identify the trailing edge of the hot material in a
similar manner, but using the column mass at time $t=10$ s, when the injection
of hot material ceases and the injection of cold material begins.

Figure 4 shows emission along the loop in the Fe XII and Fe XIV lines at
$t=10$ and $70$ s, evaluated from Run1. The orange region is the hot spicule
tip, while the red region is the shock-heated material ahead of it. The shock
front is the dot-dashed black vertical line. The dark orange curve is
temperature in units of $10^{5}$ K, with the scale on the left. The blue curve
is the logarithm of density, with the scale on the right. The yellow and green
curves are the logarithms of Fe XII and Fe XIV intensity, respectively, with
the scale on the left. The variation of intensity is enormous; a difference of
10 corresponds to 10 orders of magnitude. The intensity is what would be
observed by the Extreme ultraviolet Imaging Spectrometer (EIS; Culhane et al.
(2007)) onboard Hinode (Kosugi et al., 2007) if the emitting plasma had a
line-of-sight depth equal to the EIS pixel dimension, i.e., if observing an
EIS pixel cube. This can be interpreted as normalized emissivity.

At $t=10$ s, the emission in both lines comes primarily from the injected hot
plasma (orange region). On the other hand, at $t=70$ s it comes primarily from
the shock heated plasma (red region). The transition happens very early on.
Shortly after the injection of the hot material stops ($t=10$ s), emission
from the shock heated material starts dominating the total emission from the
loop. This is evident in the time evolution plot of the loop-integrated
emission in Figure 5. Shown are the intensities that would be observed by EIS,
assuming that the loop has a cross section equal to the pixel area and that
all of the loop plasma is contained within a single pixel. This corresponds to
a loop that has been straightened along the line of sight and crudely
represents a line of sight passing through an arcade of similar, out of phase
loops. The black curve shows the evolution of the total emission contributed
by the spicule and pre-existing plasma. Subtracting the spicule component (red
curve) from the total gives the evolution of the emission coming solely from
the pre-existing (non-spicule) loop material (green curve). Soon after the hot
tip of the spicule completes its entry into the loop (at $t=10$ s), the
emission from the spicule falls off rapidly. This is because the hot material
at the spicule tip cools rapidly as it expands in the absence of any external
heating. It is far too faint to make a significant contribution to the
observed coronal emission, as emphasized earlier by Klimchuk (2012) and
Klimchuk & Bradshaw (2014).

Figure 4: Logarithm of Fe XII and Fe XIV intensity (yellow and green,
respectively) and temperature (dark orange) as a function of position along
the loop at t = 10 and 70 s in Run1. Units are given in the legend. The blue
curve is the logarithm of density with the scale on the right. The marked
light orange region indicates the hot tip of the spicule, while the red region
indicates the compressed post-shock plasma. The dot-dashed line (ahead of the
red region) marks the location of the shock front.

For better comparison with observations, we construct synthetic spectral line
profiles. The methodology is explained in Appendix C. To construct these
profiles, we imagine that the loop lies in a vertical plane and is observed
from above. We account for the semi-circular shape when converting velocities
to Doppler shifts. We then integrate the emission over the entire loop and
distribute it uniformly along the projection of the loop onto the solar
surface. We assume a cross section corresponding to an EIS pixel, and thereby
obtain a spatially averaged EIS line profile for loop. Finally, a temporal
average is taken over the time required for the shock to travel to the other
end of the loop ($\approx 190$ s in this case). Such spatially and temporally
averaged line profiles from a single loop (e.g., Figure 6) is equivalent to an
observation of many unresolved loops of similar nature but at different stages
of their evolution (Patsourakos & Klimchuk, 2006; Klimchuk & Bradshaw, 2014).

Asymmetric coronal line profiles with blue wing enhancement are manifestations
of mass transport in the solar corona. Type II spicules are often suggested to
be associated with such a mass transport mechanism (De Pontieu et al., 2009,
2011; Martínez-Sykora et al., 2017). However, the extreme non-Gaussian shapes
of the simulated Fe XII and Fe XIV line profiles (Figure 6) are significantly
different from observed shapes (De Pontieu et al., 2009; Tian et al., 2011;
Tripathi & Klimchuk, 2013). Also, the very large blue shifts are inconsistent
with observations. Observed Doppler shifts of coronal lines tend to be slower
than $5$ km s-1 in both active regions (Doschek, 2012; Tripathi et al., 2012)
and quiet Sun (Chae et al., 1998; Peter & Judge, 1999). In contrast, a shift
of $150$ km s-1 is evident in the simulated spectral lines (Figure 6).

Our simulation is not reliable after the shock reaches the right boundary of
the model. Because of rigid wall boundary conditions, it reflects in an
unphysical manner. One might question whether the emission after this time
could dramatically alter the predicted line profiles. We estimate the
brightness of this neglected emission using the loop temperature profile
shortly before the shock reaches the chromosphere at $t=190$ s. The
temperature peaks at the shock, and there is strong cooling from thermal
conduction both to the left (up the loop leg) and, especially, to the right
(down the loop leg). We estimate the cooling timescale according to:

$\tau_{cond}=\frac{21}{2}\frac{k_{B}n_{e}l^{2}}{\kappa_{0\parallel}T^{5/2}}~{},$
(6)

where, $k_{B}$ is the Boltzmann constant, $\kappa_{0\parallel}$ is the
coefficient of thermal conductivity along the field lines, $T$ is the
temperature at the shock, $n_{e}$ is the electron number density behind the
shock, and $l$ is the temperature scale length. We do this separately using
the scale lengths on both sides of the shock, obtaining $\tau_{cond}=1290$ s
and 7 s for the left and right sides, respectively. Radiative cooling is much
weaker and can be safely ignored. We estimate the integrated emission after
$t=190$ s by multiplying the count rate at that time by the longer of the two
timescales, thereby obtaining an upper limit on the neglected emission in our
synthetic line profiles. The result is $10565$ DN pix-1 for Fe XII and $2206$
DN pix-1 for Fe XIV. These are about $0.97$ and $2.76$ times the temporally
integrated emission before this time, for Fe XII and Fe XIV, respectively. The
factors are much smaller using the shorter cooling timescale. Even the large
factors do not qualitatively alter our conclusions. The profile shapes and
Doppler shifts would still be much different from observed. The conclusions we
draw below are also not affected by neglecting the emission after the shock
reaches the right footpoint.

Figure 5: Loop integrated intensity in Fe XII and Fe XIV as a function of time
for Run1. The red dashed curve is the contribution only from the injected
spicule material, while the green dashed curve is the contribution of the loop
material excluding the injected spicule. The black curve is the sum of both
the spicule and non-spicule intensities. Time integrated intensities over the
190 s required for the shock to reach the right footpoint are indicated. Even
though the emission from the injected spicule material decreases rapidly, it
is so much brighter than the shock heated pre-existing material that it
contributes more to the time integration.

Figure 6: Spatially and temporally averaged synthetic spectral line profiles
of Fe XII and Fe XIV for Run1. The line profiles are averaged over the full
length of the loop and over the 190 s needed for the shock to reach the other
end of the loop. Spectral line profiles are highly blue-shifted, and the
shapes do not agree with the observed asymmetric line profiles.

###  3.3 Comparison with observations

We now estimate the spicule occurrence rate that would be required to explain
the observed coronal intensities from active regions and quiet Sun. We have
already seen that, in the absence of any external (coronal) heating, the hot
material at the tip of the spicule cools down rapidly. However, we are
concerned here with the total emission, including that from pre-existing
material that is heated as the spicule propagates along the loop. Consider a
region of area $\mathcal{A}_{reg}$ on the solar surface, large enough to
include many spicules. If the spatially averaged occurrence rate of spicules
in this region is $\mathcal{R}$ (cm-2 s-1), then one may expect
$\mathcal{N}_{reg}=\mathcal{R}\tau\mathcal{A}_{reg}$ spicules to be present at
any moment, where $\tau$ is the typical spicule lifetime. Since we are
averaging over large areas, the orientations of the spicule loops does not
matter, and we can treat the loops as straightened along the line of sight, as
done for Figure 5. If $\mathcal{I}_{sp}$ (DN s-1 pix-1) is the temporally
averaged intensity of such a loop (the full loop intensity divided by 190 s in
Fig. 5), then the expected intensity from a corona that only contains spicule
loops is
$\mathcal{I}_{obs}=\mathcal{N}_{reg}\mathcal{I}_{sp}\mathcal{A}_{sp}/\mathcal{A}_{reg}=\mathcal{I}_{sp}\mathcal{R}\tau\mathcal{A}_{sp}$,
where $\mathcal{A}_{sp}$ is the cross-sectional area of the loop.

The typical intensities ($\mathcal{I}_{obs}$) observed by EIS in active
regions and quiet Sun are, respectively, $162$ and $34$ DN s-1 pix-1 in Fe XII
(195 Å) and $35$ and $4$ DN s-1 pix-1 in Fe XIV (274 Å) (Brown et al., 2008).
On the other hand, the temporally averaged intensities from our simulation
($\mathcal{I}_{sp}$) are $56.36$ and $4.22$ DN s-1 pix-1 for Fe XII and Fe
XIV, respectively. Considering $\tau$ to be $190$ s, the time it takes for the
shock to travel across the loop, we derive an occurrence rate ($\mathcal{R}$)
of spicules as a function of their cross-sectional area ($\mathcal{A}_{sp}$).
Results are shown in Figure 7 for the two lines.

Following our earlier logic, we may also argue that at any given time there
are $\mathcal{N}_{\odot}=\mathcal{R}\tau\mathcal{A}_{\odot}$ spicules on the
solar disk, where $\mathcal{A}_{\odot}$ is the area of the solar disk. Using
the estimated value of the occurrence rate of the spicules ($\mathcal{R}$),
and taking $\tau$ to be $190$ s as before, the number of spicules on the solar
disk is related to the other quantities as per
$\mathcal{N}_{\odot}=(\mathcal{I}_{obs}/\mathcal{I}_{sp})(\mathcal{A}_{\odot}/\mathcal{A}_{sp})$.
This formula represents $\mathcal{N}_{\odot}$ as a function of the spicule
cross-sectional area $\mathcal{A}_{sp}$ (Figure 8). Considering the fact that
the typical observed widths of spicules lie between $200-400$ km (Pereira et
al., 2011), we find that the full disk equivalent number of spicules required
to explain the observed intensities exceeds $10^{7}$ in the quiet Sun and
$10^{8}$ in active regions, as indicated by the green shaded region in Figure
8. However, observational estimations for the number of spicules on the disk
vary between $10^{5}$ (Sterling & Moore, 2017) and $2\times 10^{7}$ (Judge &
Carlsson, 2010). There is a large discrepancy. Far more spicules than observed
would be required to produce all the observed coronal emission. For the quiet
Sun, $100$ times more spicules would be needed, while for active regions,
$10-10^{3}$ times more would be needed. These are lower limits based on Run1.
Our other simulations imply even greater numbers of spicules (see Table 2).

We should mention here that the larger the height the spicule rises, the
longer the time it compresses the ambient material, and thus the brighter the
time averaged emission. The spicules in our simulations with $150$ km s-1
injection speed reach a height of about $23$ Mm, which is much larger than the
typically observed spicule height ($\sim 10$ Mm). Therefore, we are likely to
overestimate the emission coming from spicule loops, and so the discrepancy
between the required and observed number of spicules is even greater. It
should also be noted that the values estimated by Sterling & Moore (2017);
Judge & Carlsson (2010) consider both type I & II spicules. The discrepancy
thus increases further if one considers type II spicules alone.

Figure 7: The required occurrence rate of spicules to explain the typically
observed Fe XII and Fe XIV intensities for active regions and quiet Sun
(curves are derived from Run1). The occurrence rate is plotted against the
cross-sectional area of spicules.

Figure 8: Required number of spicules at any given time on the Sun as a
function of spicule cross-sectional area, estimated from the spicule
occurrence rate (derived using Run1). Spicules having widths of $200-400$ km
are fairly common. The required number of spicules having widths in this range
is marked by the shaded region in both the plots.

Analysis of our simulated observations thus suggests that spicules contribute
a relatively minor amount to the emission and thermal energy of the corona.
Through the generation of shocks, they may heat the local plasma, but that too
cools down rapidly due to expansion and thermal conduction. Therefore,
synthetic spectra derived from our simulation show a high discrepancy with
observed spectra. However, this does not rule out the possibility of spicules
contributing significantly to the coronal mass. The ejected spicule material
may still get heated in the corona through some other heating mechanism – a
source that exceeds the initial thermal and kinetic energy of the spicule.
However, observational evidence of such a process is still lacking. Analyzing
the excess blue wing emission of multiple spectral lines hotter than $0.6$ MK,
Tripathi & Klimchuk (2013) have concluded that the upward mass flux is too
small to explain the mass of the active region corona. Their observations
indicate that spicules hotter than $0.6$ MK are not capable of providing
sufficient mass to the corona.

So far, we have allowed our spicules to propagate within a warm ($T=0.5$ MK),
relatively low density loop in order to determine whether they, by themselves,
can explain the observed hot emission. Our simulations indicate that this is
not viable. Therefore, there must be some other heating mechanisms at play
that produce the hot, dense plasma. Setting aside the issue of heating
mechanisms, in the following section we simply test the response of a spicule
in a hot and dense flux tube.

###  3.4 Spicule propagation in a hot loop

We have considered a static equilibrium loop with apex and footpoint
temperatures of approximately $2$ and $0.02$ MK, respectively. A spicule with
a tip temperature of $2$ MK followed by a cold, dense material with
temperature $0.02$ MK is injected with a velocity of $150$ km s-1 from the
bottom boundary, similar to our previous spicules. The velocity profile of the
injected spicule is the same as shown in Figure 2. The injected spicule
generates a shock that takes about $180$ s to traverse the loop.

The spatio-temporal averaged spectral line profiles are obtained following the
method described in Appendix C. However, in this case, because of the high
background temperature, the loop itself emits significantly in the Fe XII and
Fe XIV coronal lines. We consider the situation where the line of sight passes
through many loops. Some contain spicules and some are maintained in the hot
equilibrium state. We adjust the relative proportions to determine what
combination is able to reproduce the observed red-blue (RB) profile
asymmetries, which are generally $<0.05$ (Hara et al., 2008; De Pontieu et
al., 2009; Tian et al., 2011). For an asymmetry of $\approx 0.04$, we find
that the ratios of spicule to non-spicule strands are $1:150$ for the Fe XII
line and $1:72$ for the Fe XIV line. Again the conclusion is that spicules are
a relatively minor contributor to the corona overall, though they are
important for the loops in which they occur.

Figure 9: Spatially and temporally averaged synthetic spectral line profiles of Fe XII and Fe XIV for the case when a spicule (with tip temperature of $2$ MK) is injected in a hot loop with apex temperature $2$ MK. The line profiles are averaged over $180$ s (the time taken by the shock to traverse the loop). Table 2: Summarizing the number of spicules (width $\sim 300$ km) required to explain the quiet Sun and active region intensities as predicted from the test runs. Run | Ttip | v | Loop integrated counts | Quiet Sun | Active Region
---|---|---|---|---|---
| (MK) | (km s-1) | (DN s-1 pix-1) | (Required number of spicules) | (Required number of spicules)
|  |  | Fe XII | Fe XIV | Fe XII | Fe XIV | Fe XII | Fe XIV
1 | 2 | 150 | 0.6405 | $4.8\times 10^{-2}$ | $4.02\times 10^{7}$ | $6.31\times 10^{7}$ | $1.92\times 10^{8}$ | $5.52\times 10^{8}$
2 | 2 | 50 | 0.2336 | $1.47\times 10^{-2}$ | $1.1\times 10^{8}$ | $2.06\times 10^{8}$ | $5.25\times 10^{8}$ | $1.8\times 10^{9}$
3 | 1 | 150 | $6.7\times 10^{-2}$ | $1.3\times 10^{-3}$ | $3.86\times 10^{8}$ | $2.26\times 10^{9}$ | $1.84\times 10^{9}$ | $1.98\times 10^{10}$
4 | 1 | 50 | $5.9\times 10^{-3}$ | $7.81\times 10^{-6}$ | $4.3\times 10^{9}$ | $3.8\times 10^{11}$ | $2.06\times 10^{10}$ | $3.4\times 10^{12}$
5 | 0.02 | 150 | $4.4\times 10^{-4}$ | $1.5\times 10^{-7}$ | $5.89\times 10^{10}$ | $2.02\times 10^{13}$ | $2.8\times 10^{11}$ | $1.77\times 10^{14}$
6 | 0.02 | 50 | $1.02\times 10^{-7}$ | $1.1\times 10^{-13}$ | $2.5\times 10^{14}$ | $2.7\times 10^{19}$ | $1.2\times 10^{15}$ | $2.4\times 10^{20}$

##  4 Summary and discussion

The solar atmosphere displays a wide variety of spicules with different
temperatures and velocities. It has been suggested that type II spicules are a
major source of coronal mass and energy (De Pontieu et al., 2007, 2009, 2011).
In this work, we numerically investigate the role of spicules in producing
observed coronal emissions. In particular, we examine whether, in the absence
of any external heating, the hot tips of the spicules and the shock-heated
ambient plasma can explain the observed coronal emission. For this, we inject
spicules with different temperatures and velocities into a coronal loop in
static equilibrium. We choose a relatively cool equilibrium so that the loop
does not itself produce appreciable emission in the absence of a spicule. Each
of our injected spicules consists of a hot tip followed by a cold body. We
consider three different temperatures for the hot tips, viz., $2$, $1$ and
$0.02$ MK, while the cold, dense chromospheric plasma that follows the tip has
a temperature of $0.02$ MK. Six different simulations are run by injecting
each of these spicules with an initial velocity of either $50$ km s-1 or $150$
km s-1 (see Table 1). We also have constructed spectral line profiles and
estimated the spicule occurrence rate required to explain the observed
intensities from the quiet Sun and active regions. Our main results are
summarized as follows.

#### Shock formation during spicule propagation

All six runs described above suggest the formation of shocks due to the
injection of spicule material into the coronal flux tubes. The shocks are
stronger when the temperature differences and therefore pressure differences
with the ambient plasma are higher. Table 1 shows the variation of the
compression ratio (measure of shock strength) with changing temperature of the
spicule tip. The nature of the shock depends on the tip temperature. Spicules
with a hotter tip produce a pressure-driven shock that propagates with a speed
larger than the injection speed. Spicules with a cold tip (i.e.,
$T_{tip}=0.02$ MK) produce a piston-driven shock which propagates with a speed
close to the injection speed. The intensities and shapes of spectral line
profiles depend on the nature of the shock. The formation of shocks during
spicule injection agrees well with previous studies (Petralia et al., 2014;
Martínez-Sykora et al., 2018).

#### Rapid cooling of the hot spicule tip

Our simulations show that, in the absence of any external heating, the hot tip
of the spicule cools rapidly before reaching a substantial coronal height.
Consequently, the tip emission from coronal lines like Fe XII (195 Å) and Fe
XIV (274 Å) is short lived (Figure 5) and confined to low altitudes. The
result is consistent with earlier studies by Klimchuk (2012) and Klimchuk &
Bradshaw (2014).

#### Relative emission contributions of hot tip and shock heated plasma

Our simulations show that the pre-existing material in the loop gets heated
through shock compression and thermal conduction. However, the time-integrated
emission from this heated pre-existing material is less than that from the hot
tip, as shown in Figure 5. The tip plasma is hot for a much shorter time, but
it is inherently much brighter because of the greater densities (it is
injected in a dense state).

#### Line profile discrepancies

The shapes of our synthetic spectral line profiles show significant
discrepancies with observations. The simulated profiles are highly non-
Gaussian and far more asymmetric than observed. A strong blue shift ($\sim
150$ km s-1) of the synthetic lines is also inconsistent with the mild Doppler
shifts ($<5$ km s-1) observed in the quiet Sun and active regions.

#### Excessive number of spicules required to explain observed intensities

The spatially and temporally averaged intensities from our simulations
(Figures 6) imply that far more spicules are required to reproduce the
observed emission from the solar disk than are observed (Figure 8). The
discrepancies are up to a factor of $100$ for the quiet Sun and factors of
$10-10^{3}$ for active regions. These factors apply specifically to Run1,
where a spicule with a $2$ MK tip is ejected at a velocity of $150$ km s-1. As
listed in Table 2, the loops in our other simulations with different
combinations of tip temperature and ejection velocity are fainter, and
therefore more of them would be required to reproduce the observed disk
emission, exacerbating the discrepancy.

#### Ratio of loops with and without spicules

Under the assumption that the corona is comprised of hot loops with and
without spicule ejections, red-blue spectral line asymmetries similar to those
observed ($0.04$) require far more loops without spicules than with them. The
spicule to non-spicule loop number ratio is $1:150$ for the FeXII line and
$1:72$ for the Fe XIV line.

Our simulations indicate that spicules contribute a relatively minor amount to
the mass and energy of the corona. Such a claim had already been made by
Klimchuk (2012), where it was shown that hot tip material rapidly expanding
into the corona is unable to explain the observed coronal emission. However, a
bodily ejection of the spicule was not considered, and the emission from
ambient material effected by the expansion was not rigorously investigated
(though see Appendix B in that paper). Later, Petralia et al. (2014) argued
that the shock-heated material in front of an ejected cold spicule might be
erroneously interpreted as ejected hot material. They did not compare the
brightness of the shock-heated material with coronal observations. Our
numerical simulations improve on both of these studies. We show that neither
the expanding hot tip nor the shock-heated ambient material of a bodily
ejected spicule can reproduce coronal observations. A number of discrepancies
exist. The existence of some coronal heating mechanism - operating in the
corona itself - is required to explain the hot corona. It is not sufficient to
eject hot (or cold) material into the corona from below.

We emphasize that our conclusion does not rule out the possibility that waves
may be launched into the corona as part of the spicule ejection process, or
that new coronal currents may be created outside the flux tube in which the
ejected material resides, as suggested by Martínez-Sykora et al. (2018). Such
waves and currents would lead to coronal heating and could explain at least
some non-spicule loops. It seems doubtful, however, that this could explain
the many non-spicule loops implied by observed line profile asymmetries. It
seems that some type of heating unrelated to spicules must play the primary
role in explaining hot coronal plasma.

We thank the anonymous referee for her/his comments to improve the clarity of
the paper. SSM & AS thank Dr. Jishnu Bhattacharyya for many useful
discussions. Computations were carried out on the Physical Research
Laboratory’s VIKRAM cluster. JAK was supported by the Internal Scientist
Funding Model (competed work package program) at Goddard Space Flight Center.

##  Appendix A Static equilibrium configuration from double relaxation method

We inject spicules in a magnetic structure that is in static equilibrium. Such
an equilibrium is achieved recursively, and the final equilibrium profile is
obtained through two stages of relaxation. First, we obtain the density and
temperature profiles by solving the hydrostatic and energy balance equations
(Aschwanden & Schrijver, 2002) assuming a steady and uniform background
heating $Q_{bg}$. The CHIANTI radiative loss function $\Lambda(T)$ is used to
describe the loop’s radiation in the energy balance equation. The desired
looptop temperature is achieved by adjusting the value of $Q_{bg}$. However,
due to lack of exact energy balance, the temperature and density profiles
derived in this way do not achieve a perfect equilibrium state. Rather these
derived profiles are then used to calculate the final equilibrium loop profile
, such that the resulting temperature profile never drops below the
chromospheric temperature $T_{ch}$ ($2\times 10^{4}$ K), and the system does
not generate any spurious velocity either. In the following, we explain these
two stages in detail.

###  A.1 Heating and cooling in Relaxation-I:

Starting with the initial profiles described above, the loop is allowed to
relax under gravity with the constant background heating $Q_{bg}$. To avoid
numerical artifacts, from this stage onward, we smoothly reduce the radiative
cooling of the chromosphere to zero over a narrow temperature range between
$T_{ch}$ and $T_{min}$, where $T_{min}=1.95\times 10^{4}$ K is a conveniently
chosen temperature slightly less than $T_{ch}$. This is achieved by the
radiative loss function $\lambda(T)$, defined as

$\lambda(T)=\begin{cases}\Lambda(T),&\text{if}\,T\geq T_{ch}\\\
\left(\frac{T-T_{min}}{T_{ch}-T_{min}}\right)\Lambda(T_{ch}),&\text{if}\,T_{min}<T<T_{ch}\\\
0,&\text{if}\,T\leq T_{min}\\\ \end{cases}~{}.$ (A1)

Here $\Lambda(T)$ denotes the optically thin radiative loss function from
CHIANTI. The modified function $\lambda(T)$ is plotted in Figure 10. As the
loop relaxes, material drains from the corona and accumulates at the
footpoints. The resulting high density of the loop footpoints gives rise to
excessive cooling and brings down the footpoint temperatures below $T_{min}$,
along with generating short lived velocities. However, the loop eventually
achieves a steady-state, and we use the enhanced footpoint density at that
time ($n_{base}$) to estimate the additional heating required to keep the
chromospheric temperature above $T_{min}$. This is carried out in the next
relaxation stage.

###  A.2 Heating and cooling in Relaxation-II:

Once again, we start with the initial density and temperature profiles from
the beginning of the first stage. However, this time we apply additional
heating in the chromosphere above the constant background heating $Q_{bg}$.
This prevents the plasma from cooling below $T_{min}$ and instead lets it
hover between $T_{ch}$ and $T_{min}$. The total heating function $Q$ is given
by

$Q=\begin{cases}Q_{bg},&\text{if}\,T\geq T_{ch}\\\
\left(\frac{n}{n_{base}}\right)^{2}Q_{ch}\left(\frac{T_{ch}-T}{T_{ch}-T_{min}}\right)+Q_{bg},&\text{if}\,T_{min}<T<T_{ch}\\\
\left(\frac{n}{n_{base}}\right)^{2}Q_{ch}+Q_{bg},&\text{if}\,T\leq
T_{min}\end{cases}~{},$ (A2)

where $Q_{ch}=n_{ch}^{2}\Lambda(T_{ch})$ is the heat required to balance the
radiative losses from the footpoint plasma of the initial loop profile at
temperature $T_{ch}$ and density $n_{ch}$. Figure 10 graphically depicts the
radiative loss and heating functions that are maintained throughout the
simulation.

Figure 10: Radiative loss function and excess heating implemented to prevent
cooling of plasma below $T_{min}$. When the temperature of the loop is above
$T_{ch}$, the radiative loss is given by the CHIANTI radiative loss function
$\Lambda(T)$ (Equation A1). The only heating applied to the loop at these
temperatures is the uniform background heating $Q_{bg}$. The radiative cooling
smoothly goes to zero as the temperature approaches $T_{min}$ from above. In
the temperature range $T_{min}$ to $T_{ch}$, depending on the plasma density
and temperature, an additional heating is provided to the loop to balance the
lost energy through radiative cooling. Below the temperature $T_{min}$, an
additional heating $Q_{ch}$, proportional to $n_{ch}^{2}$, is applied to bring
the plasma back $T_{ch}$.

##  Appendix B Variation of shock speed with height

For a pressure driven shock, the shock’s speed primarily depends on the
pressure difference between the spicule’s tip and the ambient medium in which
it is propagating. Lower pressure close to the loop apex provides lesser
resistance to the shock propagation, and hence the shock speed increases. On
the other hand, high pressure close to the footpoints provides greater
resistance and thus the shock speed reduces. For a better understanding, we
track the shock front along the loop and derive its speed during its
propagation. The shock front at any instant can be identified from the density
jump moving ahead of the injected spicule material. To track it, we identify
the jump in density at each time, which is also associated with the maximum
temperature of the loop. Once the locations of the shock front along the loop
are identified, a derivative of the same gives the instantaneous shock speed
as a function of loop coordinates. Figure 11 shows the variation of shock
speed as a function of loop coordinates for three different shocks, all
ejected with velocity $150$ km s-1 but with three different tip temperatures,
viz. $2$, $1$ and $0.02$ MK. Though the shock speeds increase at the loop apex
for all three shocks, velocity amplitudes depend on the injection temperatures
and thus pressures. The larger the tip temperature, the higher the spicule tip
pressure and hence larger is the shock speed.

Figure 11: Variation of the shock speeds as a function of loop coordinates for
injection speed = $150$ km s-1 and tip temperatures = $2$, $1$ and $0.02$ MK.

##  Appendix C Forward Modelling

Spectral profiles provide a wealth of information about the plasma dynamics
along the line of sight (LOS). Adapting the method outlined in Patsourakos &
Klimchuk (2006), synthetic spectral line profiles are constructed at each
numerical grid cell using the cell’s density, velocity and temperature. At any
given time, $t$, and location along the loop, $s$, the line profile is

$I(s,t)=\frac{I_{0}}{\sqrt{\pi}v_{\text{width}}}\exp\left[\frac{-(v-v_{\text{shift}})^{2}}{v_{\text{width}}^{2}}\right]~{},$
(C1)

where $I_{0}$ is the amplitude, $v_{\text{shift}}$ is the Doppler shift, and
$v_{\text{width}}$ is the thermal line width. The amplitude is given by

$I_{0}(s,t)=n_{e}^{2}G(T)ds~{},$ (C2)

where $n_{e}$, $T$ and $ds$ denote the electron number density, temperature,
and length of the cell. The contribution function $G(T)$ for the line is taken
from the CHIANTI atomic data base (Landi et al., 2013). The Doppler shift
equals the line of sight velocity of the cell,

$v_{\text{shift}}=v_{\text{los}}~{}$ (C3)

in wavelength units, and the thermal width is given by

$v_{\text{width}}=\sqrt{\frac{2k_{B}T}{m_{ion}}}~{},$ (C4)

where $m_{\text{ion}}$ is the mass of the ion.

Once the line profile at each grid point is constructed, spatial averaging is
performed by summing the profiles along the loop and dividing by its projected
length assuming that it lies in a vertical plane and is viewed from above:

$\langle I(t)\rangle_{\text{spatial}}=\frac{\pi}{2L}\sum_{s}I(s,t)\times d$
(C5)

where $L$ is the loop length and $d$ is the pixel dimension. The loop is
assumed to have a cross section of $d^{2}$. Finally the spatially averaged
line profiles are temporally averaged over a time $\tau$, which is taken to be
the travel time of the shock along the loop; this yields

$\langle I\rangle_{\text{spatial, temporal}}=\frac{1}{\tau}\sum_{t}\langle
I(t)\rangle_{\text{spatial}}$ (C6)

## References

  * Aschwanden & Schrijver (2002) Aschwanden, M. J., & Schrijver, C. J. 2002, ApJS, 142, 269
  * Athay & Holzer (1982) Athay, R. G., & Holzer, T. E. 1982, ApJ, 255, 743
  * Bradshaw & Klimchuk (2015) Bradshaw, S. J., & Klimchuk, J. A. 2015, ApJ, 811, 129
  * Brown et al. (2008) Brown, C. M., Feldman, U., Seely, J. F., Korendyke, C. M., & Hara, H. 2008, ApJS, 176, 511
  * Cargill et al. (2022) Cargill, P. J., Bradshaw, S. J., Klimchuk, J. A., & Barnes, W. T. 2022, MNRAS, 509, 4420
  * Chae et al. (1998) Chae, J., Yun, H. S., & Poland, A. I. 1998, ApJS, 114, 151
  * Culhane et al. (2007) Culhane, J. L., Harra, L. K., James, A. M., et al. 2007, Solar Physics, 243, 19
  * De Pontieu et al. (2009) De Pontieu, B., McIntosh, S. W., Hansteen, V. H., & Schrijver, C. J. 2009, ApJ, 701, L1
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S., Hansteen, V. H., et al. 2007, PASJ, 59, S655
  * De Pontieu et al. (2011) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2011, Science, 331, 55
  * Dey et al. (2022) Dey, S., Chatterjee, P., Murthy, O. V. S. N., et al. 2022, Nature Physics, 18, 595
  * Doschek (2012) Doschek, G. A. 2012, ApJ, 754, 153
  * Guarrasi et al. (2014) Guarrasi, M., Reale, F., Orlando, S., Mignone, A., & Klimchuk, J. A. 2014, A&A, 564, A48
  * Hara et al. (2008) Hara, H., Watanabe, T., Harra, L. K., et al. 2008, ApJ, 678, L67
  * Judge & Carlsson (2010) Judge, P. G., & Carlsson, M. 2010, ApJ, 719, 469
  * Klimchuk (2006) Klimchuk, J. A. 2006, Sol. Phys., 234, 41
  * Klimchuk (2012) —. 2012, Journal of Geophysical Research (Space Physics), 117, A12102
  * Klimchuk (2015) —. 2015, Philosophical Transactions of the Royal Society of London Series A, 373, 20140256
  * Klimchuk & Bradshaw (2014) Klimchuk, J. A., & Bradshaw, S. J. 2014, ApJ, 791, 60
  * Kosugi et al. (2007) Kosugi, T., Matsuzaki, K., Sakao, T., et al. 2007, Solar Physics, 243, 3
  * Landi et al. (2013) Landi, E., Young, P. R., Dere, K. P., Del Zanna, G., & Mason, H. E. 2013, ApJ, 763, 86
  * Martínez-Sykora et al. (2017) Martínez-Sykora, J., De Pontieu, B., Hansteen, V. H., et al. 2017, Science, 356, 1269
  * Martínez-Sykora et al. (2018) Martínez-Sykora, J., Pontieu, B. D., Moortel, I. D., Hansteen, V. H., & Carlsson, M. 2018, The Astrophysical Journal, 860, 116
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mikić et al. (2013) Mikić, Z., Lionello, R., Mok, Y., Linker, J. A., & Winebarger, A. R. 2013, ApJ, 773, 94
  * Patsourakos & Klimchuk (2006) Patsourakos, S., & Klimchuk, J. A. 2006, ApJ, 647, 1452
  * Patsourakos et al. (2014) Patsourakos, S., Klimchuk, J. A., & Young, P. R. 2014, ApJ, 781, 58
  * Pereira et al. (2011) Pereira, T. M., De Pontieu, B., & Carlsson, M. 2011, in AGU Fall Meeting Abstracts, Vol. 2011, SH34B–01
  * Peter & Judge (1999) Peter, H., & Judge, P. G. 1999, ApJ, 522, 1148
  * Petralia et al. (2014) Petralia, A., Reale, F., Orlando, S., & Klimchuk, J. A. 2014, A&A, 567, A70
  * Pneuman & Kopp (1977) Pneuman, G. W., & Kopp, R. A. 1977, A&A, 55, 305
  * Pneuman & Kopp (1978) —. 1978, Sol. Phys., 57, 49
  * Roberts (1945) Roberts, W. O. 1945, ApJ, 101, 136
  * Samanta et al. (2019) Samanta, T., Tian, H., Yurchyshyn, V., et al. 2019, Science, 366, 890
  * Secchi (1877) Secchi, P. 1877, Le soleil. 2e partie (Gauthier-Villars)
  * Skogsrud et al. (2015) Skogsrud, H., Rouppe van der Voort, L., De Pontieu, B., & Pereira, T. M. D. 2015, ApJ, 806, 170
  * Sterling & Moore (2017) Sterling, A. C., & Moore, R. L. 2017, in AGU Fall Meeting Abstracts, Vol. 2017, SH43A–2791
  * Tian et al. (2011) Tian, H., McIntosh, S. W., Pontieu, B. D., et al. 2011, The Astrophysical Journal, 738, 18
  * Tripathi & Klimchuk (2013) Tripathi, D., & Klimchuk, J. A. 2013, ApJ, 779, 1
  * Tripathi et al. (2012) Tripathi, D., Mason, H. E., & Klimchuk, J. A. 2012, ApJ, 753, 37
  * Viall et al. (2021) Viall, N. M., De Moortel, I., Downs, C., et al. 2021, in Solar Physics and Solar Wind, ed. N. E. Raouafi & A. Vourlidas, Vol. 1, 35
  * Withbroe (1983) Withbroe, G. L. 1983, ApJ, 267, 825