††thanks: ∗E-mail<EMAIL_ADDRESS>
,

# Detached Circumstellar Matter as an Explanation for Slowly-Rising
Interacting Type Ibc Supernovae

Yuki Takei1,2,3,∗ Daichi Tsuna4,2 1Yukawa Institute for Theoretical Physics,
Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto, Kyoto 606-8502,
Japan  2Research Center for the Early Universe, Graduate School of Science,
University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan  3Astrophysical Big Bang
Laboratory, RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan  4TAPIR,
Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA

###### Abstract

Some hydrogen-poor (Type Ibc) supernovae (SNe) are known to have massive
circumstellar matter (CSM) that are well detached from the star. Using the
open-source code CHIPS, we construct a grid of models of SN Ibc interacting
with detached CSM. We find that interaction with detached CSM can produce a
slowly rising phase in the light curve seen in some interacting SN Ibc, which
is difficult to reproduce by interaction with CSM of a density profile
motivated from wind or eruptions that are continuous down to the star. We also
show that SNe having double peaks in their light curves with timescales of
months (e.g., SN 2022xxf) can be explained by radioactive decay of 56Ni/56Co,
followed by interaction with detached CSM.

###### keywords:

circumstellar matter — supernovae: general — stars: mass-loss

##  1 Introduction

Some massive stars are known to undergo episodic mass loss just before core-
collapse to form a dense circumstellar matter (CSM) (e.g., Chugai et al.,
2004; Pastorello et al., 2007; Mauerhan et al., 2013; Moriya et al., 2014;
Bruch et al., 2021; Strotjohann et al., 2021). Supernovae (SNe) with dense CSM
show bright light curves and prominent narrow emission lines, powered by
collision of the SN ejecta with the CSM. Studying the CSM is essential for
understanding massive star evolution, as it provides direct insights into the
dramatic pathways leading to its death.

For hydrogen-poor stripped envelope SNe, the structure of the dense CSM is
found to be diverse. The representatives of such interacting stripped-envelope
SNe has been Type Ibn (e.g., Pastorello et al., 2008), and the more recently
discovered Type Icn (e.g., Gal-Yam et al., 2022). The signs of CSM interaction
from the very early phase of the SN, and the typically fast light curve rise
of a few days, indicate confined CSM close to the star ($\lesssim 10^{15}$
cm).

More recently a class of double-peaked SNe with a luminous second peak in the
optical has been identified (Kuncarayakti et al., 2022, 2023; Das et al.,
2023). In the broad-lined Type Ic SN 2022xxf (Kuncarayakti et al., 2023)
representative of this class, the second peak slowly emerged a month after the
first peak presumably powered by the radioactive decay of 56Ni in the ejecta.
Some of these SNe show signatures of CSM interaction around the second peak,
implying a massive detached CSM far from the progenitor ($\gtrsim
10^{15}$–$10^{16}$ cm) in contrast to typical Type Ibn/Icn with confined dense
CSM. There are also some SNe Ib which show re-brightening likely due to
interaction with detached hydrogen-rich CSM, such as SN 2014C and SN 2019yvr
(Milisavljevic et al., 2015; Margutti et al., 2017; Ferrari et al., 2024).

While these recent events have presented strong evidence for the existence of
detached CSM in some SNe, detailed simulations on the diversity of the light
curves expected from these CSM are yet to be conducted. In this work we
simulate a grid of bolometric light curves from detached CSM of various radii,
using an open-source code CHIPS (Takei et al., 2022, 2024). Using a model for
the density profile of the detached CSM inspired from the model of Tsuna &
Takei (2023), we show that detached CSM can lead to interesting SN light curve
features difficult to reproduce by typically assumed density profiles of
continuous CSM down to the star, such as slow rises and long-duration double-
peaked structures.

This work is constructed as follows. In Section 2, we present our model for
constructing the detached CSM and calculating light curves powered by CSM
interaction, and investigate their basic properties. In Section 3 we compare
our model to observations of interacting SN Ibc with slowly-rising light curve
features, namely a Type Ibn SN 2010al and the aforementioned SN 2022xxf.
Section 4 provides discussion, including some of the limitations of our model
that can be improved in future works.

##  2 Our Model

Here we present our framework for constructing the detached CSM and modeling
the bolometric light curves. We first summarize our light curve modeling code
CHIPS for interacting SN Ibc (Takei et al., 2024), and the main
characteristics of the light curves obtained in that work. The model by Takei
et al. (2024) considered CSM that are continuous down to the star, and their
results offer a motivation for this work to consider detached CSM. Then we
describe the modeling for the detached CSM newly considered in this work, and
conduct a parameter study of the CSM to inspect the diversity in the light
curves.

###  2.1 Light Curve Properties for Continuous CSM

Takei et al. (2024) conducted a parameter study of interacting Type Ibc SNe
using the open-source code CHIPS, aiming to reproduce the light curves of Type
Ibn/Icn SNe. The light curve peak is powered by shocks formed via interaction
of the SN ejecta and the dense CSM, where the latter is assumed to be
generated from mass eruption from the star. CHIPS is capable of constructing
CSM structures by one-dimensional radiation hydrodynamical simulations of mass
eruption due to energy injection in the envelope (Kuriyama & Shigeyama, 2020),
and simulating the SN light curves due to collision of this CSM and the SN
ejecta.

From their simulations of mass eruption, the density profile of the CSM at
core-collapse is found to generally follow a double power-law continuous down
to the star (solid line in Figure 1; see also Tsuna et al., 2021b; Ko et al.,
2022),

$\displaystyle\rho_{\mathrm{CSM}}(r)=\rho_{*}\left[\frac{(r/r_{*})^{1.5/y_{*}}+(r/r_{*})^{n_{\mathrm{out}}/y_{*}}}{2}\right]^{-y_{*}},$
(1)

where $r_{*},\,\rho_{*}$ denote the radius and density at the transition
point, $n_{\mathrm{out}}\approx 9$–$10$ is the power-law index of the outer
CSM, and $y_{*}\approx 2$–$4$ determines the curvature at the transition
point. Depending on the injection energy, the mass of the (unbound) CSM was
found to have a range $\approx 0.01$–$0.2~{}M_{\odot}$. Takei et al. (2024)
then conducted light curve modeling of the subsequent SNe by radiative
transfer calculations in the CHIPS code, with the explosion energy $E_{\rm
ej}$ as an additional parameter. The code assumes local thermal equilibrium
(LTE) and steady state in the shocked region, to numerically model the shock
propagation for an arbitrarily CSM density profile with low computational cost
(for details see Takei & Shigeyama, 2020). They found that a typical $\sim
10^{51}$ erg explosion in a CSM of $\sim 0.01$–$0.1M_{\odot}$, generated from
the above mass eruption at $\lesssim 1$ yr before core-collapse, is sufficient
to power the bright light curve peak observed in typical Type Ibn/Icn SNe (see
their Figure 8).

An important shortcoming of their work was that the rise in the light curves
was rapid with timescales of mostly a few days or less. This timescale is
consistent with SN Ibn/Icn found in high-cadence surveys (e.g., Ho et al.,
2023) but fails to reproduce a subset of Type Ibn with long rise times of
$\gtrsim 10$ days (Pastorello et al., 2015a, b, 2016). This is mainly due to
the rise time in the light curves being set by the diffusion time in the dense
CSM. The required CSM mass to reproduce a rise time $t_{\rm rise}$ by
diffusion is (Moriya, 2015; Tsuna et al., 2024)

$\displaystyle M_{\mathrm{CSM}}$ $\displaystyle\approx$
$\displaystyle\frac{4\pi ct_{\rm rise}^{2}v_{\rm sh}}{\kappa}$ (2)
$\displaystyle\sim$ $\displaystyle
1.4\,M_{\odot}\left(\frac{t_{\mathrm{rise}}}{10\,{\rm
day}}\right)^{2}\left(\frac{v_{\mathrm{sh}}}{10^{4}\,{\rm
km\,s^{-1}}}\right)\left(\frac{\kappa}{0.1\,{\rm
cm^{2}\,g^{-1}}}\right)^{-1},$

where $c$ is the speed of light, $v_{\rm sh}$ is the velocity of the shock
sweeping the CSM, and $\kappa\approx 0.1\ {\rm cm^{2}\ g^{-1}}$ is the opacity
in the unshocked CSM expected for hydrogen-poor gas (e.g., Kleiser & Kasen,
2014; Piro & Morozova, 2014). This is an order of magnitude larger than the
largest CSM mass obtained in Takei et al. (2024). This mass may be achieved by
a near-complete envelope ejection with even larger injection energies than
considered in Takei et al. (2024), or possibly by a much more massive star
with initial mass of $\sim 100\,M_{\odot}$ that undergoes e.g. pulsational
pair-instability (Woosley et al., 2002; Woosley, 2019). However, the huge
dissipation of the SN ejecta’s kinetic energy by such massive CSM is expected
to make the light curve much brighter than other SN Ibn. These slow-rising SN
Ibn have peak luminosities similar to other SN Ibn samples, which appears
inconsistent with this expectation.

Recently Moriya (2023) raised an alternative possibility for the long rise,
arguing that some Type IIn SNe with very long rises can be realized for CSM
with a flat or rising density structure ($s>-1.5$ where $\rho_{\rm CSM}\propto
r^{s}$). A density rising with radius is naturally expected for detached CSM,
where the faraway dense CSM surrounds a low-density bubble111A model for
detached CSM (Tsuna & Takei, 2023) predicts that the dense CSM has to expand
to large radii before it can be detached from the star, which may also be
favorable for producing long rise times.. In the following we conduct light
curve modeling using a simplified density profile expected from a detached
CSM, and show that this indeed helps to produce light curves having long rises
of several tens of days, as well as multiple months-long peaks if significant
amounts of 56Ni are present in the SN ejecta.

###  2.2 Constructing the Detached CSM

Figure 1: Example density profiles of detached CSM ($x_{\rm comp}=3$) with
$R_{\rm bubble}=5\times 10^{15}$, $10^{16}$ cm (see Section 2.2 and equation 6
for the parameters). The original CSM profile without detachment (equation 1),
created by eruption of a $25~{}M_{\odot}$ helium star model, is plotted as a
solid line.

The mechanisms for forming a detached CSM still remain uncertain. Such
detached CSM may form due to binary interaction, as is discussed in the
context of SN 2014C (Sun et al., 2020; Brethauer et al., 2022). Recently Tsuna
& Takei (2023) suggested a model where the CSM is generated by partial mass
eruption(s) from the progenitor (see their Figure 1). In their model the
relative strengths between the fallback of the bound CSM and feedback from the
leftover progenitor (such as ram pressure from the stellar wind) can produce
both confined CSM of $\lesssim 10^{15}$ cm seen in SNe Ibn/Icn, and detached
CSM surrounding a low-density wind bubble extending to $R_{\rm
bubble}\gtrsim(10^{15}$–$10^{16})$ cm.

For both of these models for detached CSM, the circumstellar environment
consists of an inner low-density wind from the leftover star and an outer
dense CSM. Between these two we expect a (likely shocked) interface, generated
by collision of the faster wind that catches up with the slower CSM shell.
Such structures can be regarded as a scaled-down version of pulsar wind
nebulae, where the fast pulsar wind sweeps the slower and denser SN ejecta,
and creates a shocked interface in between that powers multi-wavelength
emission (e.g., Kennel & Coroniti, 1984; Chevalier & Fransson, 1992; Suzuki &
Maeda, 2021).

In addition to uncertainties in the mechanism for generating detached CSM, its
location and structure are likely complicated by various physical processes in
the shell. First, radiative cooling can be important in the shocked CSM222Note
that we consider the pre-explosion CSM shocked by the stellar wind, and not by
the SN ejecta. owing to its high density, which would compress the shell to
much higher densities than the unshocked part of the CSM. Second, the
interface between the shocked wind and shocked CSM is likely subject to
Rayleigh-Taylor instabilities, which would mix the interface and produce
asymmetric and clumpy density structures in the shocked region. Rayleigh-
Taylor mixing may make the density profile in the shell less sharp, as seen in
multi-dimensional simulations for pulsar-wind nebulae (e.g., Jun, 1998, their
Figure 7).

In this work, we adopt a simple one-dimensional model for the density profile
of the detached CSM with minimal model parameters: the outer radii of the
shell-bubble interface $R_{\mathrm{bubble}}$, and the compression ratio
$x_{\mathrm{comp}}(>1)$ defined below. We also assume for simplicity that the
(angle-averaged) density profile of the shell obeys a power-law in radius
$\rho\propto r^{s}$.

Mathematically, the CSM we assume here has a density profile given by

$\displaystyle\rho_{\mathrm{CSM,\,detach}}(r)=$
$\displaystyle\left\\{\begin{array}[]{ll}\dot{M}/4\pi
v_{\mathrm{w}}r^{2}&({\rm unshocked\ wind};\
r\leq(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}),\\\ qr^{s}&({\rm shell};\
(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}<r\leq R_{\mathrm{bubble}}),\\\
\rho_{\mathrm{CSM}}(r)&({\rm unshocked\ CSM};\
R_{\mathrm{bubble}}<r),\end{array}\right.$ (6)

where $\dot{M},\,v_{\mathrm{w}}$ denote the mass-loss rate and velocity of the
wind from the leftover star. We adopt the density profile of equation (1) for
the unshocked CSM at $r>R_{\rm bubble}$, and update the profile within $R_{\rm
bubble}$ to construct the detached shell.

We obtain the profile of the unshocked CSM $\rho_{\rm CSM}(r)$ using the mass
eruption calculations of the CHIPS code. To simulate mass eruption, CHIPS
requires an initial stellar model made by MESA (Paxton et al., 2011, 2013,
2015, 2018, 2019; Jermyn et al., 2023). We prepare a helium star model as a
progenitor of SN Ib, and a helium-poor star for SN Ic.

In the model of Tsuna & Takei (2023), the detachment of the star is considered
to be driven by the ram pressure of the stellar wind, which tends to be
stronger for more massive stars with larger luminosity (Eddington ratio; e.g.,
Sander & Vink, 2020). For the SN Ib progenitor, we adopt a helium star model
of ZAMS mass of $25\,M_{\odot}$, which is one of the most luminous progenitors
created in Takei et al. (2024) using version 12778 of MESA. For a SNe Ic
progenitor, we newly create a helium-deficient star with zero-age main
sequence (ZAMS) mass of $35\,M_{\odot}$ that is more massive than the models
previously generated in Takei et al. (2024). As was done in the
$25\,M_{\odot}$ model we adopt the test suite example_make_pre_ccsn, with the
standard “Dutch” wind mass-loss scheme for both hot and cold winds and mixing-
length parameter of $\alpha_{\rm MLT}=3$ ($1.5$) in hydrogen-rich(poor)
regions. For the SN Ic model, the hydrogen and helium layers are quickly
removed after core hydrogen and helium depletion respectively, by the option
`mass_change`. The resultant mass and radius of this star just before core-
collapse are $\approx 7.3\,M_{\odot}$ and $0.3\,R_{\odot}$
respectively333While we use the star at core-collapse as initial condition for
the eruption, we have checked that the stellar radius changes only by
$\lesssim 10\%$ in the final decade of its life, and mass by a much smaller
amount. Therefore, to first order the hydrodynamics of the mass eruption is
expected to be the insensitive to when the eruption is assumed to happen (i.e.
choice of $t_{\rm inj}$)..

Figure 2: Comparison of light curves changing two key parameters:
$R_{\mathrm{bubble}},\,E_{\mathrm{ej}}$ using $25\,M_{\odot}$ He-star model.
$M_{\mathrm{Ni}}$ is fixed to $0\,(0.4)\,M_{\odot}$ in the top (bottom)
panels. The fixed parameters are shown on the top of each panel. The
dependence of light curves on $R_{\mathrm{bubble}},\,E_{\mathrm{ej}}$ are
shown in left, and right panels, respectively. For comparison, the light
curves of an SN interacting with the continuous CSM are indicated with red
dashed lines, using the same parameters of $E_{\mathrm{ej}}=3\times
10^{51}\,{\rm erg},\,t_{\mathrm{inj}}=1.2\,{\rm yr}$.  Figure 3: Same as
Figure 2, but on the dependence on $t_{\mathrm{inj}}$ for two values of
$R_{\rm bubble}$, $4\times 10^{15}$ cm and $10^{16}$ cm.

We then simulate mass eruption by injecting energy at the base of the outer
layer of the star, as done in Takei et al. (2024). The input energy is given
in the code as a dimensionless parameter $f_{\mathrm{inj}}$, being scaled with
the binding energy of this outer layer ($-7.5\times 10^{49}\,{\rm erg}$ for
the $25\,M_{\odot}$ model and $\approx-8.5\times 10^{50}\,{\rm erg}$ for
$35\,M_{\odot}$ model). In this work we choose energies of
$f_{\mathrm{inj}}=0.7$ ($35\,M_{\odot}$ model),$\,0.8$ ($25\,M_{\odot}$ model)
to be injected over a period of $10^{-2}\,{\rm s}$, and follow the evolution
of the erupted envelope until $t=t_{\mathrm{inj}}$ after energy injection. The
resulting masses of the erupted CSM for the two models are both $\approx
0.2\,M_{\odot}$. To simplify the parameter space, we vary $t_{\rm inj}$ (or
CSM extent) but keep $f_{\rm inj}$ (or total CSM mass) fixed to the above
values. For fixed $t_{\rm inj}$, the effect of reducing $f_{\rm inj}$ is to
reduce the CSM mass, which generally make the light curve from interaction
dimmer and evolve slightly faster (see also Figure 8 of Takei et al., 2024).

Once we simulate $\rho_{\rm CSM}(r)$ in the absence of detachment, we convert
this to a detached CSM profile $\rho_{\rm CSM,detach}$ in equation (6), with
$R_{\rm bubble}$ and $x_{\rm comp}$ as additional input parameters. The
constants $q,\,s$ that characterize the detached shell are determined by two
equations, the density at the shocked wind and mass conservation of the CSM.
For the former, we assume that the density at the shocked wind at
$r=(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}$, which is subject to much
weaker radiative cooling, is simply increased from the unshocked wind with
compression ratio $x_{\rm comp}$,

$\displaystyle
q\left[(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}\right]^{s}=\frac{x_{\rm
comp}(\dot{M}/4\pi
v_{w})}{[(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}]^{2}}.$ (7)

For the latter, the total mass of the CSM is assumed to be the same as the
case if the profile had been continuous down to the star. In other words, the
shell has swept up mass occupied by the dense CSM had the wind been absent,
i.e.

$\displaystyle\int_{(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}}^{R_{\mathrm{bubble}}}qr^{s}4\pi
r^{2}dr$
$\displaystyle=\int_{0}^{(1-x_{\mathrm{comp}}^{-1})R_{\mathrm{bubble}}}\rho_{\rm
CSM}(r)4\pi r^{2}dr,$ (8)

where $\rho_{\rm CSM}(r)$ corresponds to that in equation (1).

Example density profiles of our detached CSM model are shown in Figure 1,
together with the profile for continuous CSM as comparison. Because the
unshocked CSM is much denser than the unshocked wind, we find $s$ in the shell
to be positive, i.e. a density increasing with radius. As we show below, this
can lead to slowly rising light curves as originally suggested in Moriya
(2023).

While the velocity of the unshocked CSM (bottom of equation 6) is obtained
from our simulations, the velocity in the unshocked wind and shell are more
uncertain. For simplicity, these velocities are fixed to $3000$ km s-1
independent of radius. While the light curves would depend little on the
velocity of the (low-density) unshocked wind, it would depend more on the
velocity for the shell $v_{\rm shell}$, as the ejecta’s kinetic energy
dissipated by the interaction with the shell roughly scales as $(v_{\rm
sh}-v_{\rm shell})^{2}$ (e.g., Murase et al., 2014). A lower $v_{\rm shell}$
with enhanced kinetic energy dissipation would make the light curve brighter
and last slightly longer, as we see later in Figure 4.

###  2.3 Light Curves and Their Parameter Dependences

In this work we explore the dependence of the light curves on physical
parameters related to the explosion ($E_{\mathrm{ej}},\,M_{\mathrm{Ni}}$) and
the detached CSM ($t_{\mathrm{inj}},\,R_{\mathrm{bubble}}$). Here $M_{\rm Ni}$
is the mass of radioactive 56Ni synthesized in the explosion, and we include
the time-dependent emission from radioactive heating of the ejecta following
Takei et al. (2024). The other parameters are here fixed for simplicity to
$x_{\mathrm{comp}}=3$, and $\dot{M}/v_{\mathrm{w}}=10^{-7}M_{\odot}\,{\rm
yr^{-1}}/({\rm km\>s^{-1}})$, unless otherwise mentioned. The detailed choice
of these two parameters do not qualitatively affect the conclusions of this
section.

In Figure 2 and 3 we plot the bolometric light curves varying $10^{15}\,{\rm
cm}\leq R_{\mathrm{bubble}}\leq 2\times 10^{16}\,{\rm cm},\,10^{51}\,{\rm
erg}\leq E_{\mathrm{ej}}\leq 10^{52}\,{\rm erg},\,0.6\,{\rm yr}\leq
t_{\mathrm{inj}}\leq 2\,{\rm yr},$ and
$M_{\mathrm{Ni}}=\\{0,\,0.4\\}\,M_{\odot}$. The values for $E_{\rm ej}$ and
$M_{\rm Ni}$ cover the two extremes, from Type Ibn/Icn consistent with small
or no 56Ni (e.g., Moriya & Maeda, 2016; Pellegrino et al., 2022) to broad-
lined Type Ic with high $E_{\rm ej}$ and $M_{\rm Ni}$ (Taddia et al., 2019).

Generally, CSM interaction is found to power emission up to
$10^{42}$–$10^{43}$ erg s-1, which is comparable to or higher than that from
56Ni decay. The light curves from detached CSM emerge later and slower
compared to the case of continuous CSM (i.e. the density profile simply
follows equation 1 to the center), shown as red dashed lines. The radius of
the detached CSM $R_{\rm bubble}$ mainly governs when this peak appears, and
can create a diversity in the light curve shapes as discussed below.

The left panels of Figure 2 show the dependence of light curves on the
location of the CSM shell $R_{\mathrm{bubble}}$. For the top panel without
nickel decay, the rise time is longer for larger $R_{\mathrm{bubble}}$, and
can be several tens of days for some models. This is because the width of the
shell becomes larger with $R_{\mathrm{bubble}}$, and it takes a longer time
for the shock to cross the shell where the density rises with radius. The peak
luminosity monotonically decreases with $R_{\mathrm{bubble}}$, because the
shell would have a lower density for larger $R_{\mathrm{bubble}}$ (see also
Figure 1).

In the lower left panel, the light curves display the most interesting
transition with $R_{\rm bubble}$. For small $R_{\rm bubble}$ the first peak is
powered pre-dominantly by the interaction, and the effect of the nickel decay
is generally buried by the CSM interaction. The 56Ni is seen as a late tail
after the shock completely sweeps the dense CSM and interaction ceases. For a
large $R_{\rm bubble}$, the first peak is instead powered by 56Ni decay. The
second peak emerges when the interaction starts to dominate over 56Ni decay
($\gtrsim 50\,{\rm days}$ after explosion for these parameters). The sharp
break at late time for large $R_{\rm bubble}$ is powered by interaction with
the part of the unshocked CSM whose density drops steeply with radius (the
$r>r_{*}$ part of equation 1).

As can be seen from the lower right panel, whether the second peak emerges
also depends on $E_{\mathrm{ej}}$. For larger $E_{\mathrm{ej}}$, the shock
front arrives at the CSM shell faster for a fixed $R_{\rm bubble}$. If the
epoch of CSM interaction is similar to when the luminosity from the
radioactive decay reaches its peak, CSM interaction can also bury the 56Ni
component.

We finally mention the dependence of the light curves on $t_{\mathrm{inj}}$,
plotted in Figure 3. In the left panels, we plot the dependence on
$t_{\mathrm{inj}}$ adopting $R_{\mathrm{bubble}}=4\times 10^{15}\,{\rm cm}$,
while we use $R_{\mathrm{bubble}}=10^{16}\,{\rm cm}$ in the right panels. As
can be seen from the figure, the peak luminosity and rise time vary with
$t_{\mathrm{inj}}$ when $R_{\mathrm{bubble}}$ is small, while they are
insensitive to $t_{\mathrm{inj}}$ if $R_{\mathrm{bubble}}$ is large. This is
because the mass of the detached shell $M_{\mathrm{shell}}$ is determined by
the mass of the continuous CSM enclosed within $R_{\mathrm{bubble}}$ (equation
8). For the profile of continuous CSM (equation 1), most of the CSM mass is
contained in the vicinity of $r_{*}$, where $r_{*}\propto t_{\rm inj}$. Thus
for fixed total CSM mass $M_{\rm CSM}$, the shell mass $M_{\mathrm{shell}}$ is
nearly constant ($\approx M_{\rm CSM}$) in the range
$r_{*}<R_{\mathrm{bubble}}$, while $M_{\mathrm{shell}}$ decreases with
$t_{\mathrm{inj}}$ for $r_{*}>R_{\mathrm{bubble}}$. The latter effect is seen
in the left panels, because $r_{*}\approx R_{\mathrm{bubble}}=4\times 10^{15}$
cm for $t_{\mathrm{inj}}\sim 1\,{\rm yr}$. For $t_{\mathrm{inj}}\gtrsim
1\,{\rm yr}$, $M_{\rm shell}$ decreases with increasing $t_{\rm inj}$, and
thus the peak luminosity gradually decreases with $t_{\rm inj}$. On the other
hand in the right panel ($R_{\mathrm{bubble}}=10^{16}\,{\rm cm}$),
$r_{*}<R_{\rm bubble}$ is always satisfied for $t_{\mathrm{inj}}=0.6$–$2\,{\rm
yr}$, and thus the rise time and peak luminosity are similar for all models.

Table 1: Estimated parameters for SNe 2010al and 2022xxf. SN name | Progenitor model | $M_{\mathrm{ej}}$ | $E_{\mathrm{ej}}$ | $f_{\mathrm{inj}}\,(M_{\mathrm{CSM}})$ | $t_{\mathrm{inj}}$ | $M_{\mathrm{Ni}}$ | $R_{\mathrm{bubble}}$ | $x_{\mathrm{comp}}$
---|---|---|---|---|---|---|---|---
|  | $[M_{\odot}]$ | $[10^{51}\,{\rm erg}]$ | $([M_{\odot}])$ | [yr] | $[M_{\odot}]$ | $[10^{15}\,{\rm cm}]$ |
SN 2010al | He star $25\,M_{\odot}$ | 4.8 | 1.2 | 0.8 (0.20) | 0.4 | 0 (assumed) | 1.5 | 2
SN 2022xxf | CO star $35\,M_{\odot}$ | 5.5 | 8 | 0.7 (0.24) | 1.2 | 0.6 | 28 | 3

##  3 Comparison with Observations

Figure 4: Comparison of our model with the observed light curves of SN 2010al
(top) and SN 2022xxf (bottom). For comparison, in each plot the corresponding
light curve with CSM of same mass that is continuous down to the star is
plotted as dashed lines, with the same explosion epoch as solid lines. The
dotted line in the bottom panel shows the case where the velocity of the
detached shell $v_{\rm shell}$ is reduced to 0.

In this section, we compare our results with the observational data of SN
2010al that has a slowly-rising ($\gtrsim 10\,{\rm days}$) phase in its light
curve (Pastorello et al., 2015a), and SN 2022xxf that has the second peak
suggested to be powered by the interaction with a massive CSM (Kuncarayakti et
al., 2023). The adopted progenitor models and best-fit parameters for these
SNe are summarized in Table 1.

As mentioned in Section 2, some samples of SNe Ibn with long initial rise have
been found (Pastorello et al., 2015a, b, 2016), which is difficult to explain
by diffusion in the CSM. The top panel of Figure 4 shows the comparison of our
model with one of the events SN 2010al. As seen from this figure, the detached
CSM can explain the slow rise of $\gtrsim 10\,{\rm day}$ with the moderate CSM
mass, while the light curve rapidly reaches its peak within $\sim 1\,{\rm
day}$ for a CSM with the same mass but continuous down to the star following
equation (1).

We note that the possibility of a flatter CSM profile at inner radii for SN
2010al was also suggested in Maeda & Moriya (2022). The radii where the flat
CSM emerges are estimated in their work to be $\lesssim 2\times 10^{15}$ cm
(their Figure 7), consistent with our results. If this detached CSM is
generated by ram pressure feedback from the progenitor’s stellar wind (Tsuna &
Takei, 2023), we speculate that this SN may have exploded from a progenitor
whose wind (just before core-collapse) was much stronger than progenitors of
fast-rising SN Ibn. A detached CSM of required mass and radius may also be
realized by binary interaction via common envelope ejection (Wu & Fuller,
2022, see their Figure 5), although there are considerable uncertainties in
the parameters that govern the ejection of the CSM.

The comparison of our model with the bolometric light curve of SN 2022xxf
(Appendix of Kuncarayakti et al., 2023) is plotted in the bottom panel. We can
successfully reproduce the first peak by the radioactive decay, and the second
by the interaction with a detached CSM shell at (a few) $\times 10^{16}$ cm.
The large explosion energy of $8\times 10^{51}\,{\rm erg}$ and massive
progenitor mass are consistent with the required nickel mass of
$0.6\,M_{\odot}$. These values for the SN ejecta are roughly in the range
measured for broad-lined Type Ic SNe (e.g., Taddia et al., 2019).

We note that the luminosity and duration of the second peak can slightly
change depending on the uncertain value of $v_{\rm shell}$. The dotted line is
the case for the same detached CSM but $v_{\rm shell}$ reduced to $0$, which
shows a brighter and slightly longer second peak. The observed velocity of
$1000$–$3000$ km s-1 for the detached CSM of SN 2022xxf (Kuncarayakti et al.,
2023) may indicate that the real light curve likely lies in between these two
assumptions.

Furthermore, we find that the wind bubble has to be very extended to power CSM
interaction until $\sim 100$ days after explosion. If the stellar wind is
responsible for sweeping the CSM to create this shell, the wind velocity
estimated from our CSM model of $t_{\rm inj}=1.2$ yr is $v_{\rm w}\sim
R_{\mathrm{bubble}}/t_{\mathrm{inj}}\sim 7000\,{\rm km\>s^{-1}}$. While
terminal wind speeds close to this value are observed in some galactic WO-type
Wolf-Rayet stars (Drew et al., 2004; Sander et al., 2012), the required wind
velocity is extreme. This may be alleviated if we consider different
progenitors and/or energy injections where the same amount of CSM mass can be
ejected with a lower velocity, which will make $t_{\rm inj}$ longer for the
same CSM density profile. Alternatively, the wind may be launched to such high
speeds by mechanisms other than radiation pressure, such as binary
interaction.

##  4 Summary and Discussion

In this work we modeled SNe Ibc interacting with CSM detached from the
progenitor, using the open-source code CHIPS. We showed that the detached CSM
can yield a wide variety of peculiar light curve features seen in interacting
SNe Ibc, such as the slow rise and/or multiple peaks. As a demonstration we
also reproduced the bolometric light curves of SN 2010al and SN 2022xxf. In
what follows, we discuss the possible caveats of our model, and avenues for
future modeling and observations.

The largest uncertainty of our model is the assumption of LTE in the shocked
region, which forces the temperatures of radiation and gas to be equal. This
requires the shocked gas to be efficiently cooled, which becomes more
difficult for detached CSM that have lower gas densities than confined CSM
with the same mass. Furthermore, the high-energy photons (mainly in X-rays)
emitted from the shocked gas has to be converted to optical, which is
generally more difficult for CSM with lower densities (Chevalier & Irwin,
2012; Svirski et al., 2012; Tsuna et al., 2021a).

If the gas in the shocked region can cool faster than the dynamical time, the
high-energy photons may be reprocessed to optical (Maeda & Moriya, 2022).
While the CSM for SN 2010al is still compact that this assumption holds (as
demonstrated in Maeda & Moriya 2022), the situation is less clear for SN
2022xxf with much larger CSM radii.

To diagnose whether this is true, we compare the cooling timescale
$t_{\mathrm{cool}}$ at when the CSM is crossing the shell, with the time that
is needed to cross the CSM shell of width $\Delta R_{\mathrm{bubble}}$,
$t_{\mathrm{cross}}$. These timescales are roughly estimated as
$t_{\mathrm{cool}}\approx U_{\mathrm{int}}/4\pi\eta^{\mathrm{ff}}$ and
$t_{\mathrm{cross}}\approx x_{\rm
comp}^{-1}R_{\mathrm{bubble}}/v_{\mathrm{sh}}$, where $U_{\mathrm{int}}$
denotes the internal energy density at the immediate downstream of the shock
front, and $\eta^{\mathrm{ff}}\,[{\rm erg\>s^{-1}\>cm^{-3}\>sr^{-1}}]$ is the
frequency-integrated free-free emissivity. Assuming that the CSM is the ideal
gas with adiabatic index of $5/3$ composed of fully-ionized oxygen,
$U_{\mathrm{int}},\,\eta^{\mathrm{ff}}$ can be expressed as (in cgs units)
(Rybicki & Lightman, 1979),

$\displaystyle U_{\mathrm{int}}=\frac{9}{32}\rho v_{\mathrm{sh}}^{2},$ (9)
$\displaystyle 4\pi\eta^{\mathrm{ff}}=7.3\times
10^{-25}n_{\mathrm{O}}^{2}T_{\mathrm{e}}^{1/2}\bar{g}_{\mathrm{B}},$ (10)

where $\rho,n_{\mathrm{O}}$ and $T_{\mathrm{e}}$ are respectively the density,
number density of oxygen and the electron temperature at the downstream, and
$\bar{g}_{\mathrm{B}}\sim 1$ is the frequency and velocity average of the
Gaunt factor. For the case of SN 2022xxf, substituting the above equations
into $t_{\mathrm{cool}}$ yields,

$\displaystyle t_{\mathrm{cool}}$ $\displaystyle\sim$ $\displaystyle 4.8\times
10^{3}\,{\rm day}\left(\frac{T_{\mathrm{e}}}{7\times 10^{8}\,{\rm
K}}\right)^{-1/2}$ (11) $\displaystyle\times$
$\displaystyle\left(\frac{\rho}{10^{-16}\,{\rm
g\>cm^{-3}}}\right)^{-1}\left(\frac{v_{\mathrm{sh}}}{2\times 10^{4}\,{\rm
km\>s^{-1}}}\right)^{2},$

where we derive $T_{\mathrm{e}}\approx 7\times 10^{8}\,{\rm K}$ assuming that
the time for energy transfer from ions to electrons by Coulomb relaxation is
limited to the duration of interaction (Chevalier & Fransson, 2006), here
$t_{\mathrm{cross}}\approx 60\,{\rm days}$. Since this cooling time is much
longer than $t_{\mathrm{cross}}$, it appears difficult for optical photons to
be emitted during the interaction.

The order-of-magnitude estimates above are based on a spherically symmetric
CSM with a homogeneous density structure. Clumpy structures caused by
Rayleigh-Taylor instability and/or radiative cooling mentioned in Section 2.2
can lead to locally enhanced gas cooling. This would result in a much more
efficient conversion of high-energy photons to optical bands, for the same CSM
mass and dissipation luminosity. Since the cooling time scales as $t_{\rm
cool}\propto\rho^{-2}$, the density of clumps for SN 2022xxf is required to be
$\gtrsim 10$ times higher than what would be for a homogeneous CSM.

A radiative shock due to highly clumped CSM is also suggested in a Type IIn SN
2005ip, to explain the long-lasting bright H$\alpha$ luminosity of $\gtrsim
10^{40}\,{\rm erg\>s^{-1}}$ for over 1000 days and the persistent
intermediate-width lines (Smith et al., 2009, 2017, see also Chugai & Danziger
1994 for SN 1988Z). Smith et al. (2009) inferred clumping in the CSM of a
factor of $\sim 10$ from the observed persistent intermediate-width lines of
SN 2005ip, which is similar to what is required for SN 2022xxf. While
discussing the detailed effects of clumping is beyond the scope of this work,
we plan to explore the detailed structure of the clumpy CSM in future work, by
studying wind-CSM interaction via multi-dimensional hydrodynamical simulations
(Tsuna & Huang, in prep.). We believe incorporating such modeling would
realize more accurate multi-wavelength light curve predictions from optical to
X-rays.

Finally, these interacting SNe with detached CSM can power bright non-thermal
emission, via particle acceleration across the collisionless shock between the
ejecta and CSM. Interaction with detached CSM is believed to be responsible
for bright radio sources years after core-collapse (e.g., Dong et al., 2021;
Stroh et al., 2021), and there is also an intriguing Type Ibn SN consistent
with being the optical counterpart of a high-energy neutrino event (Stein et
al., 2023). Recently, the Astrophysical Multimessenger Emission Simulator
(AMES) code has been developed to calculate multimessenger signals such as
radio, $\gamma$-rays, and high-energy neutrinos (Murase, 2018, 2023). In
future work we plan to incorporate the output from this work in AMES, to
predict diverse multi-wavelength and multi-messenger signals from interactions
with detached CSM.

## Acknowledgments

We thank Toshikazu Shigeyama for comments, and Kohta Murase for discussions.
DT is supported by the Sherman Fairchild Postdoctoral Fellowship at Caltech.

## References

  * Brethauer et al. (2022) Brethauer, D., et al. 2022, ApJ, 939, 105
  * Bruch et al. (2021) Bruch, R. J., et al. 2021, ApJ, 912, 46
  * Chevalier & Fransson (1992) Chevalier, R. A., & Fransson, C. 1992, ApJ, 395, 540
  * Chevalier & Fransson (2006) —. 2006, ApJ, 651, 381
  * Chevalier & Irwin (2012) Chevalier, R. A., & Irwin, C. M. 2012, ApJ, 747, L17
  * Chugai et al. (2004) Chugai, N. N., Blinnikov, S. I., Cumming, R. J., Lundqvist, P., Bragaglia, A., Filippenko, A. V., Leonard, D. C., Matheson, T., & Sollerman, J. 2004, MNRAS, 352, 1213
  * Chugai & Danziger (1994) Chugai, N. N., & Danziger, I. J. 1994, MNRAS, 268, 173
  * Das et al. (2023) Das, K. K., et al. 2023, arXiv e-prints, arXiv:2306.04698
  * Dong et al. (2021) Dong, D. Z., et al. 2021, Science, 373, 1125
  * Drew et al. (2004) Drew, J. E., Barlow, M. J., Unruh, Y. C., Parker, Q. A., Wesson, R., Pierce, M. J., Masheder, M. R. W., & Phillipps, S. 2004, MNRAS, 351, 206
  * Ferrari et al. (2024) Ferrari, L., et al. 2024, MNRAS, 529, L33
  * Gal-Yam et al. (2022) Gal-Yam, A., et al. 2022, Nature, 601, 201
  * Ho et al. (2023) Ho, A. Y. Q., et al. 2023, ApJ, 949, 120
  * Jermyn et al. (2023) Jermyn, A. S., et al. 2023, ApJS, 265, 15
  * Jun (1998) Jun, B.-I. 1998, ApJ, 499, 282
  * Kennel & Coroniti (1984) Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 694
  * Kleiser & Kasen (2014) Kleiser, I. K. W., & Kasen, D. 2014, MNRAS, 438, 318
  * Ko et al. (2022) Ko, T., Tsuna, D., Takei, Y., & Shigeyama, T. 2022, ApJ, 930, 168
  * Kuncarayakti et al. (2022) Kuncarayakti, H., et al. 2022, ApJ, 941, L32
  * Kuncarayakti et al. (2023) —. 2023, A&A, 678, A209
  * Kuriyama & Shigeyama (2020) Kuriyama, N., & Shigeyama, T. 2020, A&A, 635, A127
  * Maeda & Moriya (2022) Maeda, K., & Moriya, T. J. 2022, ApJ, 927, 25
  * Margutti et al. (2017) Margutti, R., et al. 2017, ApJ, 835, 140
  * Mauerhan et al. (2013) Mauerhan, J. C., et al. 2013, MNRAS, 430, 1801
  * Milisavljevic et al. (2015) Milisavljevic, D., et al. 2015, ApJ, 815, 120
  * Moriya (2015) Moriya, T. J. 2015, ApJ, 803, L26
  * Moriya (2023) —. 2023, MNRAS, 524, 5309
  * Moriya & Maeda (2016) Moriya, T. J., & Maeda, K. 2016, ApJ, 824, 100
  * Moriya et al. (2014) Moriya, T. J., Maeda, K., Taddia, F., Sollerman, J., Blinnikov, S. I., & Sorokina, E. I. 2014, MNRAS, 439, 2917
  * Murase (2018) Murase, K. 2018, Phys. Rev. D, 97, 081301
  * Murase (2023) —. 2023, arXiv e-prints, arXiv:2312.17239
  * Murase et al. (2014) Murase, K., Thompson, T. A., & Ofek, E. O. 2014, MNRAS, 440, 2528
  * Pastorello et al. (2015a) Pastorello, A., et al. 2015a, MNRAS, 449, 1921
  * Pastorello et al. (2008) —. 2008, MNRAS, 389, 113
  * Pastorello et al. (2007) —. 2007, Nature, 447, 829
  * Pastorello et al. (2016) —. 2016, MNRAS, 456, 853
  * Pastorello et al. (2015b) —. 2015b, MNRAS, 449, 1941
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) —. 2015, ApJS, 220, 15
  * Paxton et al. (2018) —. 2018, ApJS, 234, 34
  * Paxton et al. (2019) —. 2019, ApJS, 243, 10
  * Pellegrino et al. (2022) Pellegrino, C., et al. 2022, ApJ, 938, 73
  * Piro & Morozova (2014) Piro, A. L., & Morozova, V. S. 2014, ApJ, 792, L11
  * Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, Radiative processes in astrophysics
  * Sander et al. (2012) Sander, A., Hamann, W. R., & Todt, H. 2012, A&A, 540, A144
  * Sander & Vink (2020) Sander, A. A. C., & Vink, J. S. 2020, MNRAS, 499, 873
  * Smith et al. (2017) Smith, N., et al. 2017, MNRAS, 466, 3021
  * Smith et al. (2009) —. 2009, ApJ, 695, 1334
  * Stein et al. (2023) Stein, R., et al. 2023, GRB Coordinates Network, 34837, 1
  * Stroh et al. (2021) Stroh, M. C., et al. 2021, ApJ, 923, L24
  * Strotjohann et al. (2021) Strotjohann, N. L., et al. 2021, ApJ, 907, 99
  * Sun et al. (2020) Sun, N.-C., Maund, J. R., & Crowther, P. A. 2020, MNRAS, 497, 5118
  * Suzuki & Maeda (2021) Suzuki, A., & Maeda, K. 2021, ApJ, 908, 217
  * Svirski et al. (2012) Svirski, G., Nakar, E., & Sari, R. 2012, ApJ, 759, 108
  * Taddia et al. (2019) Taddia, F., et al. 2019, A&A, 621, A71
  * Takei & Shigeyama (2020) Takei, Y., & Shigeyama, T. 2020, PASJ, 72, 67
  * Takei et al. (2024) Takei, Y., Tsuna, D., Ko, T., & Shigeyama, T. 2024, ApJ, 961, 67
  * Takei et al. (2022) Takei, Y., Tsuna, D., Kuriyama, N., Ko, T., & Shigeyama, T. 2022, ApJ, 929, 177
  * Tsuna et al. (2021a) Tsuna, D., Kashiyama, K., & Shigeyama, T. 2021a, ApJ, 914, 64
  * Tsuna et al. (2024) Tsuna, D., Matsumoto, T., Wu, S. C., & Fuller, J. 2024, arXiv e-prints, arXiv:2401.02389
  * Tsuna & Takei (2023) Tsuna, D., & Takei, Y. 2023, PASJ, 75, L19
  * Tsuna et al. (2021b) Tsuna, D., Takei, Y., Kuriyama, N., & Shigeyama, T. 2021b, PASJ, 73, 1128
  * Woosley (2019) Woosley, S. E. 2019, ApJ, 878, 49
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015
  * Wu & Fuller (2022) Wu, S. C., & Fuller, J. 2022, ApJ, 940, L27


# Dependence and mixing for perturbations of copula-based Markov chains

Martial Longla<EMAIL_ADDRESS>University of Mississippi, Department of
mathematics  Mathias Muia Nthiani<EMAIL_ADDRESS>University of
Mississippi, Department of mathematics  Fidel Djongreba Ndikwa
<EMAIL_ADDRESS>University of Maroua, Department of mathematics

###### Abstract

This paper explores the impact of perturbations of copulas on dependence
properties of the Markov chains they generate. We use an observation that is
valid for convex combinations of copulas to establish sufficient conditions
for the mixing coefficients $\rho_{n}$, $\alpha_{n}$ and some other measures
of association. New copula families are derived based on perturbations of
copulas and their multivariate analogs for $n$-copulas are provided in
general. Several families of copulas can be constructed from the provided
framework.

Key words: Perturbation, central limit theorem, Copulas, Mixing rates,
Dependence coefficients.

Mathematical Subject Classification (2000): 62G08, 62M02, 60J35

##  1 Introduction

Many authors have contributed to constructing new families of copulas to model
dependence among variables or factors in economics, finance, risk management
and other applied fields. Following the ideas of Durante and al. (2013), we
consider the perturbation method that adds to a copula an extra term called
perturbation. This modification is often used to increase flexibility of the
dependence structure of the considered models. Here, we look at other classes
of modifications and their impact on the dependence structure as studied by
Komornik and al. (2017). We consider the long run impact of these
perturbations on the dependence structure and the measures of association. The
case is presented for $\rho$-mixing and $\alpha$-mixing. This is a
continuation of the work done by Longla and al. (2021) in which $\psi$-mixing,
$\beta$-mixing and $\phi$-mixing were explored. The paper is structured as
follows. Introduction in Section 1 is divided into several parts, each of
which is dedicated to a specific topic of interest. Copulas are introduced in
Section 1.1, perturbations of copulas are defined in Section 1.2 with some
results on perturbations, Section 1.3 is dedicated to the perturbation of
components of a vector and their impact on the copulas, Section 1.4 is
dedicated to the mixing coefficients of interest, Section 1.5 visits the
impact of perturbations on some dependence coefficients (measures of
association). In Section 2, we present the impact of perturbations on the
mixing structure of copula-based Markov chains, addressing $\alpha-$mixing in
Section 2.1 and $\rho$-mixing in Section 2.2. In Section 3, we provide the
proofs of the main results of this work.

###  1.1 Copula theory

The importance of copulas is due to the fact that they are tools that link
univariate distributions to the joint distribution of the random variables.
They were first presented in Sklar (1959), and over the recent years they have
become very important for analyzing temporal dependence of time series. For
$n\geq 2$, an $n$-copula is an $n$-variate distribution function on
$[0,1]^{n}$ with uniform marginals. One of the most used types of copula is
the $2$-copula, that is mostly used for copula-based Markov chains. By
definition, a $2$-copula is a bivariate distribution function that has uniform
marginal distributions on $I=[0,1]$. They have gained a lot of popularity in
financial, economics and risk management time series because they are used to
generate Markov chains with important needed properties. For example,
Ibragimov (2009) illustrates applications of copulas in Econometric theory.
The definition of a $2$-copula and related topics can be found in Nelsen
(2006). Formally, a function $C:[0,1]^{2}\rightarrow[0,1]$ is called bivariate
copula if it satisfies the following conditions:

  1. i.

$C(0,x)=C(x,0)=0$ (meaning that $C$ is grounded);

  2. ii.

$C(x,1)=C(1,x)=x,\forall x\in I$ (meaning each coordinate is uniform on $I$);

  3. iii.

$C(a,c)+C(b,d)-C(a,d)-C(b,c)\geq 0,\forall\ [a,b]\times[c,d]\subset I^{2}.$

Darsaw and al. (1992) derived the transition probabilities for stationary
Markov chains with uniform marginals on $[0,1]$ as $P(X_{n}\in
A|X_{n-1}=x)=C_{,1}(x,y),\forall n\in\mathbb{Z}$ and
$A=\left(-\infty,y\right]$. Here, $C_{,i}(x,y)$ denotes the derivative with
respect to the $i^{th}$ variable. This property has been used by many authors
to establish mixing properties of copula-based Markov chains. We can cite
Longla (2015), Longla (2014), Longla and Peligrad (2012) who provided some
results for reversible Markov chains and Beare (2010) who presented results
for $\rho$-mixing. We are assuming in the sequel that variables in a copula-
based Markov chain have the uniform distribution on $[0,1]$.

It’s been shown in the literature (see Darsow and al. (1992) and the
references therein) that if $(X_{1},\cdots,X_{n})$ is a Markov chain with
consecutive copulas $(C_{1},\cdots,C_{n-1})$, then for any integer $k$, the
fold product given by

$C(x,y)=C_{k}*C_{k+1}(x,y)=\int^{1}_{0}C_{k,2}(x,t)C_{k+1,1}(t,y)dt$

is the copula of $(X_{k},X_{k+2})$ and the $\star$-product given by

$C(x,y,z)=C_{k}\star C_{k+1}(x,y,z)=\int_{0}^{y}C_{k,2}(x,t)C_{k+1,1}(t,z)dt$

is the copula of $(X_{k},X_{k+1},X_{k+2})$. The $n$-fold product of $C(x,y)$
denoted $C^{n}(x,y)$ is defined by the recurrence $C^{1}(x,y)=C(x,y)$,
$C^{n}(x,y)=C^{n-1}*C(x,y),\quad\text{for $n>1.$}$ This product is the joint
distribution of $(X_{k},X_{k+n})$ when all copulas are equal along the chain.
The most popular copulas are $\Pi(u,v)=uv$ (the independent copula), the
Hoeffding lower and Upper bounds $W(u,v)=\max(u+v-1,0)$ and $M(u,v)=\min(u,v)$
respectively. Convex combinations of copulas
$\\{C_{1}(x,y),\cdots,C_{k}(x,y)\\}$ defined by
$\displaystyle\\{C(x,y)=\sum_{j=1}^{k}a_{j}C_{j}(x,y),0\leq
a_{j},\sum_{j=1}^{k}a_{j}=1\\}$ are also copulas and the following holds.

###### Proposition 1.1.1

convex combinations and fold product

  1. 1.

The fold product of convex combinations of copulas is a convex combination of
copulas.

  2. 2.

The $n$-fold product of any convex combination of copulas is a convex
combination of copulas.

The proof of this proposition relies on the fact that the fold product defines
a bilinear operator. For any copulas
$A_{1}(x,y),A_{2}(x,y),B_{1}(x,y),B_{2}(x,y)$ be copulas. It is easy to verify
that

$A_{1}*(a_{1}B_{1}+a_{2}B_{2})(x,y)=a_{1}A_{1}*B_{1}(x,y)+a_{2}A_{1}*B_{2}(x,y)$
$(a_{1}B_{1}+a_{2}B_{2})*A_{1}(x,y)=a_{1}B_{1}*A_{1}(x,y)+a_{2}B_{2}*A_{1}(x,y)$

So,
$(a_{1}A_{1}+a_{2}A_{2})*(b_{1}B_{1}+b_{2}B_{2})(x,y)=a_{1}b_{1}A_{1}*B_{1}(x,y)+a_{1}b_{2}A_{1}*B_{2}(x,y)+a_{2}b_{1}A_{2}*B_{1}(x,y)+a_{2}b_{2}A_{2}*B_{2}(x,y).$
if all initial coefficients are positive with $a_{1}+a_{2}=1$ and
$b_{1}+b_{2}=1$, then the last set of coefficients is made of positive numbers
with
$a_{1}b_{1}+a_{1}b_{2}+a_{2}b_{1}+a_{2}b_{2}=a_{1}(b_{1}+b_{2})+a_{2}(b_{1}+b_{2})=1$.
Thus, the fold product is a convex combination. The second part follows by
multiplying the convex combination by itself. This result by itself is
elementary, but plays a big role in establishing many results or observations
in the sequel. From the proof of Proposition 1.1.1, it follows that for a
convex combination $\displaystyle C(x,y)=\sum_{i=1}^{k}a_{i}C_{i}(x,y)$ and
any integer $n$ the following formula holds.

$C^{n}(x,y)=\sum_{j=1}^{k^{n}}{b_{j}}\times{{}_{1}C_{j}}*\cdots*{{}_{n}C_{j}(x,y)},$
(1.1)

where $\displaystyle\sum_{j=1}^{k^{n}}b_{j}=1,\quad b_{j}\geq 0$ , each of the
copulas ${}_{i}C_{j}(x,y)=C_{j_{i}}(x,y)$ for some $j_{i}\in\\{1,\cdots,k\\}$
and the sum is over all possible products of $n$ copulas selected from the
original $k$ copulas with replacement. The notation ${}_{i}C_{j}$ indicates
that the copula $C_{j_{i}}$ was selected in the given $j^{th}$ element of
$B=\\{C_{1},\cdots,C_{k}\\}^{n}$. If we denote $c^{n}(x,y)$ the density of the
copula $C^{n}(x,y)$, then the following holds.

###### Lemma 1.1.2

If a copula $C(x,y)$ is such that its density $c(x,y)$ satisfies $c(x,y)\geq
K$ on a set of Lebesgue measure 1 for some real number $K>0$, then for all
positive integer $n$ we have $c^{n}(x,y)\geq K^{n}$ and $C^{n}(x,y)\geq
xyK^{n}$.

###  1.2 Perturbations of copulas

Copulas are well known tools for understanding dependence among random
variables. They offer a choice of appropriate models for dependence between
random variables independently from selection of marginal distributions. Often
in model analysis it is necessary to understand the dependence structure of
random variables for inference purposes. Commonly used measures of dependence
include Kendall’s tau ($\tau$), Spearman’s rho ($\rho$), Blomqvist’s beta
($\beta$), Gini’s gamma ($\gamma$), among others. Sometimes linear correlation
fails to be reliable in measuring dependence between random variables;
especially where non-linear monotonic transformations of random variables are
involved. It is at this point that the concordance and dependence measures
(e.g. Kendal’s $\tau$, Spearman’s $\rho$) perform best as they are invariant
under monotonic transformations. We consider the copula $C$ and the related
copula $C_{H}$ defined by $C_{H}(u,v)=C(u,v)+H(u,v),$ where
$H:[0,1]^{2}\to\mathbb{R}$ is a continuous function. The function $H$ is the
perturbation factor and the copula $C_{H}$ is called a perturbation of $C$.
This perturbation method introduced by Durante and al. (2013) and its special
case with

$H(x,y)=H_{\theta}(u,v):=\theta(u-C(u,v))(v-C(u,v)),\theta\in[0,1]$

was investigated by Mesiar and al. (2015). We consider the class of
perturbations defined in general by

$C_{\theta}(u,v)=C(u,v)+\theta[C_{1}(u,v)-C(u,v)],\quad\theta\in[0,1],\quad
C_{1}\mbox{ is a copula,}$ (1.2)

then classes of perturbations with $C_{1}(u,v)=\Pi(u,v)=uv$ and
$C_{1}(u,v)=M(u,v)$ are respectively as follows

$C_{\theta,\Pi}(u,v)=C(u,v)+\theta[uv-C(u,v)]=(1-\theta)C(u,v)+\theta\Pi(u,v),\quad\theta\in[0,1],$
(1.3) $C_{\theta,M}(u,v)=C(u,v)+\theta[M(u,v)-C(u,v)]=(1-\theta)C(u,v)+\theta
M(u,v),\quad\theta\in[0,1].$ (1.4)

This class of perturbations modifies all copulas except when $C=C_{1}$.
Technically, perturbations of this kind happen when the bivariate population
of interest with copula $C(u,v)$ has been contaminated by another population
with copula $C_{1}(u,v)$. We consider the long run effect of the perturbations
given by (1.3) and (1.4). By Longla and al. (2021), we have

$C^{n}_{\theta,M}(u,v)=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}C^{i}(u,v)+\theta^{n}M(u,v)$ (1.5)
$C^{n}_{\theta,\Pi}(u,v)=(1-\theta)^{n}C^{n}(u,v)+(1-(1-\theta)^{n})\Pi(u,v).$
(1.6)

Based formulas (1.5) and (1.6), via simple computations, the following is
established.

###### Proposition 1.2.1

For any $\theta\in(0,1)$,

  1. 1.

If $\displaystyle\lim_{n\to\infty}C^{n}(u,v)=C_{o}(u,v)$, then
$\displaystyle\lim_{n\to\infty}C^{n}_{\theta,M}(u,v)=C_{o}(u,v)$

  2. 2.

$\displaystyle\lim_{n\to\infty}C^{n}_{\theta,\Pi}(u,v)=\Pi(u,v).$

The proof of Proposition 1.2.1 relies on Proposition 1.0.1 of Longla and al.
(2021). It makes more sense looking at this proposition to consider formula
(1.6) as a perturbation of the independence copula by means of the copula
$C(u,v)$. As an example, we can take a copula $C(u,v)$ from the two parameters
Mardia family of copulas (a subfamily of this class of copulas is the one
parameter Frechet copula family ). Copulas from this family are defined as
follows.

$C=aM(u,v)+bW+(1-a-b)\Pi,\quad\mbox{for some}\quad a,b,a+b\in[0,1].$ (1.7)

When $a=\theta^{2}(1+\theta)/2$ and $b=\theta^{2}(1-\theta)/2$ the Mardia
copula turns into a one parameter copula called Frechet copula. Longla (2014)
has shown that for this class of copulas,

$C^{n}=\frac{1}{2}[(a+b)^{n}+(a-b)^{n}]M+\frac{1}{2}[(a+b)^{n}-(a-b)^{n}]W+[1-(a+b)^{n}]\Pi.$

Therefore,

$\lim_{n\to\infty}C^{n}(u,v)=\begin{cases}\frac{1}{2}M(u,v)+\frac{1}{2}W(u,v),&\mbox{if}\quad
a+b=1,a\neq 1,b\neq 1\\\ \mbox{Undefined},&\mbox{if}\quad b=1,a=0,\\\
M(u,v),&\mbox{if}\quad a=1,b=0,\\\ \Pi(u,v),&\mbox{if}\quad a+b<1.\end{cases}$
(1.8)

Convergence in this case is uniform convergence. This can be obtained via
simple computations that we omit here. This limit by itself shows that copulas
from the Mardia family generate Markov chains for which $(X_{0},X_{n})$ are
either asymptotically independent (when $a+b<1$), or asymptotically almost
surely equal (when $a=1,b=0$) or asymptotically almost surely $X_{n}=X_{0}$
with probability $1/2$ or $X_{n}=1-X_{0}$ with probability $1/2$ (when
$a+b=1,a\neq 1,b\neq 1$). For $b=1,a=0$ there is no limit of the copula
sequence. This limit also shows that Frechet copulas either converge to
$\frac{1}{2}M(u,v)+\frac{1}{2}W(u,v)$ (for $|\theta|<1$) or converge to
$M(u,v)$ (for $\theta=1$) or have no limit when $\theta=-1$. The last case is
because $W^{2}(u,v)=M(u,v)$, and $M*W(u,v)=M*W(u,v)=W(u,v)$.

As mentioned by Durante and al. (2013) perturbations of copulas have been
studied under simplified assumptions on $C$ and $H$. In particular, in the
literature it is often assumed that $C(u,v)=\Pi(u,v).$ The model is
interpreted as a contamination of independence. We will also consider in this
paper another type of perturbation that was studied in Longla and al. (2021)
and was considered for measures of association in the works of Durante and al
(2013), Mesiar and al (2015) and Komornik and al (2017). This perturbation
method acts at the level of the random variables, polluting the underlying
variable $(X)$ with some noise variable $(Z)$. In the cited works, the
variable $X$ is mostly bivariate and the variable Z is also bivariate,
independent of $X$ and with possible dependence between its components. In
this work, we will consider a vector $X=(X_{1},\cdots,X_{k})$ and a
perturbation vector $Z=(Z_{1},\cdots,Z_{k})$ for some positive integer $k$.
The goal is to derive the copula of $X+Z$ under various conditions on $X$ and
$Z$ and provide some examples for $X$ having copula or distribution
$M(u_{1},\cdots,u_{k})$ or $\Pi(u_{1},\cdots,u_{k})$.

###  1.3 Perturbation of $X$ via $X+Z$ and its copulas

Assume we have a random $k$-dimensional vector $X=(X_{1},\cdots,X_{k})$ of
copula $C(u)$ and marginal distributions $F_{1},\cdots,F_{k}$. Consider a
noise vector $Z=(Z_{1},\cdots,Z_{k})$ independent of $X$. We are interested in
extracting the copula of $X+Z$ from its joint distribution function. Assume
$\tilde{F}_{k}$ is the cumulative distribution function of the random variable
$X_{k}+Z_{k}$ and $u=(u_{1},\cdots,u_{k})$.

###### Theorem 1.3.1

If $Z=(Z_{1},\cdots,Z_{s},0\cdots,0)$ for $s\leq k$ and $(Z_{1},\cdots,Z_{s})$
has independent components with marginal distributions $G_{1},\cdots,G_{s}$,
then the copula of $X+Z$ is

$C_{1}(u)=\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}H(u,t)dG_{1}(t_{1})\cdots
dG_{s}(t_{s}),$ (1.9)

where
$H(u,t)=C(F_{1}(\tilde{F}_{1}^{-1}(u_{1})-t_{1}),\cdots,F_{s}(\tilde{F}_{s}^{-1}(u_{s})-t_{s}),u_{s+1},\cdots,u_{k})$.

A consequence of Theorem 1.3.1 is the following:

###### Corollary 1.3.2

If $Z=(Z_{1},\cdots,Z_{k})$ has independent components with marginal
distributions $G_{1},\cdots,G_{k}$, then the copula of $X+Z$ is

$C_{2}(u)=\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}C(F_{1}(\tilde{F}_{1}^{-1}(u_{1})-t_{1}),\cdots,F_{k}(\tilde{F}_{k}^{-1}(u_{k})-t_{k}))dG_{1}(t_{1})\cdots
dG_{k}(t_{k}).$ (1.10)

###### Theorem 1.3.3

If $Z=(W,\cdots,W,0\cdots 0)$ where $W$ has distribution $G$ and $Z_{i}=0$ for
all $i>s$, then the copula of $X+Z$ is

$C_{3}(u)=\int_{-\infty}^{\infty}C(F_{1}(\tilde{F}_{1}^{-1}(u_{1})-t),\cdots,F_{s}(\tilde{F}_{s}^{-1}(u_{s})-t),u_{s+1},\cdots,u_{k})dG(t).$
(1.11)

A consequence of Theorem 1.3.3 is the following

###### Corollary 1.3.4

If $Z=(Z_{1},\cdots,Z_{1})$ where $Z_{1}$ has distribution $G$, then the
copula of $X+Z$ is

$C_{4}(u)=\int_{-\infty}^{\infty}C(F_{1}(\tilde{F}_{1}^{-1}(u_{1})-t),\cdots,F_{k}(\tilde{F}_{k}^{-1}(u_{k})-t))dG(t).$
(1.12)

The proofs of these statements are based on simple computations using
definitions of probabilities via conditioning an Sklar’s theorem as in Longla
and al. (2021).

###### Remark 1.3.5

Notice that in each of the two theorems, the position of the variables that
are modified doesn’t change the details of the proof. In fact, if
$\sigma:\\{1,\cdots,n\\}\to\\{1,\cdots,n\\}$ is a permutation such that
$\sigma(u)=v$, then the copula of $X_{\sigma}=X+Z_{\sigma}$, where components
of $Z$ have been permuted via $\sigma$ to obtain $Z_{\sigma}$ is

$C_{i,\sigma}(u)=C_{i}(v),\quad i=1,3.$

###  1.4 Mixing coefficients

The concept of mixing is important in the study of statistical properties of
deterministic dynamical systems (ergodic theory). The most important mixing
coefficients include but are not limited to $\phi$-mixing, $\rho$-mixing,
$\beta$-mixing, $\alpha$-mixing, $\rho*$-mixing, $\psi^{\prime}$-mixing and
$\psi$-mixing. The definitions of these can be found in Bradley (2007). Other
works on mixing include but are not limited to Miller (1994), Miller (1995)
Gaposhkin (1991) and the references therein. Bradley (1997) presents a proof
that every $\psi^{\prime}$-mixing Markov chain is $\rho^{*}$-mixing. He also
includes a detailed list of relationships between different mixing
coefficients. Longla (2015) showed that Markov chains generated by copulas are
lower $\psi$-mixing when the density of the absolutely continuous part of the
generating copula is bounded away from zero on a set of Lebesgue measure 1.
Longla and Peligrad (2012) defined the mixing coefficients for an absolutely
continuous copula and an absolutely continuous invariant distribution for the
states of a stationary Markov chain that they generate. They provided a result
on absolute regularity for mixtures of copulas. Archimedean copulas with non-
strict generators have been studied by Longla (2014). A condition was
established for exponential $\rho$-mixing of copula-based Markov chains in the
case of non-strict Archimedean copulas. A non-strict Archimedean copula is
defined by some convex function $\varphi:[0,1]\to[0,1]$ such that
$\varphi^{-1}(u)=0\quad\forall u>1$ via
$C(u,v)=\varphi^{-1}(\varphi(u)+\varphi(u)).$ The $\alpha$-mixing condition is
also known as strong mixing condition in the literature (see Bradley (2007)
and Bradley (2005)). We borrow the following historical background from
Bradley (2005): The strong mixing condition was introduced by Rosenblatt
(1956). The $\psi^{\prime}$-mixing condition was introduced by Bradley (1983),
following the work of Blum and al. (1963) on a related coefficient
($\psi$-mixing, not considered here). The $\rho$-mixing condition was
introduced by Kolmogorov and Rozanov (1960). These coefficients are defined in
the literature as follows (see Bradley (2007))

$\alpha(\mathscr{A},\mathscr{B})=\sup_{B\in\mathscr{B},A\in\mathscr{A}}|P(A\cap
B)-P(A)P(B)|,$ $\rho(\mathscr{A},\mathscr{B})=\sup_{f\in
L^{2}(\mathscr{A}),g\in
L^{2}(\mathscr{B})}corr(f,g),\quad\quad\psi^{\prime}(\mathscr{A},\mathscr{B})=\inf_{B\in\mathscr{B},A\in\mathscr{A}}\frac{P(B\cap
A)}{P(A)P(B)},$

where $\mathscr{A}=\sigma(X_{i},i\leq 0)$, $\mathscr{B}=\sigma(X_{i},i\geq n)$
and $P$ is the defined probability measure. For Markov chains generated by an
absolutely continuous copula (see Longla (2013) or Longla (2015)) these
coefficients are

$\psi^{\prime}_{n}:=\psi^{\prime}(\mathscr{A},\mathscr{B})=\inf_{A,B\in\mathscr{B}}\frac{|\int_{A}\int_{B}c_{n}(x,y)dydx|}{\lambda(A)\lambda(B)},$
$\rho_{n}=\sup_{f,g}\\{\int^{1}_{0}\int^{1}_{0}c_{n}(x,y)f(x)g(y)dxdy:||g||_{2}=||f||_{2}=1,\mathbb{E}(f)=\mathbb{E}(g)=0\\},$
$\alpha_{n}:=\alpha(\mathscr{A},\mathscr{B})=\sup_{A,B\in\mathscr{B}}|\int_{A}\int_{B}(c_{n}(x,y)-1)dydx|,$

where $c_{n}$ is the copula density of the random variable $(X_{0},X_{n})$,
$\mathscr{B}$ is the Borel sigma field and $\lambda$ is the Lebesgue measure
on $(0,1)$. A stochastic process is said to be an $\alpha$-mixing, if
$\alpha_{n}\rightarrow 0$; $\psi^{\prime}$-mixing if
$\psi^{\prime}_{n}\rightarrow 1$ and $\rho$-mixing if $\rho_{n}\rightarrow 0$.
The process is exponentially mixing, if the convergence rate is exponential.

###  1.5 Perturbations of copulas and dependence coefficients

If $C$ is the copula of $(X,Y)$, then the Spearman’s rho is
$\displaystyle\rho_{S}(C)=12\int_{0}^{1}\int_{0}^{1}C(u,v)dudv-3;$ Kendall’s
$\tau$ is
$\displaystyle\tau(C)=1-4\int_{0}^{1}\int_{0}^{1}C_{,1}(u,v)C_{,2}(u,v)dudv;$
Blomqvist’s $\beta$ is $\displaystyle\beta(C)=4C(1/2,1/2)-1$ and is called the
medial correlation coefficient. The Gini’s $\gamma$ is
$\displaystyle\gamma(C)=4\int_{0}^{1}C(u,u)+C(u,1-u)du;$ the upper tail
dependence coefficient is $\displaystyle\lambda_{U}(C)=\lim_{u\to
1}\frac{1-2u+C(u,u)}{1-u}$ and the lower tail dependence is
$\displaystyle\lambda_{L}(C)=\lim_{u\to 0}\frac{C(u,u)}{u}$ (see Nelsen (2006)
for more on these coefficients). From formulas (1.5) and (1.6), it follows:

###### Theorem 1.5.1

For any copula-based Markov chain generated by a copula $C(u,v)$,

  1. 1.

$\displaystyle\rho_{S}(C^{n}_{\theta,M})=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}\rho_{S}(C^{i})+\theta^{n}$ and
$\rho_{S}(C^{n}_{\theta,\Pi})=(1-\theta)^{n}\rho_{S}(C^{n})$.

  2. 2.

$\displaystyle\beta(C^{n}_{\theta,M})=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}\beta(C^{i})+\theta^{n}$ and
$\beta(C^{n}_{\theta,\Pi})=(1-\theta)^{n}\beta(C^{n})$.

  3. 3.

$\displaystyle\gamma(C^{n}_{\theta,M})=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}\gamma(C^{i})+3\theta^{n}$ and
$\gamma(C^{n}_{\theta,\Pi})=(1-\theta)^{n}\gamma(C^{n})+2(1-(1-\theta)^{n})$.

  4. 4.

$\displaystyle\lambda_{L}(C^{n}_{\theta,M})=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}\lambda_{L}(C^{i})+\theta^{n}$ and
$\lambda_{L}(C^{n}_{\theta,\Pi})=(1-\theta)^{n}\lambda_{L}(C^{n})$.

  5. 5.

$\displaystyle\lambda_{U}(C^{n}_{\theta,M})=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}\lambda_{U}(C^{i})+\theta^{n}$ and
$\displaystyle\lambda_{U}(C^{n}_{\theta,\Pi})=(1-\theta)^{n}\lambda_{U}(C^{n})$.

If we denote $\xi(C)$ any of the considered measures of association other than
$\lambda_{L}$ or $\lambda_{U}$, then the following holds when $\theta\neq 1$.

###### Remark 1.5.2

Theorem 1.5.1 shows that as $n\to\infty$,
$\xi(C^{n}_{\theta,\Pi})\to\xi(\Pi)$. When the independence copula is
perturbed by any copula, the measure of association is asymptotically
unchanged. It can also be seen that $\xi(C^{n}_{\theta,M})$, as $n\to\infty$,
depends solely on properties of $C(u,v)$. Using Longla (2021), we can conclude
that the limit is equal to $\lim\xi(C^{n})$ when the later exists.

##  2 Perturbations of copulas and mixing coefficients

###  2.1 On $\alpha$-mixing for copula-based Markov chains

The convergence rate of this coefficient is related to the central limit
theorem via several conditions in the literature (see Merlevede and Peligrad
(2020) and the references therein). From Lemma 3 of Longla and Peligrad
(2012), the following holds.

###### Corollary 2.1.1

Any convex combination of copulas, each of which generates $\alpha$-mixing
Markov chains, generates $\alpha$-mixing Markov chains.

Following Corollary 2.1.1, after deriving the $n$-th transition kernel for
copula-based Markov chains generated by the perturbation copula, we can
establish the following.

###### Theorem 2.1.2

When a copula is perturbed using the independence copula $\Pi(u,v)$, the
$\alpha$-mixing rate of convergence is multiplied by an exponential factor.
For a perturbation that uses the copula $M(u,v)$,
$\alpha_{n}({C_{\theta,M}})\to 0$ when the initial copula generates
$\alpha$-mixing Markov chains. In a case when $\lim\alpha_{n}(C)=c\neq 0$,
$C_{\theta,M}$ doesn’t generate $\alpha$-mixing.

###### Remark 2.1.3

As strange as it might sound, this theorem shows that a perturbation of a non
$\alpha$-mixing generating copula by means of an $\alpha$-mixing generating
copula turns it into an $\alpha$-mixing generating copula. In this specific
case, the limit of $C_{\theta,M}^{n}(u,v)$ is actually not influenced by
$M(u,v)$ for any $\theta\in(0,1)$ and $C_{\theta,\Pi}^{n}(u,v)$ converges to
$\Pi(u,v)$. This also means that in the case of the perturbation of the
independence copula, long run realizations of the Markov chain generated by
the perturbation copula are independent of the initial state.

The following result for ergodic sequences is Theorem 3.4. of Bradley (2005).

###### Theorem 2.1.4

Suppose $(X_{n},n\in Z)$ is a strictly stationary Markov chain which is
ergodic and aperiodic. If $\alpha_{n}<1/4$ for some $n\geq 1$, then
$\alpha_{n}\to 0$ (but not necessarily exponentially fast) as $n\to\infty$.

Based on Theorem 2.1.4, the following holds for convex combinations of
copulas.

###### Theorem 2.1.5

For any set of copulas $C_{1}(u,v)\cdots C_{k}(u,v)$, if there exists a subset
of copulas $C_{k_{1}}\cdots C_{k_{s}},$ $s\leq k\in\mathbb{N}$ such that
$\alpha(\hat{C})<1/4\quad\text{for}\quad\hat{C}=C_{k_{1}}*\cdots*C_{k_{s}},$
then $\alpha_{s}(C)<1/4$ and any strictly stationary ergodic and aperiodic
Markov chain generated by $C=a_{1}C_{1}+\cdots+a_{k}C_{k}\quad\text{for }\quad
0<a_{1},\dots,a_{k}<1\quad\text{is}\quad\alpha-\text{mixing}.$

###  2.2 On $\rho$-mixing coefficient

Bradley (2005) states the following in his Theorem 3.3:

###### Theorem 2.2.1

Suppose $(X_{n},n\in Z)$ is a (not necessarily stationary) Markov chain. Then
each of the following statements holds:

If $\rho_{n}<1$ for some $n\geq 1$, then $\rho_{n}\to 0$ at least
exponentially fast as $n\to\infty$.

Theorem 2.2.1 requires $\rho_{n}<1$, therefore the following is easily
established.

###### Theorem 2.2.2

For any set of copulas $C_{1}(u,v)\cdots C_{k}(u,v)$, if there exists a subset
of copulas $C_{k_{1}}\cdots C_{k_{s}},$ $s\leq k\in\mathbb{N}$ such that
$\rho(\hat{C})<1\quad\text{for}\quad\hat{C}=C_{k_{1}}*\cdots*C_{k_{s}},$ then
$\rho_{s}(C)<1$ and any Markov chain generated by
$C=a_{1}C_{1}+\cdots+a_{k}C_{k}\quad\text{for }\quad
0<a_{1},\dots,a_{k}<1\quad\text{is exponential}\quad\rho-\text{mixing}.$

##  3 Proofs, conclusions and remarks

###  3.1 Proof of Lemma 1.1.2

Without loss of generality, we can assume that $(X,Y)$ has cumulative
distribution $C(x,y)$. This means that $X$ and $Y$ are uniform on $(0,1)$.
Thus,

$P(X_{0}<x,X_{n}<y)=C^{n}(x,y)=C^{n-1}*C(x,y).$

If by induction we assume $c^{n-1}(x,y)\geq K^{n-1}$ on a set of Lebesgue
measure 1, then $c^{1}(x,y)=c(x,y)\geq K$ and
$C^{n}(x,y)=\int_{0}^{1}C^{n-1}_{,2}(x,t)C_{,1}(t,y)dt$ imply

$c^{n}(x,y)=\int_{0}^{1}c^{n-1}(x,t)c(t,y)dt\geq\int_{0}^{1}K^{n-1}Kdt=K^{n}.$

Integrating both sides of this inequality gives $C^{n}(x,y)\geq K^{n}xy$.

###  3.2 Proof of Theorem 1.3.1

It is enough to conduct this proof for $k=3$ and $s=2$. For random variables
$X_{1},X_{2},X_{3}$ with distributions $F_{1},F_{2},F_{3}$, joint distribution
$F$ and copula $C(u,v,w)$, let $\tilde{F}_{1},\tilde{F}_{2}$ be the cumulative
distributions of $X_{1}+Z_{1}$ and $X_{2}+Z_{2}$. Under the assumption that
the vectors $X$ and $Z$ are independent, if we denote
$u_{1}=\tilde{F}_{1}(x_{1}),u_{1}=\tilde{F}_{2}(x_{2})$ and
$u_{3}=F_{3}(x_{3})$, then we have by Sklar’s theorem

$C_{1}(u_{1},u_{2},u_{3})=P(X_{1}+Z_{1}\leq x_{1},X_{2}+Z_{2}\leq
x_{2},X_{3}\leq x_{3}).$

Using independence of $X$ and $Z$ and computing by conditioning, we obtain

$C_{1}(u_{1},u_{2},u_{3})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}P(X_{1}\leq
x_{1}-t_{1},X_{2}\leq x_{2}-t_{2},X_{3}\leq x_{3})dG_{1}(t_{1})dG_{2}(t_{2})$

$=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(x_{1}-t_{1},x_{2}-t_{2},x_{3})dG_{1}(t_{1})dG_{2}(t_{2}).$

The last equality gives the needed fomula, with
$x_{i}=\tilde{F}_{i}^{-1}(u_{1})$ using Sklar’s theorem. The proof of Theorem
1.3.3 is similar to that of Theorem 1.3.1 and the corollaries are special
cases.

###  3.3 Proof of Theorem 1.5.1

Copulas of interest are convex combinations of copulas. The dependence
coefficients here are linear functionals of convex combinations of copulas and
are equal to $0$ for $\Pi$ and $1$ for $M$, except for
$\gamma(C_{\theta,\Pi})=2$ and $\gamma(C_{\theta,M})=3$. So, Theorem 1.5.1
holds.

###  3.4 Proof of Theorem 2.1.2 , Theorem 2.1.5 and Theorem 2.2.2

Let $(X_{k},k\in\mathbb{N})$ be a copula-based Markov chain generated by
$C(u,v)$ and $(\tilde{X}_{k},k\in\mathbb{N})$ a Markov chain generated by the
perturbation. Recall that $P^{n}(A\cap B)=P(X_{1}\in A,X_{n+1}\in B)$ and
$\tilde{P}^{n}(A\cap B)=P(\tilde{X}_{1}\in A,\tilde{X}_{n+1}\in B)$. When the
perturbation uses $\Pi(u,v)$, formula (1.6) leads to

$\tilde{P}^{n}(A\cap B)-P(A)P(B)=(1-\theta)^{n}(P^{n}(A\cap B)-P(A)P(B))$.
Therefore, $\alpha_{n}(C_{\theta,\Pi})=(1-\theta)^{n}\alpha_{n}(C).$ The rest
follows from the fact that $\alpha_{n}(C)$ is bounded.

When the perturbation uses $M(u,v)$, formula (1.5) leads to

$\displaystyle\tilde{P}^{n}(A\cap B)-P(A)P(B)=\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}(P^{n}(A\cap B)-P(A)P(B))+$

$+(1-\theta)^{n}(M(A\cap B)-P(A)P(B))$.

Let $c$ be the limit of $\alpha_{n}$ if it exists. Based on the fact that the
supremum of a sum is less than the sum of suprema and the limit of
$(1-\theta)^{n}$ is 0, we conclude that
$\displaystyle|\alpha_{n}(C_{\theta,M})-c|\leq\sum_{i=1}^{n}{n\choose
i}\theta^{n-i}(1-\theta)^{i}|\alpha_{i}(C)-c|+(1-\theta)^{n}|\alpha(M)-c|\to\lim_{n\to\infty}|\alpha_{n}(C)-c|=0$.

Therefore, when $\lim_{n\to\infty}\alpha_{n}(C)=0$, we obtain
$\lim_{n\to\infty}\alpha_{n}(C_{\theta,M})=0$.

Theorem 2.1.5 is a direct consequence of Theorem 2.1.4 and formula (1.6).
Here, one of the terms of the sum being strictly less than $1/4$ in absolute
value, the entire sum is strictly less than $1/4$. The proof benefits from the
fact that we have a Markov chain and the fold product helps derive the needed
joint distribution. As for Theorem 2.2.2, the same arguments work with the
difference that the correlation is compared to 1 and the idea still works.

###  3.5 Conclusion and remarks

When dealing with Markov chains in the copula framework, the work with mixing
coefficients benefits from the fold product of copulas and this helps get more
insights as compare to the general case. We have shown the impact of
perturbations on the $\alpha$-mixing and $\rho$-mixing structure of Markov
chains generated by copulas and will be looking to extend this work to
$\psi^{\prime}$-mixing and beyond.

## References

  * [1] B.K. Beare (2010). Copulas and Temporal Dependence. Econometrica 78, Issue 1, 395–410;
  * [2] J.R. Blum, D.L. Hanson and L.H. Koopmans (1963). On the strong law of large numbers for a class of stochastic processes. Z. Wahrscheinlichkeitstheorie und Verw. Gebiete 2, 1–11;
  * [3] R.C. Bradley (2007). Introduction to Strong Mixing Conditions. Vol. 1,2, Kendrick Press;
  * [4] R.C. Bradley (2005). Basic Properties of Strong Mixing Conditions. A Survey and Some Open Questions. Probability surveys 2, 107-144;
  * [5] R.C. Bradley (1983). On the $\psi$-mixing condition for stationary random sequences. Transactions of the American Mathematical Society, 276(1) 55–66;
  * [6] R. Cogburn (1960). Asymptotic properties of stationary sequences. Univ. Calif. Publ. Statist. 3 99–146;
  * [7] W. F. Darsow, B. Nguyen, E. T. Olsen (1992). Copulas and Markov processes. Illinois journal of mathematics 36(4) 600–642;
  * [8] V. Gaposhkin (1991). Moment bounds for integrals of $\rho$-mixing fields. Theory Probab. Appl. 36, 249–260;
  * [9] I.A. Ibragimov (1959). Some limit theorems for stochastic processes stationary in the strict sense. Dokl. Akad. Nauk SSSR 125 711–714;
  * [10] R. Ibragimov (2009). Copula-based characterizations for higher-order Markov processes. Econometric Theory 25 819–846;
  * [11] A.N. Kolmogorov and Yu.A. Rozanov (1960). On strong mixing conditions for stationary Gaussian processes. Theor. Probab. Appl.5 204–208;
  * [12] J. Komornik, M. Komornikova, J. Kalicka (2017). Dependence measures for perturbations of copulas. Fuzzy Sets and Systems 324 100–116;
  * [13] F. Durante, J.F. Sanchez, M.U. Flores (2013). Bivariate copulas generated by perturbations. Fuzzy Sets and Systems 228 137–144;
  * [14] M. Longla, F. Djongreba Ndikwa, M. Muia Nthiani (2021). Perturbations of copulas and Mixing properties. Journal of the Korean Statistical Society doi: 10.1007/s42952-021-00133-5;
  * [15] M. Longla (2015). On mixtures of copulas and mixing coefficients. Journal of Multivariate analysis 139, 259–265;
  * [16] M. Longla (2014). On dependence structure of copula-based Markov chains. ESAIM: Probability and Statistics 18, 570–583;
  * [17] M. Longla (2013). Remarks on the speed of convergence of mixing coefficients and applications. Statistics and probability letters 83(10); 2439–2445;
  * [18] M Longla, M Peligrad (2012). Some aspects of modeling dependence in copula-based Markov chains Journal of Multivariate Analysis 111, 234–240;
  * [19] R. Mesiar, M. Komornikova, J. Komornik (2015). Perturbation of bivariate copula. Fuzzy Sets and Systems 268 127–140;
  * [20] C. Miller (1994). Three theorems on $\rho^{*}$-mixing random fields. J. Theor. Probab. 7, 867–882;
  * [21] C. Miller (1995). A CLT for the periodograms of a $\rho^{*}$-mixing random field. Stochastic Process. Appl 60, 313–330;
  * [22] R.B. Nelsen. An Introduction to Copulas, second edition, Springer Series in Statistics, Springer-Verlag, New York;
  * [23] M. Rosenblatt (1956). A central limit theorem and a strong mixing condition. Proc. Natl. Acad. Sci. USA 42 43–47;
  * [24] A. Sklar (1959). Fonctions de répartition à $n$ dimensions et leurs marges.Publ. Inst. Statist. Univ. Paris, 8: 229–231;


[a]Simon De Kockere

# Simulation of the propagation of CR air shower cores in ice

Krijn de Vries  Nick van Eijndhoven

###### Abstract

Currently new radio detection techniques are being explored to detect
astrophysical neutrinos beyond the PeV scale interacting in polar ice. Due to
the long attenuation length of radio waves in a medium, it can be expected
that such instruments will also be sensitive to the radio emission of cosmic
ray air showers. Furthermore, cosmic ray air showers hitting a high-altitude
layer of ice will initiate an in-ice particle cascade, also leading to radio
emission. We present the first results of detailed simulations of the in-ice
continuation of these cosmic-ray-induced particle cascades, using a
combination of the CORSIKA Monte Carlo code and the Geant4 simulation toolkit.
We give an overview of the general features of such particle cascades and
present a parameterization in terms of Xmax of the longitudinal and lateral
particle distributions. We discuss the feasibility of observing the in-ice
particle cascades, both through the detection of the Askaryan radio emission
as well as by using the RADAR reflection technique. Based on these results we
find that the expected signals from the continuation of in-ice cosmic-ray
induced particle cascades will be very similar to neutrino signals. This means
a thorough understanding of these events is necessary in the search for
neutrino candidates, while it also promises an interesting in-situ natural
calibration source.

##  1 Introduction

A high-energy neutrino interacting in a dense medium like ice will initiate a
particle cascade with an excess of negative charges. This negative charge
excess will lead to coherent radio emission, as predicted by Askaryan [1]. Due
to the long attenuation length of radio in ice, this mechanism provides a
valuable way of studying astrophysical neutrinos of the highest energies,
efficiently extending the energy range of Cherenkov light neutrino
observatories like the IceCube Neutrino Observatory. This technique is
currently being explored at the Polar regions by several experiments, e.g. [2,
3, 4, 5].

However, also cosmic-ray induced air showers create coherent radio emission
reaching the surface, and can thus be expected to form an important background
for radio neutrino observatories [6, 7]. Furthermore, at high altitudes of 2 -
3 km, corresponding to the altitudes of [2] and [5], cosmic-ray induced air
showers will typically still have a very energy dense core. The propagation of
these energy dense shower cores through the ice will lead to Askaryan radio
emission, forming an additional background component.

Only a few studies on the propagation of these shower cores through ice and
its corresponding radio emission exist [8, 9, 10, 11], and to fully understand
this background component a more detailed analysis is required. Moreover, not
only will these shower cores lead to background signals, if well understood
they can also serve as a free _in-situ_ calibration source. Here we present
the first results of a detailed simulation of the propagation of cosmic ray
air shower cores in ice and its corresponding radio emission, using a
combination of the CORSIKA Monte Carlo code and the Geant4 simulation toolkit.

##  2 Cosmic ray air showers

The development of a cosmic ray air shower is a statistical process, resulting
in the fluctuation of properties from shower to shower. For illustration, we
will discuss the properties of a single proton induced air shower generated
with the CORSIKA 7.7100 Monte Carlo code [12]. More precisely, we will discuss
its properties around an altitude of 2.4 km, a typical value for the surface
of a high-altitude Polar ice sheet. The proton has a primary energy
$E_{p}=10^{17}$ eV and is generated at a zenith angle $\theta=0^{\circ}$. We
use the QGSJETII-04 high energy hadronic interaction model, the GHEISHA 2002d
low energy hadronic interaction model and a MSIS-E-90 atmospheric model for
South Pole on December 31, 1997. We apply thinning for electromagnetic
particles falling below 10 GeV, using a maximum thinning weight of 10. For
hadrons (without $\pi^{0}$’s) and muons a kinetic energy cut-off of 0.3 GeV
was used. For electrons, photons and $\pi^{0}$’s a kinetic energy cut-off of
0.003 GeV was used.

Figure 1 shows the longitudinal profile of the particle cascade in air. The
cascade reaches a maximum number of $e^{+}+e^{-}$ at a depth of $X_{max}$ =
680 g/cm2, which is a good average value, as can be seen in e.g. Figure 2 of
[13]. From Figure 1 it can be seen that the air shower reaches the altitude of
2.4 km right after this shower maximum. At this point in the shower, gamma
radiation, electrons and positrons clearly outnumber the muonic and hadronic
part. In terms of energy however, the difference is less significant, and the
muons and hadrons cannot simply be ignored.

Figure 1:  The number of particles (_left)_ and the distribution of the energy
(_right_) of a proton induced air shower, simulated with the CORSIKA Monte
Carlo code. The proton has a primary energy $E_{p}=10^{17}$ eV and is
generated at a zenith angle $\theta=0^{\circ}$. The shower reaches a maximum
number of $e^{+}+e^{-}$ at a depth $X_{max}=680$ g/cm2. The vertical black
line indicates the depth value corresponding to an altitude of 2.4 km.

In Figure 3 we show the average kinetic energy per particle at an altitude of
2.4 km, calculated over radial intervals of $\Delta r=0.1$ m. For energies of
80 MeV and above, electromagnetic particles will on average be able to
interact, leading to pair-production, continuing the development of the
particle cascade. Below this energy ionization losses start to dominate, and
the particle cascade dies out. So as indicated in the figure, although the
shower reaches this altitude right after shower maximum, we still expect the
shower core to develop further, while in total the particle cascade will start
to die out. Figure 3 shows the total energy within a given radius at the same
altitude of 2.4 km. Here we clearly see the shower has a very energy dense
core. The first 100 cm of the particle cascade already contains about 15% of
the initial energy of the primary proton.

Figure 2:  The average kinetic energy per particle at an altitude of 2.4 km.
For the electromagnetic part ($\gamma$, $e^{+}$, $e^{-}$) the average was
calculated over radial intervals of $\Delta r=0.1$ m. For the muonic and
hadronic part radial intervals of $\Delta r=0.5$ m were used. The black
horizontal line indicates the critical energy of 80 MeV.

Figure 3:  The total energy within a given radius at an altitude of 2.4 km.

##  3 The in-ice core

For the simulation of the propagation of cosmic ray air shower cores in ice,
we combine the CORSIKA Monte Carlo code with the Geant4 10.5 simulation
toolkit [14].

Using CORSIKA we calculate the momentum, position and arrival time of all
secondary particles reaching an altitude of 2.4 km. This information is passed
on to a Geant4 module, which then propagates each secondary particle through
ice.

The Geant4 module simulates a volume consisting of multiple horizontal layers
of pure ice, each with a constant density and a thickness of 1 cm. The density
of the layers follows the density profile measured at the Antarctic Taylor
Dome ice cap, given by $\rho(z)=0.460+0.468\cdot(1-e^{-0.02z})$, with $\rho$
the density in g/cm3 and $z$ the depth in m. We include the
G4EmStandardPhysics constructor, providing us with the standard Geant4
electromagnetic physics, the G4DecayPhysics constructor for the decay physics
of long-lived hadrons and leptons and the G4RadioactiveDecayPhysics
constructor for radioactive decay physics. We use the default cut-off length
of 1 mm for gammas, electrons, positrons and protons, meaning that particles
of these types created during the Geant4 simulation are only included in the
simulation when they can travel distances larger than this cut-off length.
This replaces a hard kinetic energy cut-off.

Figure 4 shows the energy deposited in the ice by the reference shower
discussed in Section 2. We see again that the energy is highly concentrated
close to the core of the particle cascade, within a radius of the order of 10
cm, resembling a neutrino induced particle cascade. Furthermore, as discussed
above, we can also see that the shower core is still developing when entering
the ice.

Figure 4:  The deposited energy density in a vertical 1 cm wide slice of the
ice going though the center of the particle cascade (_left_) and the deposited
energy density in the ice in function of radius and depth (_right_). The black
line marks the area for which $\omega_{p}>100$ MHz (see Section 4).

In Figure 6 we show the number of particles in function of depth of the
reference shower when propagated through ice, and compare with the case where
it continues developing in air instead until reaching sea level. As can be
seen, the transition from air into ice does not influence the development of
the electromagnetic part of the cascade, and we simply see a continuation of
the in-air profile to larger depth values. This means that standard air shower
parameterizations such as the Gaisser-Hillas parameterization [15] can be used
to describe the electromagnetic part of the cascade in the ice.

Figure 6 shows a snapshot of the particle cascade at a time $t=20.5$ ns after
the moment of impact of the shower core front on the ice, corresponding with a
cascade front depth of 300 g/cm2 with respect to the ice surface. As can be
seen, the particles are concentrated in a thin disk, with a thickness of the
order of 1 cm. As such, the lateral dimension will be the relevant dimension
when studying the corresponding coherent radio emission.

Figure 5:  The number of particles in function of depth of the simulated air
shower. The dashed lines show the case where the particle cascades propagates
through air until reaching sea level. The solid lines show the cascade
development when propagated through ice at an altitude of 2.4 km, indicated by
the vertical black line. Here we exceptionally use the hard energy cut-offs of
CORSIKA in the Geant4 simulation as well.

Figure 6:  The number of particles per cm3 in function of radius and depth at
a time $t=20.5$ ns after the moment of impact of the shower core front on the
ice.

We describe the lateral particle distribution of the in-ice cascade by
constructing the $w_{1}(r)$ distribution. By definition, $w_{1}(r)dr$
represents the fractional number of charged particles within the interval
$[r,r+dr[$ at a given time, normalized such that the integral from $r=0$ to
$r=1$ m gives 1. Here, $r$ represents the radius in the shower axis frame. As
the cascade develops through time, this distribution will also be time
dependent. Figure 8 shows the $w_{1}(r)$ distribution for the reference shower
at different points in time, indicated by the depth of the cascade front with
respect to the ice surface. They show a peak around $r=3$ cm, followed by a
long tail.

To include a wider range of primary energies and zenith angles, we create 10
different shower sets, each covering primary energies of $10^{16}$ eV up to
$10^{18}$ eV in steps of half a decade111Thinning was applied in CORSIKA for
showers with primary energies $E_{p}$ $\geq$ $10^{17}$ eV on electromagnetic
particles falling below $10^{-7}$ $E_{p}$ with a thinning weight smaller than
$10^{-7}$ $E_{p}[\text{GeV}]$., and zenith angles of $0^{\circ}$, $15^{\circ}$
and $30^{\circ}$. This results in 150 showers in total, including the
reference shower from Section 2. Within each set, the random seeds used to run
the CORSIKA simulations are identical. In order to derive a simple
parameterization of the $w_{1}(r)$ distributions, we group the showers
together based on their $X_{max}$ value, ignoring primary energy and zenith
angle. For every group we then construct the $w_{1}(r)$ distributions for a
given cascade front depth, and calculate the mean $w_{1}(r)$ distribution. The
result is shown in Figure 8, for a cascade front depth of 450 g/cm2 with
respect to the ice surface. As can be seen, the peak value of the $w_{1}(r)$
distribution correlates with the $X_{max}$ value of the particle shower. A
higher $X_{max}$ value implies a denser shower core, resulting in a sharper
peak in the $w_{1}(r)$ distribution.

Figure 7:  The $w_{1}(r)$ distributions of the reference shower at different
points in time, indicated by the depth of the cascade front with respect to
the ice surface.

Figure 8:  The average $w_{1}(r)$ distributions at a cascade front depth of
450 g/cm2 with respect to the ice surface of the simulation set covering 150
showers. In the legend the lower and upper limits of the different $X_{max}$
groups are given.

We find that the average $w_{1}(r)$ distributions can be well described by the
analytical expression

$W(r)=a\sqrt{r}e^{-\left(r/b\right)^{c}},$ (1)

with the parameters $a$, $b$ and $c$ depending on the $X_{max}$ value of the
shower. Both $\log(a)$, $b$ and $c$ follow a linear relationship with
$X_{max}$, as shown in Figure 9. This means that for a shower with known
$X_{max}$ value, the parameters $a$, $b$ and $c$ for a given cascade front
depth can be constructed, which can be combined with Equation 1 to get the
corresponding $w_{1}(r)$ distribution. The $w_{1}(r)$ distribution at a
cascade front depth of 450 g/cm2 with respect to the ice surface for the
reference shower from Section 2 obtained this way is shown in Figure 10, where
we also compare with the corresponding distribution calculated directly from
the Geant4 simulation.

Figure 9:  The fit parameters of Equation 1 for the average $X_{max}$ values
of each group, at a cascade front depth of 450 g/cm2 with respect to the ice
surface. Figure 10:  The $w_{1}(r)$ distribution at a cascade front depth of
450 g/cm2 with respect to the ice surface for the reference shower from
Section 2. The dashed red line shows the distribution calculated directly from
the Geant4 simulation. The solid black line shows the reconstruction of the
distribution using Equation 1 and the linear fits from Figure 9, evaluated at
$X_{max}=680$ g/cm2.

##  4 Applications

To find a first estimate of the Askaryan radio emission associated with the
propagation of cosmic ray air shower cores in ice, we implement the end-point
formalism [16] in the Geant4 module, using the code from the work presented in
[17] as an example. This formalism is also implemented in the CoREAS extension
of CORSIKA, which calculates the radio emission of air showers, and is shown
to agree well with experimental results [18]. However, this approach assumes a
constant index of refraction $n$ of the medium, which works well for air, but
might be an oversimplification for the top layers of natural ice. We use
$n=1.52$.

Figure 11 shows the electrical field in function of time as seen by a point-
like antenna, for the reference shower from Section 2. The antenna is
positioned 150 m deep in the ice ($z=-150$ m), at a horizontal distance of 160
m away from the point of impact of the shower core on the ice surface
($x=-160$ m), which is close to the expected Cherenkov angle. The signal is
bipolar and radial polarized, as expected for Askaryan radiation, and reaches
a magnitude well above typical detection thresholds of 10-100 $\mu$V/m. Here
it should be noted that for a quantitative statement the induced field should
be convolved with the antenna, which is foreseen in future work.

Figure 11: The electrical field in function of time as seen by a point-like
antenna at $x=-160$ m and $z=-150$ m with respect to the point of impact on
the ice surface, for the reference shower from Section 2. It was calculated
using the endpoint formalism.

As a second application we discuss the feasibility of detecting cosmic ray air
shower cores propagating in ice with RADAR reflection techniques [19]. The
shower core will create a dense plasma in the ice, which should be able to
reflect incoming radio waves. We estimate the plasma frequency $\omega_{p}$ of
the plasma created in the ice, indicating the transition between fully
coherent and incoherent scattering for a collisionless plasma. Although the
in-ice ionization plasma suffers from collisions, as a rule of thumb, signals
with a frequency $\omega<\omega_{p}$ will be reflected by the plasma.

The plasma frequency is given by
$\omega_{p}=8980\sqrt{n_{q}[\text{cm}^{-3}]}\text{ Hz}$, with $n_{q}$ the free
charge density of the plasma. We estimate $n_{q}$ based on the deposited
energy density $\rho_{E}$ in the ice, using $n_{q}=\rho_{E}/(50\text{ eV})$.
As indicated in Figure 4, the plasma reaches plasma frequencies of 100 MHz and
higher, which agrees well with realistic RADAR designs [20].

## References

  * [1] G.A. Askaryan, Sov. Phys. JETP 14/21, 441/658 (1962/1965)
  * [2] ARA collaboration, Astropart. Phys. 35, 457 (2012)
  * [3] ARIANNA collaboration, IEEE Trans. Nucl. Sci. 62, 2202 (2015)
  * [4] ANITA collaboration, Astropart. Phys. 32, 10 (2009)
  * [5] RNO-G collaboration, JINST 16, P03025 (2021)
  * [6] Frank G. Schröder, Progress in Particle and Nuclear Physics 93, 1-68 (2017)
  * [7] Tim Huege, Dave Besson, Prog. Theor. Exp. Phys. 12, 12A106 (2017)
  * [8] S. Razzaque et al., Phys. Rev. D 65, 103002 (2002)
  * [9] D. Seckel et al., Proceedings of the 30th International Cosmic Ray Conference 5, 1029 (2008)
  * [10] Amir Javaid, Monte Carlo simulation for radio detection of ultra-high energy air shower cores by ANITA-II, University of Delaware (2012)
  * [11] Krijn D. de Vries et al., Astropart.Phys. 74, 96-104 (2016)
  * [12] D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz and T. Thouw, FZKA, 6019, 90 (1998)
  * [13] S. Buitink et al., Nature 531, 70 (2016)
  * [14] S.Agostinelli et al., Nucl. Instrum. Methods Phys. Res. A, 506, 250-303 (2003)
  * [15] T.K. Gaisser and A.M. Hillas, Proc. 15th ICRC, Plovdiv, Bulgaria, 8, 353 (1977)
  * [16] Clancy W. James et al., Physical Review E, 84, 5 (2011)
  * [17] Anne Zilles, Emission of Radio Waves in Particle Showers: Validation of microscopic simulations with the SLAC T-510 experiment and their potential in the future Square Kilometer Array, Springer (2017)
  * [18] Anna Nelles et al., Astroparticle Physics, 60, 13-24 (2014)
  * [19] Krijn D. de Vries, Kael Hanson, Thomas Meures, Astropart.Phys. 60 ,25-31 (2015)
  * [20] S. Prohira et al., arXiv:2104.00459 (2021)