# Evolution of reconnection flux during eruption of magnetic flux ropes

Samriddhi Sankar Maity Indian Institute of Astrophysics, Bangalore-560034,
India  Joint Astronomy Programme and Department of Physics, Indian Institute
of Science, Bangalore, India  Piyali Chatterjee Indian Institute of
Astrophysics, Bangalore, India  Ranadeep Sarkar University of Helsinki,
Finland  Ijas S. Mytheen Amrita School of Physical Science, Amrita Vishwa
Vidyapeetham, Kollam, India

###### Abstract

Coronal mass ejections (CMEs) are powerful drivers of space weather, with
magnetic flux ropes (MFRs) widely regarded as their primary precursors.
However, the variation in reconnection flux during the evolution of MFR during
CME eruptions remains poorly understood. In this paper, we develop a 3D
magneto-hydrodynamic model using which we explore the temporal evolution of
reconnection flux during the MFR evolution using both numerical simulations
and observational data. Our initial coronal configuration features an
isothermal atmosphere and a potential arcade magnetic field beneath which an
MFR emerges at the lower boundary. As the MFR rises, we observe significant
stretching and compression of the overlying magnetic field beneath it.
Magnetic reconnection begins with the gradual formation of a current sheet,
eventually culminating with the impulsive expulsion of the flux rope. We
analyze the temporal evolution of reconnection fluxes during two successive
MFR eruptions while continuously emerging the twisted flux rope through the
lower boundary. We also conduct a similar analysis using observational data
from the Helioseismic and Magnetic Imager (HMI) and the Atmospheric Imaging
Assembly (AIA) for an eruptive event. Comparing our MHD simulation with
observational data, we find that reconnection flux play a crucial role in
determination of CME speeds. From the onset to the eruption, the reconnection
flux shows a strong linear correlation with the velocity. This simulation of a
solar eruption provides important insights into the complex dynamics of CME
initiation and progression.

##  1 Introduction

Coronal Mass Ejections are composed of clouds of magnetized plasma that are
expelled from the Sun into the heliosphere due to sudden release of free
magnetic energy stored in the twisted coronal magnetic field (Chen, 2017).
They are of interest due to scientific and technological reasons since CMEs
can drive interplanetary shocks that energize solar particles and cause
significant space weather effects in the geospace. In-situ data obtained by
satellites passing through the interplanetary CMEs (ICMEs) have established
that typical ICMEs have the structure of highly twisted MFR (Wang et al.,
2016; Hu, 2017). Many CME models therefore incorporate a magnetic flux rope–
consisting of helical field lines twisting about a central axis in the corona
as the basic underlying magnetic field structure for CME precursors (Titov &
Démoulin, 1999; Gibson & Fan, 2006; Duan et al., 2019; Liu, 2020; Chen, 2017).
Therefore, understanding their evolution in early phase is crucial for CME
studies.

Although is it widely accepted that MFRs constitute the core structure of the
CMEs, it remains unknown wheather MFRs exist in the solar corona before CME
initiation or forms during the eruption (Chen, 2011; Patsourakos et al.,
2020). Some opine that the MFR may exist before eruptions (Cheng et al.,
2011), although there is no agreement on how and where an MFR might form. An
alternate hypothesis is that the MFRs can emerge from beneath the photosphere
(Fan, 2001; Martínez-Sykora et al., 2008; Magara, 2004; Archontis et al.,
2009). Yet another proposition is that the MFRs can be built directly in the
corona via shearing of magnetic footpoints and reconnection prior to the
eruption (van Ballegooijen & Martens, 1989; Amari et al., 2003; Aulanier et
al., 2010; Chatterjee et al., 2016). Most eruptions, particularly those
originating from solar active regions (ARs), occur along magnetic polarity
inversion lines (PILs) within strong field regions. Furthermore, some flare-
productive ARs exhibit relatively short time interval between successive
eruptions, while displaying very similar structure in CME morphology. This
kind of event is referred to as homologous eruptions (Zhang & Wang, 2002).
Observations indicate that the evolution of the source regions of these
homologous events often display specific evolutionary pattern, including (a)
continuous shearing motion (Li et al., 2010; Romano et al., 2015, 2018; Sarkar
et al., 2019), (b) rotation of sunspots (Régnier & Canfield, 2006; Zhang et
al., 2008), and (c) emergence of flux (Ranns et al., 2000; Nitta & Hudson,
2001; Sterling & Moore, 2001; Dun et al., 2007; Xu et al., 2017).

To explain solar eruptions many theories have been suggested (Forbes et al.,
2006; Shibata & Magara, 2011; MacTaggart & Hood, 2009; Schmieder et al., 2013;
Aulanier, 2014; Janvier et al., 2015), but only two main models have been
successful to explain homologous eruptions using magnetohydrodynamics (MHD)
simulations. These are flux emergence model (Chatterjee & Fan, 2013) and the
breakout model (DeVore & Antiochos, 2008). The flux emergence model suggests
that new magnetic flux rises from below the solar surface, bringing with it
twisted magnetic field lines created by turbulent convection that trigger
eruptions (Fan, 2010). Chatterjee & Fan (2013) demonstrated repeated CMEs
caused by a highly twisted flux tube emerging into the solar corona. In their
simulation, a coronal flux rope partially erupted and reformed multiple times.
However, the time gap between two consecutive eruptions was insufficient, the
post-eruption field did not have enough time to return to a relax state before
the onset of the subsequent eruption. In contrast, the breakout model features
a multipolar magnetic configuration with a null point located above a central
flux system which is sheared by photospheric motion. As the shearing
increases, the magnetic reconnection begins at the null point. This
reconnection process acts as a trigger for an eruption. When the eruption
occurs, it removes the overlying magnetic field above the core that was
previously restraining the core. After the eruption, the original magnetic
structure is restored (Antiochos et al., 1999). Continuous shearing motions
can result in repetition of such mechanisms and give rise to multiple
eruptions, which are otherwise tend to be confined (DeVore & Antiochos, 2008).

Figure 1: Location of STEREO spacecraft and Earth on 2011 August 4.

Although magnetic reconnection can only be indirectly observed, it is a
critical process in the solar corona that results in the formation of closed
magnetic loops which energizes the plasma and particles leading to sudden
increase in flare radiation. The correlation suggests a connection between the
CME evolution kinematics and the reconnection (Zhang et al., 2001; Temmer et
al., 2008; Bein et al., 2012; Patsourakos et al., 2013). Magnetic reconnection
is a ubiquitous process during coronal mass ejection (CME) events, but its
precise role in the eruption sequence remains a subject of debate. It’s not
definitively established whether reconnection acts as the trigger for the
eruption or if it’s a result of the eruption process. Many have also utilized
observations of flare emission patterns in the lower solar atmosphere as a
means to indirectly study and infer the characteristics of reconnection
processes occurring higher up in the corona (Forbes & Priest, 1984; Kopp &
Poletto, 1986). Observing the early evolution and acceleration of coronal mass
ejections (CMEs) in the lower corona presents significant challenges. As a
result, many studies opt for an alternative approach. They focus on comparing
two key parameters: the CME velocity, measured at several solar radii after
the CME has undergone its peak acceleration, and the total magnetic flux that
was reconnected during the associated flare event. Note that most of the
studies have focused on at a instant after the peak phase and it is unclear
weather the evolution flux is at all related to the velocity. Numerical MHD
simulations have proven to be powerful tools for reproducing the time-
dependent, nonlinear evolution of 3D magnetic configurations and investigating
the temporal changes in flux during the early phases of magnetic flux rope
(MFR) eruptions. A key property characterizing solar flares is the amount of
magnetic flux passing through the reconnection sheet beneath the MFR, commonly
referred to as the reconnection flux (Gopalswamy et al., 2017). The magnetic
flux reconnected in the solar corona cannot be directly observed or measured.
However, one can directly quantify it through its relationship with the
magnetic flux encompasses by the flare ribbon or the foot points of the MFR
(Kazachenko et al., 2017; Qiu et al., 2004).

Figure 2: Different phases of the filament eruption on 2011 August 04 as
observed by STEREO-A. (a) Image of the Sun in STEREO-A EUVI 304 Å (̇b) COR1
image superimposed with STEREO-A EUVI 304 Å (̇c) Superimposed images of
STEREO-A EUVI 304 Å , COR1 and COR2. The yellow arrows in panels (a) and (b)
indicate the filament leading edge. The images are plotted using JHelioviewer
(https://www.jhelioviewer.org/).

In this paper, we conduct a 3-dimensional compressible MHD simulation to study
the time evolution of the reconnection flux and correlate it with the velocity
of the MFR during its early evolution. We further compare our analysis with an
observed event. The observational data and methods is mentioned in Section 2,
the numerical setup is presented in Section 3. The results are mentioned in
Section 4 and we conclude with a summary in Section 5.

Figure 3: Left panel: Flare associated brightening observed in AIA 1600 Å
channel during the flare occurred in AR 11261. Right panel: The associated HMI
$B_{r}$ magnetic field (gray color scale with saturation value $\pm$500 Gauss)
and the temporal evolution of the flare ribbons from 03:41 UT to 04:17 UT on
2011 August 04 .

##  2 Observation Data

We study the temporal evolution of re-connection flux and the associated CME
speed during the M-class flare (SOL2011-SOL2011-08-04T03:41) which occurred in
AR 11261 during 2011 August 4. During this period, STEREO-A and STEREO-B were
positioned at longitudinal separations of 101∘ and 92∘ (see Figure 1),
respectively, relative to the Sun-Earth line. This positioning allowed the two
spacecraft to provide near-limb views of Earth-directed CMEs. On the other
hand the location of AR 11261 was close to the solar disk center (between
30∘to 35∘ west) relative to Earth, offering reliable data on the solar
magnetic field of the AR. This provides a unique opportunity to study the co-
temporal evolution of re-connection flux, estimated from on-disk SDO data, and
the speed of the associated eruption, starting from the initiation height at
lower corona as observed from near-limb observations by STEREO.

Figure 4: (a) Time profile of signed cumulative reconnection flux integrated
over the positive (red) and negative (blue) magnetic polarities underlying the
flare ribbons. (b) Time profile of unsigned instantaneous reconnection flux
integrated over the positive (red) and negative (blue) magnetic polarities
underlying the flare ribbons. (c) Height time plot of the filament leading
edge measured from STEREO EUVI, COR1 and COR2 for the same event as shown in
Figure 2. (d) Velocity profile of the filament leading edge derived from Panel
(c).

We analyze the kinematics of the eruption associated with the M-class flare on
2011 August 04, by tracking the filament leading-edge height as observed by
STEREO-A. We use two instruments from the Sun Earth Connection Coronal and
Heliospheric Investigation (SECCHI; Howard et al. (2008)) onboard STEREO-A to
track the early evolution of the Earth-directed filament eruption, which was
observed as prominence eruption from STEREO-A (see Figure 2). The instruments
include the Extreme Ultraviolet Imager (EUVI; $<1.7R_{\odot}$, Wuelser et al.
(2004)) and the inner coronagraph COR1 ($1.3-4R_{\odot}$, Thompson et al.
(2003)).

We follow the method by (Kazachenko, 2023), to estimate the temporal profile
of reconnecting flux. The observations of the flare ribbons are obtained from
the Atmospheric Imaging Assembly (AIA; Lemen et al. (2012)) onboard SDO in
1600 Å channel. The associated full-disk magnetograms are obtained from the
Helioseismic and Magnetic Imager (HMI; Schou et al. (2012)). During the
recurrent eruptive M-class flare in AR 11261 on 2011 August 4, we identify the
area swept out by the flare-ribbons observed in AIA 1600 Å images as shown in
Figure 3. The half of the unsigned magnetic flux underlying the area
associated with the brightening observed in AIA 1600 Å are recorded as the
temporal profile of the re-connection flux. The kinematics of the filament
eruption and the temporal variation in associated reconnection flux are shown
in Figure 4.

Figure 5: The radiative cooling function in used in this work is a modified
version of Cook et al. (1989)

##  3 Numerical Model

We numerically solve the complete magneto-hydrodynamic equations in three-
dimensional spherical coordinates to investigate the early evolution of
magnetic flux ropes in the solar corona. For this, we employ the Pencil Code
(Pencil Code Collaboration et al., 2021), an open-source, highly modular, and
MPI-parallelized code designed for compressible MHD flows 111http://pencil-
code.nordita.org. Utilizing a sixth-order finite difference scheme and a
third-order Runge-Kutta time-stepping method—among various other options
provided by the Pencil Code—we solve the following set of compressible MHD
equations.

$\frac{\mathrm{D}\ln\rho}{\mathrm{D}t}=-\mathbf{\nabla}\cdot\mathbf{U}$ (1)

$\frac{\mathrm{D}\mathbf{U}}{\mathrm{D}t}=-\frac{\mathbf{\nabla}p}{\rho}+\frac{GM_{\odot}}{r^{2}}\hat{\mathbf{r}}+\frac{\mathbf{J}\times\mathbf{B}}{\rho}+\mathbf{F}_{\mathrm{corr}}+\rho^{-1}\mathbf{F}_{\mathrm{visc}}$
(2)

$\frac{\partial\mathbf{A}}{\partial
t}=\mathbf{U}\times\mathbf{B}-\eta\mu_{0}\mathbf{J}$ (3) $\begin{split}\rho
c_{v}{T}\frac{\mathrm{D}\ln{T}}{\mathrm{D}t}=-\rho
c_{v}{T}(\gamma-1)\nabla\cdot\mathbf{U}+\eta\mu_{0}\mathbf{J}^{2}+2\rho\nu{\sf{S}}_{ij}^{2}\\\
+\nabla\cdot\mathbf{q}_{\mathrm{cond}}-\rho^{2}\Lambda(T)+\mathcal{H}+{\mathbf{\rho}\zeta_{\mathrm{shock}}\left(\nabla\cdot\mathbf{U}\right)^{2}}\end{split}$
(4)

In the above, $\mathbf{U}$ is the velocity field; $\mathbf{B}$ is the magnetic
field; $\mathbf{A}$ is the magnetic potential; $\mathbf{J}$ is the current
density; $\rho$, $p$ and $T$ are, respectively, the plasma density, pressure
and temperature of the system. The viscous force is modelled as

$\rho^{-1}\mathrm{F_{visc}}=\nabla\cdot(2\rho\nu\mathbf{S})+\nabla\left(\rho\zeta_{shock}\nabla\cdot\mathbf{U}\right)$
(5)

where $\nu$ is kinematic viscosity and $\mathrm{S}$ is the traceless rate-of-
strain tensor. Additionally, we incorporate a semi-relativistic Boris
correction into the classical MHD momentum equation to reduce numerical
diffusion and avoid the need for excessively small time steps in MHD
simulations, achieved by using an artificially reduced speed of light
(Chatterjee, 2020).

We include radiative cooling in our model, represented in equation (4), where
$\Lambda(T)$ denotes the radiative loss function as described by (Cook et al.,
1989). The cooling function is plotted in Figure (5). Furthermore, we use a
simplified coronal heating function that varies only with height above the
solar surface. This function follows an exponential decrease as the height
increases.

$\mathcal{H}=\frac{F}{L_{H}}\frac{R_{\odot}^{2}}{r^{2}}\mathrm{exp}\left\\{-(r-R_{\odot})/L_{H}\right\\}$
(6)

where the input energy flux density is $F=9.74\times
10^{5}\,\mathrm{erg\,cm^{-2}\,s^{-1}}$ and the decay length is
$L_{H}=1.948\times 10^{10}\,\mathrm{cm}$.

We have also included field-aligned Spitzer thermal conduction in our model.
Here we have used the hyperbolic diffusion equation approach instead of
directly including Spitzer conduction in equation (4), similar to Chatterjee
(2020).

Let $\mathbf{q}_{\mathrm{cond}}$ represents the solution of the non-Fickian
transport equation and $\mathbf{q}_{\mathrm{sp}}$ denotes the conduction flux
according to the Spitzer model. We solve the following equation for the heat
flux $\mathbf{q}_{\mathrm{cond}}$.

$\frac{\partial\mathbf{q}_{\mathrm{cond}}}{\partial
t}=-\frac{\mathbf{q}_{\mathrm{cond}}-\mathbf{q}_{\mathrm{sp}}}{\tau_{sp}}+\beta(dr\cdot\nabla)^{6}\mathbf{q}_{\mathrm{cond}}$
(7)

where
$\mathbf{q}_{\mathrm{sp}}=K_{sp}T^{5/2}\hat{\mathbf{b}}(\hat{\mathbf{b}}\cdot\nabla
T)$, $\hat{\mathbf{b}}$ denotes the unit vector along the field direction and
$K_{sp}=10^{-6}$ $\mathrm{erg\,K^{-7/2}\,cm^{-1}\,s^{-1}}$. $\tau_{sp}$
denotes a finite timescale for $\mathbf{q}_{\mathrm{cond}}$ to converge to the
Spitzer heat flux and is equal to 0.1 s. Our time step varies between 0.1 -
0.3 ms.

The computational domain consists of a spherical wedge with an extent
$R_{\odot}<r<6R_{\odot}$, $5\pi/12<\theta<7\pi/12$ and
$-\pi/9.6<\phi<\pi/9.6$. It is resolved by a grid of $512\times 288\times
160$, which is non uniform in $r$ but uniform in $\theta$ and $\phi$. The grid
spacing in $r$ is $dr=0.002R_{\odot}$ at the lower boundary which gradually
increases in a logarithmic manner reaching $dr=0.003R_{\odot}$ at the upper
boundary.

Figure 6: The twisted torus right before emerging through the bottom boundary
at $r=R_{\odot}$ (in gray). The number of turns of a single field line (black)
around the circular axis of the torus is shown. Figure 7: The 3D evolution of
the magnetic field of a the twisted flux rope as it emerges into the solar
corona at the specified times (in hours). Red field lines represents field
line anchored in the ambient arcade, while blue, green, and cyan field lines
indicate field lines originating from the emerging flux region. An animation
of this evolution is available in the online paper.

We consider a high-temperature coronal plasma as an ideal gas with an
adiabatic constant of $\gamma=1.66$ . Before the simulation starts, the domain
is considered to be in hydrostatic equilibrium at a uniform temperature of
$T_{0}=10^{6}$ K with a density stratification given by

$\rho=\rho_{0}\mathrm{exp}\left\\{-\frac{R_{\odot}}{H_{p}}(1-\frac{R_{\odot}}{r})\right\\}$
(8)

where, $\rho_{0}=1.0\times 10^{-15}$ g cm-3 is the initial density and
$H_{p0}=c_{s0}/g_{r}$ denotes the initial pressure scale height. The initial
atmosphere consists of a pre-existing potential field with arcade like
geometry as given in Fan (2012) (see their Equation 12 and 13).

At the lower boundary, we impose an electromotive force given by

$\mathbf{E}\mid_{r=R_{\odot}}=-\frac{1}{c}\mathrm{\mathbf{v}_{0}}\times\mathbf{B}_{\mathrm{torus}}(R_{\odot},\theta,\phi,t)$

which bodily transports the twisted torus radially into the domain. The major
and minor radius of the torus are $0.25R_{\odot}$ and $0.042R_{\odot}$
respectively and with a field line twist rate of 0.068 rad Mm-1. The twist
rate indicates the field line winds around the tube axis per unit length along
the axis. The flux tube is initially located at $r_{0}=0.705{R_{\odot}}$, thus
the torus is initially completely beneath the surface and then ascends bodily
upward at constant speed toward the lower boundary, until it reaches a height
$r_{\mathrm{stop}}$ when the emergence is stopped and a part of the torus is
still inside of the computational domain. The velocity field at the lower
boundary is uniformly set in the region where the emerging torus intersects
tit and its zero elsewhere. The imposed emergence speed, $v_{0}$ = 2
$\mathrm{km\,s^{-1}}$, which is much smaller than the Alfven speed $v_{A}$ =
1.69 $\mathrm{Mm\,s^{-1}}$, ensuring that the emerging flux rope evolves
quasi-statically during the flux emergence phase at the lower boundary. When
the emergence is stopped, the velocity at the lower boundary is set to zero,
with no inflows or outflows and footpoints are rigidly anchored. For the
$\theta$ boundaries, a non-penetrating stress-free condition is applied for
the velocity field while the magnetic field is treated with perfectly
electrically conducting walls. The $\phi$ boundaries are periodic. For the top
boundary, a simple outward extrapolating boundary condition is used, allowing
plasma and magnetic field to flow through.

In this simulation, we drive the emergence of the torus until
$r_{\mathrm{stop}}=0.85R_{\odot}$ to study the reconnection flux during the
evolution of the twisted flux rope in the corona which produces homologous
CMEs.

##  4 Results

Our simulation produces two homologous CME eruptions; each of the eruption
characterized by a sudden rise in the kinetic energy and corresponding release
of the magnetic energy. We first briefly discussed the initiation and
evolution of the eruptions and then analyze the velocity and reconnection flux
of different CME eruptions in the simulation. Note that the magnetic
reconnection in our model using MHD approximation is of numerical origin, we
cannot resolve length scales smaller than the Larmor radius of the proton or
the mean free path. However, our goal here is to find a relation (if any
exists) between the reconnection flux and the speed of the ejecta in our
simulations.

###  4.1 Overview of eruptions

Figure 8: Total kinetic energy $E_{kin}$ (red) and total magnetic energy
$E_{mag}$ (blue) as a function of time. The dashed vertical line represents
the time when the flux emergence stops. Figure 9: Identification of current
sheet and field lines of the flux rope from different viewing angles. The
current sheet (in white) traced from the temperature isosurface with a value
$\log_{10}T=6.6$ shown in panel (a). The field lines passing through such
sheets gives the reconnection flux. The sigmoid fieldlines are shown in panels
(b) and (c). The time in all these panels are taken at $t=27.32$ hours.

The simulation starts with a magnetic flux rope (MFR) emerging from the lower
boundary and moving into an existing coronal potential arcade field. The
normal magnetic flux distributions at the lower boundary are represented by
bipolar bands $B_{r}(R_{\odot},\theta,\phi)=B_{s}(\theta)$ where the
$B_{s}(\theta)$ denotes the potential arcade field at the lower boundary, as
shown in Figure (7). The emergence process is halted once a specified amount
of twisted flux is driven into the corona. Although the imposed emergence of a
twisted flux rope through the lower boundary may not perfectly reflect real-
world conditions, it serves as a means to achieve a sequence of near-force-
free coronal flux rope equilibria with increasing amounts of locally detached,
twisted flux. In all cases, we observe the development of current layers with
a sigmoid morphology beneath the flux rope. Magnetic reconnection within these
current layers continues to add flux to the flux rope, even though the total
magnetic energy gradually decreases due to reconnections. This continuous
addition of flux facilitates successive eruptions.

Figure 10: Time evolution of velocity, (a)-(b), and reconnection flux,
(c)-(d), for both the eruptions in the simulation. The velocity is calculated
by tracking the flux rope in the central meridional plane. The reconnection
flux is calculated by the magnetic flux passing through the footpoints of the
fieldlines as well as through the current sheet shown in Fig. 9a. The red
(blue) circles are used to represent simulation data before (after) the
eruption.

The first eruption is triggered when the twist of the winding of the flux rope
field lines about the torus axis reaches approximately 1 wind between the
anchored foot points. We stopped the emergence of the twisted torus tube at
33.5 hours when the extent of the flux rope that emerged is such that the
field line’s winding about the flux rope’s axis is 1.5 between anchored
footpoints. This is estimated in a simple way by multiplying the total number
of winds around the circular torus, i.e., 5 winds, by the fraction of the
torus circumference above the solar surface. The second eruption began after
we stopped the flux emergence. The helical kink instability is expected to
develop in a line-tied coronal flux rope if the total windings of the field
line twist around the axis exceed a critical value between the line-tied ends
(Hood & Priest, 1981; Török & Kliem, 2003; Török et al., 2004). This critical
value is 1.25 (Hood & Priest, 1981) for uniformly cylindrical flux tube. The
animation (corresponding to Figure 7) clearly shows significant writhing
motion at the start of the eruption for the second eruption but not in first
eruption, indicative of the helical kink instability. Figure 8 illustrates the
evolution of kinetic and magnetic energies for both eruptions. The time at
which the kinetic energy reaches its peak is designated as the eruption time.

Furthermore, we examined whether the flux rope is also unstable to the torus
instability (Bateman, 1978; Kliem & Török, 2006; Isenberg & Forbes, 2007),
which is triggered by a vertical force imbalance between the upward hoop force
and the downward strapping force acting on the magnetic flux rope (MFR). The
external strapping field is assumed to be the potential field $B_{p}$ with the
same normal distribution at the lower boundary. The rate of decline with
height is measured by the decay index $n=-d\ln B_{h}/d\ln h$. The critical
decay index for a 3D line-tied arched flux rope has been calculated by
Isenberg & Forbes (2007) to be close to 1.7. In this work, for a decay index
of 1.7, the critical heights are 1.21 $R_{\odot}$ and 1.28 $R_{\odot}$,
respectively for the two successive eruptions. The height of the apex of the
axial field line of the coronal flux rope, estimating it to be 1.6 $R_{\odot}$
and 1.7 $R_{\odot}$ for the two successive CMEs at the time of eruption
calculated by tracking the flux rope over the course of their evolution (see
Appendix). For both eruptions, the flux rope crosses the critical threshold of
height, indicating that it is also torus unstable. This indicates the
fundamental role of kink instability and torus instability in triggering
ejective eruptions of a 3D line-tied coronal flux rope.

###  4.2 Evolution of velocity and reconnection flux

The panels (a) and (b) of Figure 10 illustrate the radial velocity of the flux
rope during both eruptions. These velocities were determined by tracking the
apex of each flux rope throughout their evolution (see Appendix). The first
ejection has a peak radial velocity of $224\,\mathrm{km\,s^{-1}}$, which is
slightly higher than the $213\,\mathrm{km\,s^{-1}}$ velocity of the second
CME. This observation is understood by the time versus kinetic energy plot,
which shows the second CME registered lower kinetic energy than its precursor.
The initial eruption occurs at a radial distance of $1.9R_{\odot}$, while the
second eruption takes place at $1.8R_{\odot}$. Due to the lower velocity of
the second eruption, there is no interaction between the two CMEs. An
interesting observation is that the velocity starts to decrease after the
eruption which is also observed in observational data. This could be
attributed to the absence of solar wind, which normally drives the magnetic
flux rope (MFR) further into the heliosphere post-eruption.

The panels (c) and (d) of Figure 10 depicts the changes in reconnection flux
over time for both eruptions. We monitored the hottest region beneath the flux
rope in the $r\textrm{-}\theta$ plane both before and after the eruption as
shown in Figure 9). This region corresponds to the current sheet, as the
current sheet typically exhibits higher temperatures compared to the
surrounding areas. By tracking the same hottest region, we can effectively
identify the current sheet and quantify the flux passing through it. The total
reconnection flux passing through the is given by

$\phi_{RC}=\int B_{n}\,dS=\int B_{r}\,dA$

where $B_{n}$ represents the magnetic flux normal to the elemental area $dS$
of the current sheet, and $B_{r}$ is the radial magnetic field in the ribbons,
which are the footpoints of newly reconnected magnetic field lines at the
lower boundary. The term $dA$ denotes the elemental area of these footpoints.
Thus reconnection flux is determined by calculating the magnetic flux swept by
the field lines passing the lower boundary. Direct measurement of $B_{n}$ and
$dS$ within the current sheet is not feasible, even in simulations, due to the
intricate structure of the reconnection sheet. However, $B_{r}$ and $dA$ are
relatively straightforward to calculate, even for a spherical wedge-shaped
domain. Therefore, this approach provides an indirect but well-defined measure
of the magnetic flux passing through the reconnection sheet. This method
allows for an effective evaluation of the magnetic reconnection process by
focusing on the measurable quantities at the footpoints, offering valuable
insights into the dynamics of the reconnection events.

Figure 11: (a) CME reconnection flux vs velocity plot for the first eruption.
The red (blue) circles are used to represent before (after) the eruption,
respectively. The solid (dashed) line is the linear fit to the points before
(after) the eruption, respectively. (b) Same as (a) but for the second
eruption.

The peak reconnection flux (RC flux) for the first and second eruptions are
$9\times 10^{20}$ Mx and $8.5\times 10^{20}$ Mx, respectively. The panels (c)
and (d) of Figure 10 indicate that the reconnection flux decreases after
reaching its peak value, with the eruptions occurring near the maximum
reconnection flux. This phenomenon might be explained by the intricacies of
calculating the RC flux. When an eruption occurs, the RC sheet breaks down,
making it challenging to track the field lines passing through the current
sheet (CS). Despite this complexity, similar patterns are observed in
observational data (see panel (b) of Figure 4, suggesting a consistent
underlying process. The observed decrease in reconnection flux post-peak might
be tied to the dynamic changes within the RC sheet during an eruption. As the
RC sheet disintegrates, accurately measuring and tracking the magnetic field
lines through the CS becomes challenging. This complexity impacts the
calculated values of RC flux, reflecting the real-time physical changes
occurring during an eruption. Moreover, these findings align with
observational data, further validating the temporal relation between the
reconnection flux and the onset of eruptions (panels b-d of Fig. 4). This
highlights the importance of considering both the peak and the subsequent
decrease in reconnection flux when studying the dynamics of solar eruptions
and magnetic reconnection events.

###  4.3 Correlation between MFR velocity and reconnection flux

The panels (a) and (b) of Figure 11 show scatter plots of the magnetic flux
rope (MFR) radial velocity ($v_{r}$) versus the reconnection flux
($\phi_{RC}$). Both before and after the eruption, the radial velocity
exhibits a linear relationship with the reconnection flux for the eruptive
events. Notably, the Pearson’s correlation coefficients before the eruption
for the two events are 0.86 and 0.96, respectively, indicating a strong
correlation between the reconnection flux and the radial velocity. This strong
linear relationship persists post-eruption, with the Pearson’s correlation
coefficients slightly decreasing to 0.82 and 0.91, respectively. This
indicates the significant influence of reconnection flux on the radial
velocity of CMEs. A similar qualitative relationship is also observed in the
case of observational data (see panel (d) of Figure 4) where the velocity of
the erupting filament and the associated reconnection flux simultaneously
increase as indicated by the cyan shaded region. Due to observational
constraints (limited time-cadence of STEREO EUVI), it was challenging to
gather more data points during the initial acceleration phase of the eruption.
After the peak acceleration phase, both the velocity and reconnection flux
exhibit a monotonically decreasing temporal profile while maintaining similar
correlated behavior. The numerical analysis from our simulation complements
the above-mentioned observational results, allowing us to determine the
temporal evolution of the reconnection flux and velocity and their
correlation.

##  5 Summary and Discussion

Through a simple and efficient MHD simulation, we demonstrated that homologous
CMEs can result from the repeated formation and partial eruption of unstable
flux ropes during the sustained emergence of highly twisted magnetic fields in
an active region. Our results show that, with the continuous emergence of
flux, a new current sheet (CS) forms above the same polarity inversion line
(PIL) after the completion of a previous eruption, leading to a new eruption.
This recurring formation and disruption of current sheet, driven by the
ongoing emergence of flux within a single bipolar configuration, results in
homologous eruptions.

We initiated our simulation with a twisted magnetic flux rope emerging from
the lower boundary into a pre-existing coronal arcade structure. Our results
revealed that the flux rope initially enters a phase of slow, quasi-static
ascent. During this phase, we observed the formation of a sigmoid-shaped
current layer beneath the flux rope. This current layer is significant as it
likely serves as the location where thin current sheets develop and magnetic
reconnection occurs (Aulanier et al., 2005; Titov, 2007; Savcheva et al.,
2012). The reconnection in the current layer effectively adds twisted flux to
the magnetic flux rope, allowing it to rise quasi-statically before undergoing
a partial eruption. Subsequent flux emergence forms a new current channel,
which exhibits similar behavior and erupts after the emergence stops.

In both eruptions, the magnetic flux rope (MFR) demonstrates torus instability
after reaching a height $\sim 1.6R_{\odot}$. The first eruption remains stable
to helical kink instability, while the second eruption is kink unstable. The
subsequent CME, in our simulation, has a slighly lower speed and lower kinetic
energy. Initially, the reconnection flux increases, followed by a gradual
decrease, after each eruption. These trends are consistent with the
observation of an eruptive event. Our simulation indicates that the
reconnection process significantly influences the evolution of the MFR. The
reconnection flux shows a linear correlation with the velocity, both before
and after the eruption. This result indicates that reconnection may play an
important role in propelling the magnetic ejecta after the loss of equilibrium
phase.

Although our numerical simulation model does not include the solar wind, the
inclusion of radiative cooling and explicit coronal heating makes it more
realistic than previously published homologous solar eruption models
(Chatterjee & Fan, 2013; Bian et al., 2022). We account for localized heating
due to the formation and dissipation of current sheets in the corona, as well
as the redistribution of heat through field-aligned conduction. The lower
boundary of our simulation domain is positioned at the coronal base rather
than the chromosphere, and therefore, our model does not simulate prominence
formation. In our computational models, we incorporated an adiabatic index of
$\gamma=5/3$, leads to more pronounced adiabatic cooling as the flux rope
expands. Additionally, we explicitly accounted for heating processes and heat
transport mechanisms. Despite these more stringent conditions, our results
show that the erupting flux rope can still attain velocities characteristic of
those observed in the lower corona.

Although we employ an analytical configuration of the magnetic flux rope (MFR)
and include various physical processes in our model, this approach may not
fully represent the realistic conditions of the solar eruptions. To address
the complexity and evolution of the magnetic configuration in a real scenario,
a detailed and accurate description of the evolving magnetic field is
essential. Our future goal is to conducting detailed analysis of solar
eruptions that are directly driven or constrained by photospheric magnetograms
is crucial for obtaining an accurate understanding.

The work was supported by Indo-US Science and Technology Forum
(IUSSTF/JC-113/2019). SSM acknowledge NOVA HPC and VINO server of IIA used to
perform the numerical simulation and analysis. P.C. acknowledges the DiRAC
Data Intensive service (DIaL2) at the University of Leicester (project id:
dp261), managed by the University of Leicester Research Computing Service on
behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at
Leicester was funded by BEIS, UKRI and STFC capital funding and STFC
operations grants. DiRAC is part of the UKRI Digital Research Infrastructure.
P.C. also acknowledges the allocation of computing resources at the PDC Center
for High Performance Computing at KTH in Stockholm, funded by the National
Academic Infrastructure for Supercomputing in Sweden. RS acknowledges support
from the project EFESIS (Exploring the Formation, Evolution and Space-weather
Impact of Sheath-regions), under the Academy of Finland Grant 350015.
Additionally, authors thank the NASA SDO team for providing valuable HMI and
AIA data. SDO is a mission under NASA’s Living with a Star program. Open
access is funded by Helsinki University Library.

## References

  * Amari et al. (2003) Amari, T., Luciani, J. F., Aly, J. J., Mikic, Z., & Linker, J. 2003, ApJ, 585, 1073, doi: 10.1086/345501
  * Antiochos et al. (1999) Antiochos, S. K., DeVore, C. R., & Klimchuk, J. A. 1999, ApJ, 510, 485, doi: 10.1086/306563
  * Archontis et al. (2009) Archontis, V., Hood, A. W., Savcheva, A., Golub, L., & Deluca, E. 2009, ApJ, 691, 1276, doi: 10.1088/0004-637X/691/2/1276
  * Aulanier (2014) Aulanier, G. 2014, in Nature of Prominences and their Role in Space Weather, ed. B. Schmieder, J.-M. Malherbe, & S. T. Wu, Vol. 300, 184–196, doi: 10.1017/S1743921313010958
  * Aulanier et al. (2005) Aulanier, G., Pariat, E., & Démoulin, P. 2005, A&A, 444, 961, doi: 10.1051/0004-6361:20053600
  * Aulanier et al. (2010) Aulanier, G., Török, T., Démoulin, P., & DeLuca, E. E. 2010, ApJ, 708, 314, doi: 10.1088/0004-637X/708/1/314
  * Bateman (1978) Bateman, G. 1978, MHD instabilities (The MIT Press)
  * Bein et al. (2012) Bein, B. M., Berkebile-Stoiser, S., Veronig, A. M., Temmer, M., & Vršnak, B. 2012, ApJ, 755, 44, doi: 10.1088/0004-637X/755/1/44
  * Bian et al. (2022) Bian, X., Jiang, C., Feng, X., Zuo, P., & Wang, Y. 2022, ApJ, 925, L7, doi: 10.3847/2041-8213/ac4980
  * Chatterjee (2020) Chatterjee, P. 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 213, doi: 10.1080/03091929.2019.1672676
  * Chatterjee & Fan (2013) Chatterjee, P., & Fan, Y. 2013, ApJ, 778, L8, doi: 10.1088/2041-8205/778/1/L8
  * Chatterjee et al. (2016) Chatterjee, P., Hansteen, V., & Carlsson, M. 2016, Phys. Rev. Lett., 116, 101101, doi: 10.1103/PhysRevLett.116.101101
  * Chen (2017) Chen, J. 2017, Physics of Plasmas, 24
  * Chen (2017) Chen, J. 2017, Physics of Plasmas, 24, 090501, doi: 10.1063/1.4993929
  * Chen (2011) Chen, P. F. 2011, Living Reviews in Solar Physics, 8, 1, doi: 10.12942/lrsp-2011-1
  * Cheng et al. (2011) Cheng, X., Zhang, J., Liu, Y., & Ding, M. D. 2011, ApJ, 732, L25, doi: 10.1088/2041-8205/732/2/L25
  * Cook et al. (1989) Cook, J. W., Cheng, C. C., Jacobs, V. L., & Antiochos, S. K. 1989, ApJ, 338, 1176, doi: 10.1086/167268
  * DeVore & Antiochos (2008) DeVore, C. R., & Antiochos, S. K. 2008, ApJ, 680, 740, doi: 10.1086/588011
  * Duan et al. (2019) Duan, A., Jiang, C., He, W., et al. 2019, ApJ, 884, 73, doi: 10.3847/1538-4357/ab3e33
  * Dun et al. (2007) Dun, J., Kurokawa, H., Ishii, T. T., Liu, Y., & Zhang, H. 2007, ApJ, 657, 577, doi: 10.1086/510373
  * Fan (2001) Fan, Y. 2001, ApJ, 554, L111, doi: 10.1086/320935
  * Fan (2010) —. 2010, ApJ, 719, 728, doi: 10.1088/0004-637X/719/1/728
  * Fan (2012) —. 2012, ApJ, 758, 60, doi: 10.1088/0004-637X/758/1/60
  * Forbes & Priest (1984) Forbes, T. G., & Priest, E. R. 1984, Sol. Phys., 94, 315, doi: 10.1007/BF00151321
  * Forbes et al. (2006) Forbes, T. G., Linker, J. A., Chen, J., et al. 2006, Space Sci. Rev., 123, 251, doi: 10.1007/s11214-006-9019-8
  * Gibson & Fan (2006) Gibson, S. E., & Fan, Y. 2006, Journal of Geophysical Research (Space Physics), 111, A12103, doi: 10.1029/2006JA011871
  * Gopalswamy et al. (2017) Gopalswamy, N., Yashiro, S., Akiyama, S., & Xie, H. 2017, Sol. Phys., 292, 65, doi: 10.1007/s11207-017-1080-9
  * Hood & Priest (1981) Hood, A. W., & Priest, E. R. 1981, Geophysical and Astrophysical Fluid Dynamics, 17, 297, doi: 10.1080/03091928108243687
  * Howard et al. (2008) Howard, R. A., Moses, J. D., Vourlidas, A., et al. 2008, Space Sci. Rev., 136, 67, doi: 10.1007/s11214-008-9341-4
  * Hu (2017) Hu, Q. 2017, Science China Earth Sciences, 60, 1466
  * Isenberg & Forbes (2007) Isenberg, P. A., & Forbes, T. G. 2007, ApJ, 670, 1453, doi: 10.1086/522025
  * Janvier et al. (2015) Janvier, M., Aulanier, G., & Démoulin, P. 2015, Sol. Phys., 290, 3425, doi: 10.1007/s11207-015-0710-3
  * Kazachenko (2023) Kazachenko, M. D. 2023, arXiv e-prints, arXiv:2310.02878, doi: 10.48550/arXiv.2310.02878
  * Kazachenko et al. (2017) Kazachenko, M. D., Lynch, B. J., Welsch, B. T., & Sun, X. 2017, ApJ, 845, 49, doi: 10.3847/1538-4357/aa7ed6
  * Kliem & Török (2006) Kliem, B., & Török, T. 2006, Phys. Rev. Lett., 96, 255002, doi: 10.1103/PhysRevLett.96.255002
  * Kopp & Poletto (1986) Kopp, R. A., & Poletto, G. 1986, in The Sun and the Heliosphere in Three Dimensions, ed. R. G. Marsden & L. A. Fisk, Vol. 123, 65, doi: 10.1007/978-94-009-4612-5_8
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17, doi: 10.1007/s11207-011-9776-8
  * Li et al. (2010) Li, Y., Lynch, B. J., Welsch, B. T., et al. 2010, Sol. Phys., 264, 149, doi: 10.1007/s11207-010-9547-y
  * Liu (2020) Liu, R. 2020, Research in Astronomy and Astrophysics, 20, 165
  * MacTaggart & Hood (2009) MacTaggart, D., & Hood, A. W. 2009, A&A, 508, 445, doi: 10.1051/0004-6361/200913197
  * Magara (2004) Magara, T. 2004, ApJ, 605, 480, doi: 10.1086/382148
  * Martínez-Sykora et al. (2008) Martínez-Sykora, J., Hansteen, V., & Carlsson, M. 2008, ApJ, 679, 871, doi: 10.1086/587028
  * Nitta & Hudson (2001) Nitta, N. V., & Hudson, H. S. 2001, Geophys. Res. Lett., 28, 3801, doi: 10.1029/2001GL013261
  * Patsourakos et al. (2013) Patsourakos, S., Vourlidas, A., & Stenborg, G. 2013, ApJ, 764, 125, doi: 10.1088/0004-637X/764/2/125
  * Patsourakos et al. (2020) Patsourakos, S., Vourlidas, A., Török, T., et al. 2020, Space Sci. Rev., 216, 131, doi: 10.1007/s11214-020-00757-9
  * Pencil Code Collaboration et al. (2021) Pencil Code Collaboration, Brandenburg, A., Johansen, A., et al. 2021, The Journal of Open Source Software, 6, 2807, doi: 10.21105/joss.02807
  * Qiu et al. (2004) Qiu, J., Wang, H., Cheng, C. Z., & Gary, D. E. 2004, ApJ, 604, 900, doi: 10.1086/382122
  * Ranns et al. (2000) Ranns, N. D. R., Harra, L. K., Matthews, S. A., & Culhane, J. L. 2000, A&A, 360, 1163
  * Régnier & Canfield (2006) Régnier, S., & Canfield, R. C. 2006, A&A, 451, 319, doi: 10.1051/0004-6361:20054171
  * Romano et al. (2018) Romano, P., Elmhamdi, A., Falco, M., et al. 2018, ApJ, 852, L10, doi: 10.3847/2041-8213/aaa1df
  * Romano et al. (2015) Romano, P., Zuccarello, F., Guglielmino, S. L., et al. 2015, A&A, 582, A55, doi: 10.1051/0004-6361/201525887
  * Sarkar et al. (2019) Sarkar, R., Srivastava, N., & Veronig, A. M. 2019, The Astrophysical Journal Letters, 885, L17, doi: 10.3847/2041-8213/ab4da2
  * Savcheva et al. (2012) Savcheva, A., Pariat, E., van Ballegooijen, A., Aulanier, G., & DeLuca, E. 2012, ApJ, 750, 15, doi: 10.1088/0004-637X/750/1/15
  * Schmieder et al. (2013) Schmieder, B., Démoulin, P., & Aulanier, G. 2013, Advances in Space Research, 51, 1967
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Shibata & Magara (2011) Shibata, K., & Magara, T. 2011, Living Reviews in Solar Physics, 8, 6, doi: 10.12942/lrsp-2011-6
  * Sterling & Moore (2001) Sterling, A. C., & Moore, R. L. 2001, J. Geophys. Res., 106, 25227, doi: 10.1029/2000JA004001
  * Temmer et al. (2008) Temmer, M., Veronig, A. M., Vršnak, B., et al. 2008, ApJ, 673, L95, doi: 10.1086/527414
  * Thompson et al. (2003) Thompson, W. T., Davila, J. M., Fisher, R. R., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4853, Innovative Telescopes and Instrumentation for Solar Astrophysics, ed. S. L. Keil & S. V. Avakyan, 1–11, doi: 10.1117/12.460267
  * Titov (2007) Titov, V. S. 2007, ApJ, 660, 863, doi: 10.1086/512671
  * Titov & Démoulin (1999) Titov, V. S., & Démoulin, P. 1999, A&A, 351, 707
  * Török & Kliem (2003) Török, T., & Kliem, B. 2003, A&A, 406, 1043, doi: 10.1051/0004-6361:20030692
  * Török et al. (2004) Török, T., Kliem, B., & Titov, V. S. 2004, A&A, 413, L27, doi: 10.1051/0004-6361:20031691
  * van Ballegooijen & Martens (1989) van Ballegooijen, A. A., & Martens, P. C. H. 1989, ApJ, 343, 971, doi: 10.1086/167766
  * Wang et al. (2016) Wang, Y., Zhuang, B., Hu, Q., et al. 2016, Journal of Geophysical Research (Space Physics), 121, 9316, doi: 10.1002/2016JA023075
  * Wuelser et al. (2004) Wuelser, J.-P., Lemen, J. R., Tarbell, T. D., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5171, Telescopes and Instrumentation for Solar Astrophysics, ed. S. Fineschi & M. A. Gummin, 111–122, doi: 10.1117/12.506877
  * Xu et al. (2017) Xu, Z., Yang, K., Guo, Y., et al. 2017, ApJ, 851, 30, doi: 10.3847/1538-4357/aa9995
  * Zhang et al. (2001) Zhang, J., Dere, K. P., Howard, R. A., Kundu, M. R., & White, S. M. 2001, ApJ, 559, 452, doi: 10.1086/322405
  * Zhang & Wang (2002) Zhang, J., & Wang, J. 2002, ApJ, 566, L117, doi: 10.1086/339660
  * Zhang et al. (2008) Zhang, Y., Liu, J., & Zhang, H. 2008, Sol. Phys., 247, 39, doi: 10.1007/s11207-007-9089-0

##  6 Appendix

The height-time plot is produced by tracking the dark region in the $r-\theta$
plane of the temperature profile, which represents the MFR. By following this
dark region consistently, we determine the MFR’s velocity during the flux
rope’s evolution.

Figure 12: (a) Height time profile of the flux rope before and after the first
eruption at $t=27.9$ hours. The height is calculated by tracking the dark
region (flux rope) in the temperature profile. (b) Same as (a), but for the
second eruption at $t=34.7$ hours.

The variation in reconnection flux over time is attributed to changes in the
area of the foot points of the field lines near the lower boundary at a height
of $1.04R_{\odot}$. This relationship is illustrated for both eruptions, with
the corresponding areas plotted over time in Figures 13. These plots provide a
detailed representation of how the swept area evolves and influences the
reconnection flux during each eruption.

Figure 13: Cumulative area of footpoints of the fieldlines passing through the
the reconnection sheet for the first (a) and second eruption (b),
respectively.


# Python4Physics:
A physics outreach program

R. A. Briceño1

1Department of Physics, University of California, Berkeley, CA 94720, USA
Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA
94720, USA
Email<EMAIL_ADDRESS>
ORCID: 0000-0003-1109-1473

T. C. Rogers2,3

2Department of Physics, Old Dominion University, Norfolk, VA 23529, USA
3Jefferson Lab, 12000 Jefferson Avenue, Newport News, VA 23606, USA
Email<EMAIL_ADDRESS>
ORCID: 0000-0002-0762-0275

We describe a summer outreach program developed to cultivate interest in
physics in particular and physical sciences more broadly among high school and
early college students using small projects in the Python programming
language. We discuss the lessons we learned in the hopes that they will be
valuable to other physicists in planning their own outreach efforts. We also
provide links to resources and materials from the Python4Physics program,
which we hope might be useful in other outreach programs.1[] We describe a
summer outreach program developed to cultivate interest in physics in
particular and physical sciences more broadly among high school and early
college students using small projects in the Python programming language. We
discuss the lessons we learned in the hopes that they will be valuable to
other physicists in planning their own outreach efforts. We also provide links
to resources and materials from the Python4Physics program, which we hope
might be useful in other outreach programs.1

##  1 Outreach to high schools

Professional research scientists frequently express a desire to cultivate a
greater interest in science in a wider population of young people. However, a
university or lab environment provides a limited set of opportunities to
interact with high school youth (and other prospective scientists) who are
truly in their most formative years. It is often during high school years that
many scientists are initially inspired to follow a path toward a career in
science, and it is clearly a missed opportunity that scientists rarely
interact directly with students at this stage. To address this, for the last
five years, we have been actively reaching out to local high schools, both
domestically and abroad, to cultivate greater interest in Science, Technology,
Engineering, and Math (STEM). To increase the likelihood that those high
school students will choose and stay in a STEM career, we created a free
course, called Python4Physics, to teach students the minimal scientific
background needed to perform computing-based research in a fun and accessible
way.

Our purpose with this document is to provide other researchers hoping to
develop their own broad-based outreach programs with resources and an account
of lessons learned from our own experiences.

##  2 Python4Physics

Our first Python4Physics course took place in person at Old Dominion
University during the summer of 2019. The basic structure of the course was
organized around lessons in the Python programming language. A small number of
introductory slides were presented, and then students gathered in groups to
solve fun but simple and open-ended problems on the computer. Examples of some
of the problems were projectile motion projects with graphing and animation,
elastic collisions in both nonrelativistic and relativistic scenarios, and the
fitting of a function in a manner akin to the identification of a particle
discovery in particle physics. All lectures were designed assuming the
students had not yet seen calculus, but we did assume some experience with
basic algebra.

By structuring the course around programming-based tasks, we aimed to initiate
students into scientific computer programming more generally, thereby
equipping them with tools that would be helpful in a wide range of potential
future scientific paths. At the same time, by emphasizing specific scientific
problems we hoped to give the students a chance to explore generally
scientific ideas in an open-ended way in the style of a practicing
professional scientist. This avoided having the course become simply another
“programming class.” We used the Python programming language because it is a
popular high-level language with a vast amount of free online resources and
support. Our guiding principle has been to maximize the motivational and
inspirational aspects of these projects rather than technical details.

In the summer of 2019, these meetings took place twice a week and lasted for
several hours. After a short introductory explanation of the task or project,
our role was simply to move around the room assisting student groups as they
worked. In this, we had help as well from volunteer undergraduate students,
graduate students, and postdoctoral researchers. In this way, the program also
provided learning opportunities for our early career researchers to develop
their own teaching, outreach, and communication skills.

In order to reduce our burden and the pressure on the students, we decided not
to provide grades on the projects. Instead, we gave them verbal feedback,
which was largely positive reaffirmation and constructive criticism. We also
provided them with our code that includes solutions to the projects. We
encouraged them to only look at our code _after_ they had made a strong
independent attempt to solve the problem.

In the summer of 2020, the global pandemic forced us to reorganize the program
around online activities, conducted via YouTube, and this impacted certain
aspects of the course itself. A positive outcome was that it allowed us to
expand the reach of the program to a much larger (100-300 students per year),
and more international audience (30-70 countries represented per year). We
also expanded the scope of topics that we discussed to include more general
scientific and computational topics outside of physics, including topics
motivated by current events such as the basics of modeling epidemics. The
popularity of Python4Physics in Latin American was sufficiently large, that
with the help of postdocs and graduate students, in the summer of 2021, we
hosted a Spanish version of the class, called Python y Ciencias.111We were
able to host the program remotely and reach large international audiences
thanks to the help of the Remote Experience for Young Engineers and Scientists
(REYES) program.

The largest drawback of hosting these classes remotely was the lack of
personal interaction with the students. It was also challenging to imitate the
hands-on environment that we had in person. This is motivating us and others
to return to an in-person, or at least hybrid, environment.222For more
information about the hybrid Python4Physics class being hosted in the summer
of 2024 by the University of California, Berkeley Physics Department,
including the format of the registration page, go the official website. The
hybrid model promises to provide a more personalized experience for those who
can come to our institutions in person, while still leaving the course open
and accessible to all students from around the world to join in remotely to
our classrooms.

##  3 Obstacles and challenges

In organizing and preparing the Python4Physics course, we encountered several
challenges, and we believe that future planners of outreach programs might
benefit from a checklist of some of these. In this section, we enumerate some
of the more significant difficulties we faced and provide suggestions for
confronting them.

###  3.1 Advertising and attracting students

Since many university researchers rarely interact with students outside of the
university system, it can be challenging to find high school participants. We
largely benefited from the experience of staff in our departments and
institutions, who were able to connect us with local high school educators.
Through these contacts, we were able to either directly or indirectly reach
out to high school participants. In several instances, we took the time to go
to those high schools and give “general audience" talks about our lines of
research in front of students. As long as it is done in an accessible and
enjoyable fashion, this personal touch can go a long way to make students
comfortable in coming to your classroom.

###  3.2 Technical resources and
accessibility

Ensuring access to computers and software is of course a challenge for any
outreach course that involves programming.

Given that one of our goals is to make STEM as accessible as possible, we have
been mindful that our class requires the students to have access to a
computer. Within the USA, many High Schools now provide Chrome Books or
similar low-cost laptops during the academic year. Depending on the school
district, these laptops may or may not be accessible to students during the
summer. Ideally, the institutions hosting an outreach activity like
Python4Physics can provide these laptops during the summer, but depending on
the school district, this may not be possible.

For our first summer with in-person students, we used a stockpile of notebook
computers that we previously used in the university for lab classes. This
worked to a degree, but we encountered some difficulties with installing the
necessary software, and a fraction of the computers failed to function for the
full course.

This is an ongoing challenge for which we have not yet found an ideal
solution. Ideally, external private investors could be sought to provide funds
to ensure access to laptops.333In preparation for the first hybrid
Python4Physics course to be hosted in the summer of 2024, we have asked our
participants if they have access to a computer in our registration page.

###  3.3 Working with minors

Because our target audience is comprised of high school students, most
participants are under 18 years of age, and this comes with a number of
important legal considerations. Our experience has been that institutions like
Old Dominion University and the University of California, Berkeley, have clear
guidelines and best practices that one must follow when hosting an in-person
event. In our case, everyone who was to interact with the students was
required to undergo a training and criminal background check. This included
our undergraduate and graduate assistants. To avoid delays, it is important to
start this process early. Depending on the institution, there may be other
additional restrictions, e.g. two or more adults that have gone to said
training/checks may need to be present with students at all times.

Other requirements for working with minors may not be as obvious to those who
are more familiar with a college campus environment. For example, depending on
the locality, it is prohibited to include minors in a _cc_ on a mass email.
For sending out course information, therefore, it is important to use _bcc_ as
a default. We suggest communicating with your local institution or
institutions very early about its own guidelines for interacting with minors
so that preparations can be made.

When it came to working with minors remotely, at the beginning of the pandemic
we found the guidelines to be less clear. As a result, in collaboration with
the REYES team, we opted to follow the safest possible option. We hosted all
classes via Zoom Webinar. Zoom Webinar allows one to have two distinct
environments where participants can be placed, e.g. “ _Panel_ " or “
_Audience_ ". Depending on where the participant is placed, they may be
restricted in how they interact with others. This allowed us to ensure that
nobody could communicate directly with students, and it allowed us to control
who may activate their cameras or communicate with others. This provided a
remarkably safe experience for students, and it made it impossible for anybody
to “ _Zoombomb_ ” others with any profanity or unwanted behavior.444 Note,
that other companies have similar Webinar platforms that may be more
affordable, e.g. BlueJeans.

###  3.4 Time sink for hosts

It is inevitable that a program of the scope of Python4Physics will take time
away from faculty who are already spread quite thin. We find that the major
time sinks include the time required in developing the material, hosting the
actual classes, and miscellaneous administrative burdens.

Although some time costs are unavoidable, we are happy to provide the
following tips and resources to help reduce the burden:

  * •

Developing material: We are making all of the material publicly available.
This also includes our code that solves numerical tasks given to the students.
These code solutions are also shared for the students to self-evaluate their
solutions. You are welcome to copy and edit it as you wish. Please find the
materials at our project website here.

  * •

Hosting class: We have found that co-hosting classes reduces the burden on
each person substantially. In the last couple of years, we have opted to
create teams of volunteer faculty/staff, postdocs and graduate students to
help host the classes. Depending on the number of lecturers, each person may
have to lecture as little as one or two classes.

  * •

Administrative burden: Although we cannot remove all administrative burden
from those of you interested in hosting this, we are happy to coordinate with
you one of your efforts to host the class and provide as much service as
possible. An example of how we may be able to help is by sharing our
registration pages. Other examples can be discussed in a case-by-case basis.

Although here we list the burdens associated with the course, we would like to
emphasize that the experience of working with aspiring young scientists
through Python4Physics has been exclusively positive and fulfilling for the
instructors.

##  4 Lessons and other suggestions

Observations we have made while conducting the Python4Physics program have led
us to create our own set of guidelines that we feel help us to provide the
most fulfilling experience possible for students. While these are only based
on our personal experiences, and others may find better ways to achieve their
outreach goals, we hope that an account of our own experiences may be useful.

First, we found that the computational and programming aspects are
particularly intimidating to many students, and special efforts to make the
programming itself less formidable are always an improvement. Furthermore, we
found that we had spent excessive amounts of time on installing software in
the earliest versions of the course. It is our experience that any measures
that reduce the discussion of specific coding or software tasks will improve
the quality of the experience for students. While _exposure_ to computers and
programming is important, it is preferable to use the limited time and
attention available to cultivate broader scientific interests than to focus on
technical details. _By motivating and helping students gain some confidence,
they will be more likely to dedicate the extra time needed to learn the
technical details needed in due time. Think of yourself more like a
motivational coach, rather than a university professor._

The need for specialized software (compilers, interpreters, etc) can be
minimized by using online coding platforms, so that access to a web browser is
all that is necessary. For teaching programming, we typically utilize a single
short lecture to explain the basic structure of a generic program. There, we
briefly introduce concepts like _if_ statements and _for_ loops, but with very
little detail. An effective way to ensure that most students may take part in
some amount of programming without being stuck or left behind is to provide
them with nearly completed codes, with only one or a few pieces that need to
be finished. Fully completed codes are also provided so that they may check
their own work. This provides the majority of students with an experience of
“coding” in the very short time that is typically available.

These remarks concerning computer programming also apply to mathematical
reasoning more generally. Whenever equations need to be used for a task, the
exact final form that is necessary should be provided, in a clear and easy-to-
use format, without the need for any manipulations.

We find that student experiences always seem to improve when we maximize the
humor, drama, and visualization that we inject into the specific projects.
Current events or popular topics in online discourse provide valuable source
material for this. For example, we assigned projects to model epidemics, and
we taught some basics of cryptography, again providing examples and sample
codes. We provided guidance for visualizing interesting shapes like the
Fibonacci spiral. Many resources for fun and whimsical computer projects
already exist, and it is worthwhile to consider these.

Students who join the program come with a very wide range of backgrounds in
programming, mathematics, and science. By organizing the course as described
above, students are able to tailor the experience to their individual levels.
Those seeking a greater challenge may choose to build programs that tackle of
their own to tackle our assigned projects from scratch, while those
encountering the topics for the first time have a framework to build upon. We
also encourage a healthy collaborative environment. Among the ways we do this,
includes asking more experienced students to help those are coding for the
very first time.

##  5 Summary

Our goal is for all students to leave with experiences they can use as
starting points for future scientific endeavors. Each year we ask students to
provide feedback of their experience. Each year, we get that no less than 70%
“ _felt more confident and enthusiastic in pursuing a career in STEM and
conducting research_.

In summary, we strongly believe that Python4Physics is a powerful method to
inspire the next generation of scientists. The class requires a minimal amount
of funds, and it is easily reproducible. We are sharing our material and
experience, with the hopes that others may build from it.

## Acknowledgments

RAB was supported in part by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics under Awards No. DE-AC02-05CH11231. TCR was
supported by the U.S. Department of Energy, Office of Science, Office of
Nuclear Physics, under Award Number DE-SC0024715. This work was also supported
by the DOE Contract No. DE- AC05-06OR23177, under which Jefferson Science
Associates, LLC operates Jefferson Lab.


# Multilevel Iteration Method for Binary Stochastic Transport Problems

Dmitriy Y. Anistratov  Department of Nuclear Engineering, North Carolina State
University Raleigh, NC<EMAIL_ADDRESS>
###### Abstract

This paper presents an iteration method for solving linear particle transport
problems in binary stochastic mixtures. It is based on nonlinear projection
approach. The method is defined by a hierarchy of equations consisting of the
high-order transport equation for materials, low-order Yvon-Mertens equations
for conditional ensemble average of the material partial scalar fluxes, and
low-order quasidiffusion equations for the ensemble average of the scalar flux
and current. The multilevel system of equations is solved by means of an
iterative algorithm with the $V$-cycle. The iteration method is analyzed on a
set of numerical test problems.

###### keywords:

radiation transport , particle transport equation, stochastic media ,
Levermore-Pomraning closure , iteration methods , multilevel methods

##  1 Introduction

We consider the Levermore-Pomraning (LP) particle transport model for binary
stochastic mixtures (BSM) in 1D slab geometry [1, 2, 3]. The materials are
randomly distributed as alternating layers with homogeneous Markov mixing
statistics. The mean width of layers of material $\ell$ is equal to
$\lambda_{\ell}$. The BSM transport equation with the LP closure in 1D slab
geometry is given by

$\mu\frac{\partial}{\partial
x}\Big{(}p_{\ell}\psi_{\ell}(x,\mu)\Big{)}+\sigma_{t,\ell}p_{\ell}\psi_{\ell}(x,\mu)=\frac{1}{2}\sigma_{s,\ell}p_{\ell}\int_{-1}^{1}\psi_{\ell}(x,\mu^{\prime})d\mu^{\prime}\\\
+\frac{|\mu|}{\lambda_{\ell^{\prime}}}p_{\ell^{\prime}}\psi_{\ell^{\prime}}(x,\mu)-\frac{|\mu|}{\lambda_{\ell}}p_{\ell}\psi_{\ell}(x,\mu)+\frac{1}{2}p_{\ell}q_{\ell}(x)\,,$
(1) $\ell=1,2\,,\quad\ell^{\prime}\neq\ell\,,\quad
x\in[0,X]\,,\quad\mu\in[-1,1]\,$

with boundary conditions (BCs)

$\psi_{\ell}\Big{|}_{\stackrel{{\scriptstyle\mu>0}}{{x=0}}}=\psi_{\ell}^{in+}\,,\quad\psi_{\ell}\Big{|}_{\stackrel{{\scriptstyle\mu<0}}{{x=X}}}=\psi_{\ell}^{in-}\,,$
(2)

where $x$ is the spatial position; $\mu$ is the directional cosine of particle
motion; $\psi_{\ell}(x,\mu)$ is the conditional ensemble average (CEA) of the
angular flux provided that $x$ is in material $\ell$. $p_{\ell}$ is the
probability that the $\ell$-th material is found at a position in the spatial
domain. For homogeneous mixing statistics, the probability is given by

$p_{\ell}=\frac{\lambda_{\ell}}{\lambda_{1}+\lambda_{2}}\,.$ (3)

The total ensemble average of the angular flux is defined by

$\big{<}\psi\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}p_{1}\psi_{1}+p_{2}\psi_{2}\,.$
(4)

The CEA of the scalar flux and current for material $\ell$ are given by

$\phi_{\ell}\stackrel{{\scriptstyle\Delta}}{{=}}\int_{-1}^{1}\psi_{\ell}d\mu\,,\quad
J_{\ell}\stackrel{{\scriptstyle\Delta}}{{=}}\int_{-1}^{1}\mu\psi_{\ell}d\mu\,.$
(5)

The problems with randomly mixed materials arise in various application, for
example, in modeling of inertial confinement fusion targets, nuclear fuel and
shielding materials, clouds in atmospheric sciences, tissue in radiation
treatment planning [4]. The source iteration scheme for the linear binary
stochastic transport equation converges very slow [5, 6]. Various synthetic
acceleration methods for this class of particle transport problems have been
developed [7, 8, 9, 10, 11]. In this paper, we apply nonlinear projection
approach to derive a multilevel iteration method [12, 13]. It is defined by a
system of equations consisting of the high-order transport equation and low-
order equations for various moments of the transport solution. The iteration
method can be interpreted as a nonlinear multigrid method over elements of the
phase space. The proposed multilevel method can be applied for radiation
transport problems in random media coupled with multiphysics equations [14,
15].

The reminder of the paper is organized as follows. In Sec. 2 the system of
low-order equations for transport problems in BSM is formulated. The iteration
algorithm is described in Sec. 3. Numerical results are presented in Sec. 4.
We conclude with a brief discussion in Sec. 5.

##  2  Hierarchy of Equations of Multilevel Method for BSM

The multilevel method for solving the BSM transport problem (1) and (2) is
formulated by a high-order equation for $\psi_{\ell}$ and low-order equations
for moments the CEA of the angular flux with exact closures. To derive the
hierarchy of equations of the method we define (i) projection operators at
each level, (ii) moments of the transport solution low-order equations are
formulated for, and (iii) prolongation operators that enable to couple
equations at different levels.

A specific feature of the transport equation (1) for each material is that its
right-hand side (RHS) contains terms with absolute value of $\mu$ while the
left-hand side (LHS) explicitly depends of $\mu$. A natural way of forming a
set of moment equations based on projection of the high-order equation with
this structure of terms is to define the low-order equations for half-range
moments

$\phi_{\ell}^{\pm}\stackrel{{\scriptstyle\Delta}}{{=}}\pm\int_{0}^{\pm
1}\psi_{\ell}d\mu$ (6)

that are the CEA of the partial scalar fluxes. The high-order equation (1) is
integrated over $-1\leq\mu\leq 1$ with weights 1 and $\mu$ to obtain the low-
order equations for $\phi_{\ell}^{\pm}$ defined by

$\frac{d}{dx}\big{(}p_{\ell}(C_{\ell}^{+}\phi_{\ell}^{+}+C_{\ell}^{-}\phi_{\ell}^{-})\big{)}+\Big{(}\sigma_{a,\ell}+\frac{C_{\ell}^{+}}{\lambda_{\ell}}\Big{)}p_{\ell}\phi_{\ell}^{+}+\Big{(}\sigma_{a,\ell}-\frac{C_{\ell}^{-}}{\lambda_{\ell}}\Big{)}p_{\ell}\phi_{\ell}^{-}=\\\
\frac{p_{\ell^{\prime}}}{\lambda_{\ell^{\prime}}}\big{(}C_{\ell^{\prime}}^{+}\phi_{\ell^{\prime}}^{+}-C_{\ell^{\prime}}^{-}\phi_{\ell^{\prime}}^{-}\big{)}+p_{\ell}q_{\ell}\,,$
(7a)
$\frac{d}{dx}\big{(}p_{\ell}(E_{\ell}^{+}\phi_{\ell}^{+}+E_{\ell}^{-}\phi_{\ell}^{-})\big{)}+\Big{(}\sigma_{t,\ell}C_{\ell}^{+}+\frac{E_{\ell}^{+}}{\lambda_{\ell}}\Big{)}p_{\ell}\phi_{\ell}^{+}+\Big{(}\sigma_{t,\ell}C_{\ell}^{-}-\frac{E_{\ell}^{-}}{\lambda_{\ell}}\Big{)}p_{\ell}\phi_{\ell}^{-}=\\\
\frac{p_{\ell^{\prime}}}{\lambda_{\ell^{\prime}}}\big{(}E_{\ell^{\prime}}^{+}\phi_{\ell^{\prime}}^{+}-E_{\ell^{\prime}}^{-}\phi_{\ell^{\prime}}^{-}\big{)}\,,$
(7b)

where the low-order equations are closed by means of the linear-fractional
factors defined by

$C_{\ell}^{\pm}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\int_{0}^{\pm
1}\mu\psi_{\ell^{\prime}}d\mu}{\int_{0}^{\pm 1}\psi_{\ell}d\mu}\,,\quad
E_{\ell}^{\pm}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\int_{0}^{\pm
1}\mu^{2}\psi_{\ell}d\mu}{\int_{0}^{\pm 1}\psi_{\ell}d\mu}\,.$ (8)

The CEA of the material partial currents are defined by

$J^{\pm}_{\ell}\stackrel{{\scriptstyle\Delta}}{{=}}\pm\int_{0}^{\pm
1}\mu\psi_{\ell}d\mu\,.$ (9)

They are cast in terms of the partial scalar fluxes as follows:

$J^{\pm}_{\ell}=C^{\pm}_{\ell}\phi^{\pm}_{\ell}\,.$ (10)

The BCs for Eqs. (7) have the following form:

$\phi_{\ell}^{+}\Big{|}_{x=0}=\phi_{\ell}^{in,+}\,,\quad\phi_{\ell}^{-}\Big{|}_{x=X}=\phi_{\ell}^{in,-}\,,$
(11)

where

$\phi_{\ell}^{in,\pm}\\!\stackrel{{\scriptstyle\Delta}}{{=}}\\!\pm\int_{0}^{\pm
1}\\!\psi_{\ell}^{in,+}d\mu\,.$ (12)

The low-order equations (7) use the same projection operators and unknowns as
the DP1 equations [16, 17, 18, 19]. They can be viewed as nonlinear DP1
equations with exact closures. This approach has been applied to formulate
Yvon-Mertens (YM) method for the linear Boltzmann equation [20]. Hereafter
these equations are referred to as the low-order YM (LOYM) equations.

The BSM transport equation (1) is now projected to formulate low-order
equations for the total ensemble averaged scalar flux $\big{<}\phi\big{>}$ and
current $\big{<}J\big{>}$, where
$\big{<}f\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}p_{1}f_{1}+p_{2}f_{2}$ is
the notation for the ensemble average. The equation (1) is summed over
materials and then integrated over $-1\leq\mu\leq 1$ with weights 1 and $\mu$
to obtain the low-order quasidiffusion (LOQD) equations for
$\big{<}\phi\big{>}$ and $\big{<}J\big{>}$ given by [21]

$\frac{d\big{<}J\big{>}}{dx}+\big{<}\\!\\!\big{<}\sigma_{a}\big{>}\\!\\!\big{>}\big{<}\phi\big{>}=\big{<}q\big{>}\,,$
(13a)
$\frac{d}{dx}\Big{(}\big{<}\\!\\!\big{<}E\big{>}\\!\\!\big{>}\big{<}\phi\big{>}\Big{)}+\big{<}\\!\\!\big{<}\sigma_{t}\big{>}\\!\\!\big{>}\big{<}J\big{>}+\big{<}\\!\\!\big{<}\eta\big{>}\\!\\!\big{>}\big{<}\phi\big{>}=0\,,$
(13b)
$\big{<}J\big{>}\Big{|}_{x=0}\\!\\!=\\!\\!\big{<}\\!\\!\big{<}C^{-}\big{>}\\!\\!\big{>}\Big{(}\big{<}\phi\big{>}-\big{<}\phi^{in,+}\big{>}\Big{)}\Big{|}_{x=0}\\!\\!+\big{<}J^{in,+}\big{>}\,,$
(14a)
$\big{<}J\big{>}\Big{|}_{x=X}\\!\\!=\\!\\!\big{<}\\!\\!\big{<}C^{+}\big{>}\\!\\!\big{>}\Big{(}\big{<}\phi\big{>}-\big{<}\phi^{in,-}\big{>}\Big{)}\Big{|}_{x=X}\\!\\!+\big{<}J^{in,-}\big{>}\,,$
(14b)

where the coefficients are defined by means of the LOYM solution as follows:

$\big{<}\\!\\!\big{<}\sigma_{a}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}\sigma_{a,\ell}p_{\ell}(\phi_{\ell}^{+}+\phi_{\ell}^{-})}{\sum_{\ell=1}^{2}p_{\ell}(\phi_{\ell}^{+}+\phi_{\ell}^{-})}\,,\quad\big{<}\\!\\!\big{<}\sigma_{t}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}\sigma_{t,\ell}p_{\ell}|C_{\ell}^{+}\phi_{\ell}^{+}+C_{\ell}^{-}\phi_{\ell}^{-}|}{\sum_{\ell=1}^{2}p_{\ell}|C_{\ell}^{+}\phi_{\ell}^{+}+C_{\ell}^{-}\phi_{\ell}^{-}|}\,,$
(15)
$\big{<}\\!\\!\big{<}\eta\big{>}\\!\\!\big{>}\\!\stackrel{{\scriptstyle\Delta}}{{=}}\\!\frac{\sum_{\ell=1}^{2}\\!\big{(}\sigma_{t,\ell}\\!-\\!\big{<}\\!\\!\big{<}\sigma_{t}\big{>}\\!\\!\big{>}\big{)}p_{\ell}(C_{\ell}^{+}\phi_{\ell}^{+}\\!+\\!C_{\ell}^{-}\phi_{\ell}^{-})}{\sum_{\ell=1}^{2}p_{\ell}(\phi_{\ell}^{+}+\phi_{\ell}^{-})}\,,\quad\big{<}\\!\\!\big{<}E\big{>}\\!\\!\big{>}\\!\stackrel{{\scriptstyle\Delta}}{{=}}\\!\frac{\sum_{\ell=1}^{2}\\!p_{\ell}(E_{\ell}^{+}\phi_{\ell}^{+}+E_{\ell}^{-}\phi_{\ell}^{-})}{\sum_{\ell=1}^{2}p_{\ell}(\phi_{\ell}^{+}+\phi_{\ell}^{-})}\,,$
(16)
$\big{<}\\!\\!\big{<}C^{\pm}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}C_{\ell}^{\pm}p_{\ell}\phi_{\ell}^{\pm}}{\sum_{\ell=1}^{2}p_{\ell}\phi_{\ell}^{\pm}}\,,\quad
J_{\ell}^{in\pm}\stackrel{{\scriptstyle\Delta}}{{=}}\int_{0}^{\pm
1}\mu\psi_{\ell}^{in\pm}d\mu\,.$ (17)

The next step in formulation of the multilevel system of low-order equations
is to define prolongation operators that enable to couple the LOQD and LOYM
equations. The term
$\frac{p_{\ell}}{\lambda_{\ell^{\prime}}}(C_{\ell}^{+}\phi_{\ell}^{+}-C_{\ell}^{-}\phi_{\ell}^{-})$
is added to the LHS and RHS of Eq. (7a) and
$\frac{p_{\ell}}{\lambda_{\ell^{\prime}}}(E_{\ell}^{+}\phi_{\ell}^{+}-E_{\ell}^{-}\phi_{\ell}^{-})$
is added to the LHS and RHS of Eq. (7b) to get

$\frac{d}{dx}\big{(}p_{\ell}(C_{\ell}^{+}\phi_{\ell}^{+}+C_{\ell}^{-}\phi_{\ell}^{-})\big{)}+\bigg{(}\sigma_{a,\ell}+C_{\ell}^{+}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{+}+\bigg{(}\sigma_{a,\ell}-C_{\ell}^{-}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{-}=\\\
\frac{1}{\lambda_{\ell^{\prime}}}\sum_{\ell^{\prime\prime}=1}^{2}p_{\ell^{\prime\prime}}\Big{(}C_{\ell^{\prime\prime}}^{+}\phi_{\ell^{\prime\prime}}^{+}-C_{\ell^{\prime\prime}}^{-}\phi_{\ell^{\prime\prime}}^{-}\Big{)}+p_{\ell}q_{\ell}\,,$
(18a)
$\frac{d}{dx}\big{(}p_{\ell}(E_{\ell}^{+}\phi_{\ell}^{+}+E_{\ell}^{-}\phi_{\ell}^{-})\big{)}+\bigg{(}\sigma_{t,\ell}C_{\ell}^{+}+E_{\ell}^{+}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{+}+\bigg{(}\sigma_{t,\ell}C_{\ell}^{-}-E_{\ell}^{-}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{-}=\\\
\frac{1}{\lambda_{\ell^{\prime}}}\sum_{\ell^{\prime\prime}=1}^{2}p_{\ell^{\prime}}\Big{(}E_{\ell^{\prime\prime}}^{+}\phi_{\ell^{\prime\prime}}^{+}-E_{\ell^{\prime\prime}}^{-}\phi_{\ell^{\prime\prime}}^{-}\Big{)}\,.$
(18b)

Taking into account Eq. (10), we obtain

$\sum_{\ell^{\prime\prime}=1}^{2}p_{\ell^{\prime\prime}}\Big{(}C_{\ell^{\prime\prime}}^{+}\phi_{\ell^{\prime\prime}}^{+}-C_{\ell^{\prime\prime}}^{-}\phi_{\ell^{\prime\prime}}^{-}\Big{)}=\big{<}J^{+}\big{>}-\big{<}J^{-}\big{>}\,.$
(19)

We define

$\big{<}\\!\\!\big{<}E^{\pm}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}E_{\ell}^{\pm}p_{\ell}\phi_{\ell}^{\pm}}{\sum_{\ell=1}^{2}p_{\ell}\phi_{\ell}^{\pm}}\,$
(20)

to get

$\sum_{\ell^{\prime\prime}=1}^{2}p_{\ell^{\prime\prime}}\Big{(}E_{\ell^{\prime\prime}}^{+}\phi_{\ell^{\prime\prime}}^{+}-E_{\ell^{\prime\prime}}^{-}\phi_{\ell^{\prime\prime}}^{-}\Big{)}=\big{<}\\!\\!\big{<}E^{+}\big{>}\\!\\!\big{>}\big{<}\phi^{+}\big{>}-\big{<}\\!\\!\big{<}E^{-}\big{>}\\!\\!\big{>}\big{<}\phi^{-}\big{>}\,.$
(21)

Using Eqs. (19) and (21) in Eq. (18), we obtain the modified LOYM equations of
the following form:

$\frac{d}{dx}\big{(}p_{\ell}(C_{\ell}^{+}\phi_{\ell}^{+}+C_{\ell}^{-}\phi_{\ell}^{-})\big{)}+\bigg{(}\sigma_{a,\ell}+C_{\ell}^{+}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{+}+\bigg{(}\sigma_{a,\ell}-C_{\ell}^{-}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{-}=\\\
\frac{1}{\lambda_{\ell^{\prime}}}\Big{(}\big{<}J^{+}\big{>}-\big{<}J^{-}\big{>}\Big{)}+p_{\ell}q_{\ell}\,,$
(22a)
$\frac{d}{dx}\big{(}p_{\ell}(E_{\ell}^{+}\phi_{\ell}^{+}+E_{\ell}^{-}\phi_{\ell}^{-})\big{)}+\bigg{(}\sigma_{t,\ell}C_{\ell}^{+}+E_{\ell}^{+}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{+}+\bigg{(}\sigma_{t,\ell}C_{\ell}^{-}-E_{\ell}^{-}\Big{(}\frac{1}{\lambda_{\ell}}+\frac{1}{\lambda_{\ell^{\prime}}}\Big{)}\bigg{)}p_{\ell}\phi_{\ell}^{-}=\\\
\frac{1}{\lambda_{\ell^{\prime}}}\Big{(}\big{<}\\!\\!\big{<}E^{+}\big{>}\\!\\!\big{>}\big{<}\phi^{+}\big{>}-\big{<}\\!\\!\big{<}E^{-}\big{>}\\!\\!\big{>}\big{<}\phi^{-}\big{>}\Big{)}\,.$
(22b)

The prolongation operators that enable one to map $\big{<}\phi\big{>}$ and
$\big{<}J\big{>}$ to $\big{<}\phi^{\pm}\big{>}$ and $\big{<}J^{\pm}\big{>}$
are defined as follows [22]:

$\big{<}\phi^{\pm}\big{>}=\big{<}\\!\\!\big{<}\beta^{\pm}\big{>}\\!\\!\big{>}\Big{(}\big{<}J\big{>}-\big{<}\\!\\!\big{<}\tilde{C}^{\pm}\big{>}\\!\\!\big{>}\big{<}\phi\big{>}\Big{)}\,,\quad\big{<}J^{\pm}\big{>}=\big{<}\\!\\!\big{<}\gamma^{\pm}\big{>}\\!\\!\big{>}\Big{(}\big{<}J\big{>}-\big{<}\\!\\!\big{<}\tilde{C}^{\pm}\big{>}\\!\\!\big{>}\big{<}\phi\big{>}\Big{)}\,,$
(23)

where

$\big{<}\\!\\!\big{<}\beta^{\pm}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}p_{\ell}\phi_{\ell}^{\pm}}{\sum_{\ell=1}^{2}(C_{\ell}^{\pm}-C_{\ell}^{\mp})p_{\ell}\phi_{\ell}^{\pm}}\,,\quad\big{<}\\!\\!\big{<}\gamma^{\pm}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}C_{\ell}^{\pm}p_{\ell}\phi_{\ell}^{\pm}}{\sum_{\ell=1}^{2}(C_{\ell}^{\pm}-C_{\ell}^{\mp})p_{\ell}\phi_{\ell}^{\pm}}\,,$
(24)

and

$\big{<}\\!\\!\big{<}\tilde{C}^{\pm}\big{>}\\!\\!\big{>}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\sum_{\ell=1}^{2}C_{\ell}^{\pm}p_{\ell}(\phi_{\ell}^{+}+\phi_{\ell}^{-})}{\sum_{\ell=1}^{2}p_{\ell}(\phi_{\ell}^{+}+\phi_{\ell}^{-})}\,.$
(25)

To formulate the high-order transport equation of the multilevel method, we
cast the RHS of Eq. (1) in terms of $\phi_{\ell}$ to obtain the equation for
$\psi_{\ell}$ given by

$\mathcal{L}_{\ell}p_{\ell}\psi_{\ell}=\frac{1}{2}\sigma_{s,\ell}p_{\ell}\phi_{\ell}+\frac{|\mu|}{\lambda_{\ell^{\prime}}}p_{\ell^{\prime}}\psi_{\ell^{\prime}}+\frac{1}{2}p_{\ell}q_{\ell}\,,\
\ell^{\prime}\neq\ell\,,$ (26)
$\mathcal{L}_{\ell}\stackrel{{\scriptstyle\Delta}}{{=}}\mu\frac{\partial}{\partial
x}+\Big{(}\sigma_{t,\ell}+\frac{|\mu|}{\lambda_{\ell}}\Big{)}\,.$ (27)

##  3  Multilevel Iteration Scheme

The multilevel iteration method for transport problems in BSM is defined by
several stages which involve solving the following equations at different
levels:

  * 1.

the high-order transport equation (26) for $\psi_{\ell}$,

  * 2.

the LOYM equations (7) and (11) for $\phi^{\pm}_{\ell}$,

  * 3.

the LOQD equations (13) and (14) for $\big{<}\phi\big{>}$ and
$\big{<}J\big{>}$,

  * 4.

the modified LOYM equations (22) and (11) for $\phi^{\pm}_{\ell}$.

This method is a nonlinear multigrid scheme over elements of the phase space.
On each transport iteration, the hierarchy of equations is solved according to
the $V$-cycle illustrated in Figure 1. The detailed iteration scheme is
presented in Algorithm 1, where $s$ is the index of transport iterations.

Figure 1:  $V$-cycle

The high-order transport equation (26) for material $\ell$ depends on the CEA
of the angular flux of another material $\ell^{\prime}$. The LOYM equations
(7) for material $\ell$ also contain terms with the partial scalar fluxes and
factors of material $\ell^{\prime}$. In this paper, fully iterative approach
is studied. The multilevel method executes relaxation cycles at the level of
the material high-order equations using Gauss-Seidel (GS) iterations. The
scattering term is defined by the material scalar fluxes obtained from the
low-order problems. $n$ is the index of inner GS iterations. The number of GS
cycles equals $n_{max}$. The same approach is applied to solve coupled
material LOYM equations (7). At this level, only one GS iteration over the
material LOYM equations is performed. The modified LOYM equations (22) for
materials are decoupled due to the prolongation operator that defines the RHS.
There are no inner iterations at this stage.

To discretize the high-order transport equation (26) the linear discontinuous
(LD) finite element method is used. The LOYM and LOQD equations are
discretized consistently with the LD scheme for Eq. (26) [22].

$s=0$

while  _$||\big{
<}\phi\big{>}^{(s)}-\big{<}\phi\big{>}^{(s-1)}||_{\infty}>\epsilon||\big{<}\phi\big{>}^{(s)}||_{\infty}$
_ do

$s=s+1$ , $n=0$ , $\psi_{2}^{(0,s+1/3)}=\psi_{2}^{(s)}$

Level 1: high-order transport problem for materials

while  _$n\leq n_{max}$ _ do

$n=n+1$

$\mathcal{L}_{1}p_{1}\psi_{1}^{(n,s+1/3)}=|\mu|\lambda_{2}^{-1}p_{2}\psi_{2}^{(n-1,s+1/3)}+0.5\sigma_{s,1}p_{1}\phi_{1}^{(s)}+0.5p_{1}q_{1}$

$\mathcal{L}_{2}p_{2}\psi_{2}^{(n,s+1/3)}=|\mu|\lambda_{1}^{-1}p_{1}\psi_{1}^{(n,s+1/3)}\quad+0.5\sigma_{s,2}p_{2}\phi_{2}^{(s)}+0.5p_{2}q_{2}$

$\psi_{\ell}^{(s+1/3)}\leftarrow\psi_{\ell}^{(n_{max},s+1/3)}$

Compute $C_{\ell}^{\pm(s+1/3)}$ and $E_{\ell}^{\pm(s+1/3)}$ ($\ell=1,2$),
$\phi_{2}^{\pm(s+1/3)}$

Level 2: low-order problem for materials

Perform one GS material iteration to solve LOYM eqs. (7) for
$\phi_{\ell}^{\pm(s+2/3)}$ with $\phi_{2}^{\pm(s+1/3)}$ as an initial guess

Compute $\\!\big{<}\\!\\!\big{<}\beta^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+2/3)}$,
$\big{<}\\!\\!\big{<}\gamma^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+2/3)}$,
$\big{<}\\!\\!\big{<}\tilde{C}^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+2/3)}$,
$\big{<}\\!\\!\big{<}E^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+2/3)}$,
$\big{<}\\!\\!\big{<}\sigma_{a}\big{>}\\!\\!\big{>}^{(s+2/3)}$,
$\big{<}\\!\\!\big{<}\sigma_{t}\big{>}\\!\\!\big{>}^{(s+2/3)}$,
$\big{<}\\!\\!\big{<}\eta\big{>}\\!\\!\big{>}^{(s+2/3)}$,
$\big{<}\\!\\!\big{<}E\big{>}\\!\\!\big{>}^{(s+2/3)}$,
$\big{<}\\!\\!\big{<}C^{\pm}\big{>}\\!\\!\big{>}^{(s+2/3)}$

Level 3: low-order problem for total ensemble average moments

Solve LOQD eqs. (13) for $\big{<}\phi\big{>}^{(s+1)}$ and
$\big{<}J\big{>}^{(s+1)}$

Compute $\big{<}\phi^{\pm}\big{>}^{(s+2/3)}$ and
$\big{<}J^{\pm}\big{>}^{(s+2/3)}$

Level 2: low-order problem for materials

Solve LOYM eqs. (22) for $\phi_{\ell}^{\pm(s+1)}$

Compute $\\!\big{<}\\!\\!\big{<}\beta^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+1)}$,
$\big{<}\\!\\!\big{<}\gamma^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+1)}$,
$\big{<}\\!\\!\big{<}\tilde{C}^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+1)}$,
$\big{<}\\!\\!\big{<}E^{\pm}\big{>}\\!\\!\big{>}^{\\!(s+1)}$

Compute $\big{<}\phi^{\pm}\big{>}^{(s+1)}$ and $\big{<}J^{\pm}\big{>}^{(s+1)}$

Prolongation

$\phi_{\ell}^{(s+1)}\leftarrow(\phi_{\ell}^{-(s+1)}+\phi_{\ell}^{+(s+1)})\big{<}\phi\big{>}^{(s+1)}\big{/}\sum_{\ell^{\prime}=1}^{2}p_{\ell^{\prime}}\big{(}\phi_{\ell^{\prime}}^{-(s+1)}+\phi_{\ell^{\prime}}^{+(s+1)}\big{)}$

$\tilde{\phi}^{\pm(s+1)}_{\ell}\leftarrow\phi_{\ell}^{\pm(s+1)}\big{<}\phi^{\pm}\big{>}^{(s+1)}\big{/}\sum_{\ell^{\prime}=1}^{2}p_{\ell^{\prime}}\phi_{\ell^{\prime}}^{\pm(s+1)}$

$\psi_{\ell}^{(s+1)}\leftarrow\pm\psi_{\ell}^{(s+1/3)}\tilde{\phi}_{\ell}^{\pm(s+1)}\big{/}\int_{0}^{\pm
1}\psi_{\ell}^{(s+1/3)}d\mu^{\prime}$, $\mu\gtrless 0$

Algorithm 1 The multilevel iteration method for the BSM transport equation

##  4  Numerical Results

To analyze the multilevel iteration method we use four sets of numerical test
problems [2]. The data for the tests are presented in Table 1. The problems
are defined for a slab $0\leq x\leq 100$ with vacuum boundaries. The spatial
mesh is uniform with 100 cells. The double $S_{4}$ Gauss-Legendre quadrature
set is used. The sets of tests differ from each other by values of
$\lambda_{\ell}\sigma_{t,\ell}$. The problems in each set have materials with
different scattering ratios as well as total cross sections.

Table 1:  Data for Test Problems Test | $\sigma_{t,1}$ | $\frac{\sigma_{s,1}}{\sigma_{t,1}}$ | $\lambda_{1}$ | $q_{1}$ | $\sigma_{t,2}$ | $\frac{\sigma_{s,2}}{\sigma_{t,2}}$ | $\lambda_{2}$ | $q_{1}$ | $\lambda_{1}\sigma_{t,1}$ | $\lambda_{2}\sigma_{t,2}$
---|---|---|---|---|---|---|---|---|---|---
A1 | $\frac{10}{99}$ | 0 | $\frac{99}{100}$ | 1 | $\frac{100}{11}$ | 1 | $\frac{11}{100}$ | 1 | 0.1 | 1
A2 | $\frac{10}{99}$ | 1 | $\frac{99}{100}$ | 1 | $\frac{100}{11}$ | 0 | $\frac{11}{100}$ | 1 | 0.1 | 1
A3 | $\frac{10}{99}$ | 0.9 | $\frac{99}{100}$ | 1 | $\frac{100}{11}$ | 0.9 | $\frac{11}{100}$ | 1 | 0.1 | 1
B1 | $\frac{10}{99}$ | 0 | $\frac{99}{10}$ | 1 | $\frac{100}{11}$ | 1 | $\frac{11}{10}$ | 1 | 1 | 10
B2 | $\frac{10}{99}$ | 1 | $\frac{99}{10}$ | 1 | $\frac{100}{11}$ | 0 | $\frac{11}{10}$ | 1 | 1 | 10
B3 | $\frac{10}{99}$ | 0.9 | $\frac{99}{10}$ | 1 | $\frac{100}{11}$ | 0.9 | $\frac{11}{10}$ | 1 | 1 | 10
C1 | $\frac{2}{101}$ | 0 | $\frac{101}{20}$ | 1 | $\frac{200}{101}$ | 1 | $\frac{101}{20}$ | 1 | 0.1 | 10
C2 | $\frac{2}{101}$ | 1 | $\frac{101}{20}$ | 1 | $\frac{200}{101}$ | 0 | $\frac{101}{20}$ | 1 | 0.1 | 10
C3 | $\frac{2}{101}$ | 0.9 | $\frac{101}{20}$ | 1 | $\frac{200}{101}$ | 0.9 | $\frac{101}{20}$ | 1 | 0.1 | 10
D1 | $\frac{1}{99}$ | 0 | $\frac{99}{100}$ | 1 | $\frac{10}{11}$ | 1 | $\frac{11}{10}$ | 1 | 0.01 | 1
D2 | $\frac{1}{99}$ | 1 | $\frac{99}{100}$ | 1 | $\frac{10}{11}$ | 0 | $\frac{11}{10}$ | 1 | 0.01 | 1
D3 | $\frac{1}{99}$ | 0.9 | $\frac{99}{100}$ | 1 | $\frac{10}{11}$ | 0.9 | $\frac{11}{10}$ | 1 | 0.01 | 1
Table 2:  Numerically Estimated Spectral Radii $n_{max}$ | A1 | A2 | A3 | B1 | B2 | B3 | C1 | C2 | C3 | D1 | D2 | D3
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 0.4 | 0.41 | 0.41 | 0.23 | 0.12 | 0.12 | 0.11 | 0.06 | 0.10 | 0.46 | 0.43 | 0.46
2 | 0.16 | 0.17 | 0.16 | 0.25 | 0.12 | 0.1 | 0.14 | 0.03 | 0.11 | 0.19 | 0.20 | 0.21

Table 2 presents numerically estimated spectral radii in every test. The
results for algorithms with one and two GS inner iterations in the high-order
transport problem ($n_{max}=1,2$) are shown. To demonstrate details of
iterative behavior of the algorithms, Figures 3 and 5 show convergence rates
of the total ensemble average scalar flux

$\rho^{(s)}=\frac{\big{|}\big{|}\Delta\big{<}\phi\big{>}^{(s)}\big{|}\big{|}_{\infty}}{\big{|}\big{|}\Delta\big{<}\phi\big{>}^{(s-1)}\big{|}\big{|}_{\infty}}\,,\quad\Delta\big{<}\phi\big{>}^{(s)}=\big{<}\phi\big{>}^{(s)}-\big{<}\phi\big{>}^{(s-1)}\,$
(28)

in Tests $A_{k}$ and $C_{k}$, respectively. The algorithm with two relaxation
cycles in the high-order problem ($n_{max}=2$) converges fast in all tests. In
Tests $A_{k}$ and $D_{k}$, the spectral radii of this algorithm are smaller
than those of the algorithm with one relaxation cycle ($n_{max}=1$). In Tests
$B_{k}$ and $C_{k}$, the algorithm with one relaxation cycle is rapidly
converging. Performance of the algorithms with 1 and 2 relaxation cycles are
close to each other in these two sets of tests.

Figures 3 and 5 present convergence curves of
$||\Delta\big{<}\phi\big{>}^{(s)}||_{\infty}$ and
$||\Delta\phi_{\ell}^{(s)}||_{\infty}$ in Tests $A_{k}$ and $C_{k}$, where
$\Delta\phi_{\ell}^{(s)}=\phi_{\ell}^{(s)}-\phi_{\ell}^{(s-1)}$. The CEA of
the material scalar fluxes and total ensemble average of the scalar flux
converge with close rates. The convergence behavior in other tests is similar.

(a) Test A1

(b) Test A2

(c) Test A3

Figure 2:  Convergence rates $\rho^{(s)}$ in Tests Ak

(a) Test A1

(b) Test A2

(c) Test A3

Figure 3:  Convergence curves of $\big{<}\phi\big{>}^{(s)}$ and
$\phi_{\ell}^{(s)}$ in Tests Ak, $n_{max}=2$

(a) Test C1

(b) Test C2

(c) Test C3

Figure 4:  Convergence rates $\rho^{(s)}$ in Tests Ck

(a) Test C1

(b) Test C2

(c) Test C3

Figure 5:  Convergence curves of $\big{<}\phi\big{>}^{(s)}$ and
$\phi_{\ell}^{(s)}$ in Tests Ck, $n_{max}=2$

##  5  Conclusions

A new iteration method for the transport equation in binary stochastic media
has been presented. It is formulated by means of a hierarchy of equations
consisting of the high-order transport equation for the CEA of the angular
flux and the low-order equations for (i) the CEA of the material partial
scalar fluxes and (ii) the total ensemble average scalar flux and current. The
obtained results show that the multilevel method accelerates transport
iterations in the considered class of test problems. Further analysis is
needed to study efficiency of the iteration method. One of items is
convergence in the atomic mix limit as $\lambda_{\ell}\sigma_{t,\ell}\to 0$.

There are various possible modifications that can be applied to the multilevel
method, for example, Anderson acceleration, nonlinear Krylov acceleration and
advanced prolongation operators for the CEA of the angular flux. The low-order
equations of the multilevel method correspond to different scales of the
physical problem. This method has a potential for nonlinear radiation
transport problems in random media in which the transport equation is coupled
with multiphysics equations.

## References

  * [1] C. D. Levermore, G. C. Pomraning, D. L. Sanzo, J. Wong, Linear transport theory in a random medium, Journal of Mathematical Physics 27 (1986) 2526–2536.
  * [2] M. L. Adams, E. W. Larsen, G. C. Pomraning, Benchmark results for particle transport in a binary Markov statistical medium, Journal of Quantitative Spectroscopy Radiative Transfer 4 (1989) 253–266.
  * [3] F. Malvagi, G. C. Pomraning, Renormalized equations for linear transport in stochastic media, J. Mathematical Physics 31 (1990) 892–900.
  * [4] G. Pomraning, Linear Kinetic Theory and Particle Transport in Stochastic Mixtures, World Scientific, Singapore, New York, London, Hong Kong, 1991.
  * [5] B. Ching, T. Palmer, Analysis of iterative schemes for binary stochastic mixture transport equations, Transactions of the American Nuclear Society 82 (2000) 148.
  * [6] D. S. Miller, G. Rodrigue, F. Graziani, Convergence analysis of source iteration for the solution of coupled transport equations on a stochastic media, Transport Theory and Statistical Physics 31 (2002) 65–83.
  * [7] G.B. Zimmerman, M. L. Adams, Algorithms for Monte Carlo particle transport in binary statistical mixtures, Transactions of the American Nuclear Society 63 (1991) 287–288.
  * [8] A. K. Prinja, E. D. Fichtl, Atomic mix synthetic acceleration of dose computations in binary statistical media, Nuclear Science and Engineering 155 (2007) 441–448.
  * [9] T. S. Palmer, B. S. Ching, A “two-grid’” acceleration of binary stochastic mixture deterministic transport iterations in slab geometry, Annals of Nuclear Energy 35 (2008) 68–74.
  * [10] T. S. Palmer, A coupled diffusion synthetic acceleration for binary stochastic mixture transport iterations, Nuclear Science and Engineering 158 (2008) 40–48.
  * [11] E. D. Fichtl, J. S. Warsa, A. K. Prinja, Krylov iterative methods and synthetic acceleration for transport in binary statistical media, Journal of Computational Physics 228 (2009) 8413–8426.
  * [12] V. Y. Gol’din, On mathematical modeling of problems of non-equilibrium transfer in physical systems, in: Modern Problems of Mathematical Physics and Computational Mathematics, Nauka, Moscow, 1982, pp. 113–127, in Russian.
  * [13] D. Y. Anistratov, V. Y. Gol’din, Nonlinear methods for solving particle transport problems, Transport Theory and Statistical Physics 22 (1993) 42–77.
  * [14] D. S. Miller, G. Rodrigue, F. Graziani, Benchmarks and models for time-dependent grey radiation transport with material temperature in binary stochastic media, J. of Quantitative Spectroscopy & Radiative Transfer 70 (2001) 115–128.
  * [15] D. Y. Anistratov, Stability analysis of a multilevel quasidiffusion method for thermal radiative transfer problems, J. of Computational Physics 376 (2019) 186–209.
  * [16] J. J. Yvon, J. Nucl Energy 4 (1957) 305.
  * [17] R. Mertens, Simon Stevin Supplement 30.
  * [18] G. Bell, S. Glasstone, Nuclear Reactor Theory, Van Nostrand Reinhold Company, New York,Cincinnati, Toronto, London, Melbourne, 1970.
  * [19] S. Santandrea, R. Sanchez, Analysis and improvements of the DPN acceleration technique for the method of characteristics in unstructured meshes, Annals of Nuclear Energy 32 (2005) 163–193.
  * [20] S. S. Nikolaishvilli, On solution of one-velocity transport equation based on Yvon-Mertens approximation, Soviet Atomic Energy 20 (1966) 344.
  * [21] V. Y. Gol’din, A quasi-diffusion method of solving the kinetic equation, USSR Computational Mathematics and Mathematical Physics. 4 (1964) 136–149.
  * [22] D. Y. Anistratov, J. S. Warsa, Discontinuous finite element quasidiffusion methods, Nuclear Science and Engineering 191 (2018) 105–120.