  *

# Plasmoid drift and first wall heat deposition during ITER H-mode dual-SPIs
in JOREK simulations

D. Hu1,∗, F. J. Artola2, E. Nardon3, M. Lehnen2, M. Kong4, D. Bonfiglio5, M.
Hoelzl6, G.T.A. Huijsmans3,7 & JOREK Team111See Hoelzl et al 2021
(https://doi.org/10.1088/1741-4326/abf99f) for the JOREK Team. 1Beihang
University, No. 37 Xueyuan Road, Haidian District, 100191 Beijing, China.
2ITER Organization, Route de Vinon sur Verdon, CS 90 046,13067 Saint Paul-lez-
Durance, Cedex, France.  3CEA, IRFM, F-13108 Saint-Paul-Lez-Durance, France
4École Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC),
CH-1015 Lausanne, Switzerland  5Consorzio RFX-CNR, ENEA, INFN, Università di
Padova, Acciaierie Venete SpA. I-35127 Padova, Italy.  6Max Planck Institute
for Plasma Physics, Boltzmannstr. 2, 85748 Garching b. M., Germany  7Eindhoven
University of Technology, De Rondom 70 5612 AP Eindhoven, the Netherlands.  *
<EMAIL_ADDRESS>
###### Abstract

The heat flux mitigation during the Thermal Quench (TQ) by the Shattered
Pellet Injection (SPI) is one of the major elements of disruption mitigation
strategy for ITER. It’s efficiency greatly depends on the SPI and the target
plasma parameters, and is ultimately characterised by the heat deposition on
to the Plasma Facing Components (PFCs). To investigate such heat deposition,
JOREK simulations of neon-mixed dual-SPIs into ITER baseline H-mode and a
“degraded H-mode” with and without good injector synchronization are performed
with focus on the first wall heat flux and its energy impact. It is found that
low neon fraction SPIs into the baseline H-mode plasmas exhibit strong major
radial plasmoid drift as the fragments arrive at the pedestal, accompanied by
edge stochasticity. Significant density expulsion and outgoing heat flux
occurs as a result, reducing the mitigation efficiency. Such drift motion
could be mitigated by injecting higher neon fraction pellets’, or by
considering the pre-disruption confinement degradation, thus improving the
radiation fraction. The radiation heat flux is found to peak in the vicinity
of the fragment injection location in the early injection phase, while it
relaxes later on due to parallel impurity transport. The overall radiation
asymmetry could be significantly mitigated by good synchronization. Time
integration of the local heat flux is carried out to provide its energy impact
for wall heat damage assessment. For the baseline H-mode case with full pellet
injection, melting of the stainless steel armour of the diagnostic port could
occur near the injection port, which is acceptable, without any melting of the
first wall tungsten tiles. For the degraded H-mode cases with quarter-pellet
SPIs, which have $1/4$ total volume of a full pellet, the maximum energy
impact approaches the tolerable limit of the stainless steel with un-
synchronized SPIs, and stays well below such limit for the perfectly
synchronized ones.

##  1 Introduction

The Thermal Quench (TQ) heat flux mitigation by impurity Shattered Pellet
Injection (SPI) is an important part of the ITER Disruption Mitigation System
(DMS) strategy, the aim of which is to evenly redistribute the thermal energy
stored within the plasma onto the Plasma Facing Components (PFCs), so that
localized heat deposition on the divertor or the first wall and the
consequential material damage could be avoided during the TQ process [1]. The
TQ mitigation efficiency thus greatly depends on both the radiated energy
fraction, as well as the asymmetry of the radiation heat flux onto the first
wall.

There have been intensive studies, both experimentally [2, 3, 4, 5, 6] and
numerically [7, 8, 9, 10], on both the total radiated energy and the radiation
power density within the plasma volume after SPIs. It has been found in
general that an apparent toroidal asymmetry in the radiation power density
exists for the single SPI cases, especially in the early injection phase. The
radiation power density profile tends to relax as the TQ proceeds, and the
Toroidal Peaking Factor (TPF) of the radiation power density at the time of
the maximum radiation power tends to be smaller compared with the maximum TPF
during the whole TQ phase in both experimental observations [3] and numerical
simulations [8].

Although the aforementioned studies already provide valuable insight into the
asymmetric TQ mitigation process, the ultimate criterion of the disruption
mitigation efficiency should be the heat flux onto the first wall and its
accumulated energy impact, as the latter is directly related to the wall
material temperature rise and thus the wall damage during disruption[11]. The
recent change of the ITER first wall material to tungsten offers more
tolerance to such transient heat flux [12], however, should the accumulated
energy impact exceed a certain threshold, even tungsten could experience
significant melting [11]. Meanwhile, the stainless steel cover for the port
window suffers from a lower energy impact threshold which is about $1/3$ of
the tungsten threshold, thus is more susceptible of melting, although such
melting is acceptable as long as the melt threshold is not dramatically
exceeded for the stainless steel. [13, 14]. Furthermore, even without
exceeding the tolerable maximum energy impact, transient heat fluxes could
already produce structural modification on tungsten surface which could affect
the mechanical and thermal properties of the first wall [15, 16, 17, 18, 19].
Hence it is desirable to obtain the heat flux and energy impact distribution
onto the first wall after SPIs to assess the efficiency of the TQ mitigation
and to validate the ITER PFC heat tolerance specifications.

To take a first look into the aforementioned heat deposition behaviour in ITER
H-mode plasmas, we carry out 3D non-linear reduced MHD simulations using the
JOREK code [20] and its collisional-radiative impurity module [21]. Two sets
of equilibria are considered, a “baseline” H-mode and a “degraded” one. The
latter means to take into account the natural confinement degradation in the
disruption precursor phase, which significantly reduces the thermal energy
before the triggering of the DMS [22]. Two sets of dual-SPI configurations are
also considered, one “full pellet” with small neon mixture ratio and one
“quarter-pellet” with the same amount of neon atoms but one fourth the
hydrogen content. For the full pellet dual-SPI into baseline H-mode case, a
strong plasmoid drift along the major radial direction is found to occur due
to strong local over-pressure within the newly ablated plasmoid cloud which
induces polarization and subsequent ${\bf E}\times{\bf B}$ drift towards the
Low Field Side (LFS) [23]. Such drift and its accompanying MHD response would
undermine the assimilation when fragments are injected from the LFS mid-plane,
as well as cause strong boundary heat flux as the edge becomes stochastic.
These detrimental behaviours could be mitigated if the higher neon mixture
ratio quarter-pellet is used so that the local over-pressure is suppressed by
local radiation [24] or by considering the degraded H-mode. Indeed, the
quarter-pellet SPIs into degraded H-mode shows good impurity assimilation and
radiation fraction. Moreover, the radiation heat flux onto the first wall is
obtained by the Raysect/CHERAB code suite integrated within Integrated
Modelling & Analysis Suite (IMAS) [25, 26] and their energy impact is
calculated. It is found that the maximum energy impact for the full pellet
dual-SPI into the baseline H-mode case exceeds the tolerable limit of the
stainless steel near the injector, but stays well below that of the tungsten.
For the quarter-pellet dual-SPI into the degraded H-mode, the maximum energy
impact only exceeds the stainless steel limit slightly even when there is
$1ms$ delay between the injectors, while it stays well below said limit if the
synchronization is perfect.

The rest of the paper is arranged as the following: In Section 2, the two
target equilibria are introduced, and our basic assumptions, simulation setup
as well as the SPI configurations and parameters are presented. We also
introduce our method of calculating the accumulated energy impact from a time-
varying heat flux distribution onto the first wall. In Section 3, the plasmoid
drift behaviour during H-mode SPIs and its mitigation are investigated, and
the possible consequences of such drift are discussed. The radiation heat flux
and its energy impact onto the first wall for both the baseline and degraded
H-mode cases are shown in Section 4, and their implication to the PFC damage
is discussed. Finally, conclusion and further discussions on the TQ mitigation
efficiency with neon SPIs are presented in Section 5.

##  2 Simulation setup and calculation of the energy impact

###  2.1 Simulation setup and basic assumptions

The 3D non-linear reduced MHD version of JOREK with collisional-radiative
impurity model [20, 21] is used for this study, and our governing equations
are the same as that described in Ref. [9] without modifications, so we shall
not repeat them here. The target equilibria are one baseline H-mode chosen
from the ITER reference scenarios [27], and one degraded H-mode constructed
from the baseline H-mode equilibrium by increasing the perpendicular heat
conduction and diffusion artificially thus producing a H-L back transition,
mimicking the naturally occurring confinement degradation before the onset of
the TQ [22]. Both cases have the same toroidal magnetic field at $B_{T}=5.3T$
and the total plasma current at $I_{p}=15MA$. The initial equilibrium
profiles, including the electron temperature, density, current density and the
safety factor $q$ profiles for both target plasmas are shown in Fig. 1. The
initial thermal energy is $190MJ$ for the degraded H-mode case and $370MJ$ for
the baseline H-mode case.

|
---|---
(a) | (b)
Figure 1: The initial equilibria of (a) the baseline H-mode and (b) the
degraded H-mode, $\bar{\psi}$ is the normalized poloidal flux.
---
Figure 2: The simulation grid. The boundary is assumed to be ideally
conducting.

The grid we use is presented in Fig. 2, where the boundary approximately
represents the position of the thin wall equivalent of the ITER blanket
modules [28]. We assume the boundary to be ideally conducting in this study.
We include $12$ toroidal harmonics from $n=0$ to $n=11$ in our simulations,
and our ablative density source is artificially elongated along the toroidal
direction due to limitation in the toroidal resolution. This artificial
toroidal spread with $\Delta\phi_{NG}=0.3$ radian, corresponding to $2.46m$ to
$1.86m$ as the fragments travel from the edge to the core, enable us to keep
the Fourier components up to the first order in strength with the following
density deposition shape:

$\displaystyle
S_{n}\propto\exp{\left(-\frac{\left(R-R_{f}\right)^{2}+\left(Z-Z_{f}\right)^{2}}{\Delta
r^{2}_{NG}}\right)}\times\exp{\left(-\left(\frac{\phi-\phi_{f}}{\Delta\phi_{NG}}\right)^{2}\right)}.$
(1)

Here $R_{f}$, $Z_{f}$ and $\phi_{f}$ are the coordinates of the ablating
fragment. We choose $\Delta r_{NG}=0.08m$ in this study, but wish to emphasize
that the overlapping between the ablation clouds from different fragments
tends to create an envelop of plasmoid, so that our result is insensitive to
the exact choice of $\Delta r_{NG}$ so long as it is no larger than the plume
width. The exact value of the density source is determined by the Neutral Gas
Shielding (NGS) ablation model [29], the detailed implementation is described
in Ref. [9]. The density generated by ablation would then relax by parallel
convection driven by the pressure gradient along the field lines, following
the continuity equation.

The two temperature model described in Ref. [9] is used, where we track the
evolution of the electron and the ion temperatures separately. The Braginskii
parallel heat conduction [30] is used for each species respectively. We impose
an upper limit of this parallel heat conduction corresponding to the value of
the Braginskii electron thermal conductivity at $T_{e}\simeq 1.6keV$, beyond
which the heat conductivity would not rise any more. This is an attempt to
mimic the flux limit of the free streaming electrons and to prevent the
overestimation of Braginskii thermal conduction in the high temperature
regime. With the initial density, the perpendicular particle diffusion is set
to be constant at $6m^{2}/s$ and the perpendicular heat conduction at
$10m^{2}/s$. The choice of these coefficients are meant approach the gyro-Bohm
scaling [31] while still maintain numerical stability. We do not expect the
exact value to impact our result drastically since the stochastic transport
and the macroscopic convective transport dominate over the turbulent one after
the injection for our simulations. The perpendicular viscosity is set to
$30m^{2}/s$ and the parallel one to $150m^{2}/s$, these values are chosen to
meet the requirement of numerical instability suppression. The Spitzer
resistivity [32] with the effective charge contribution [33] included is used
for our simulation. Here, we have neglected the neo-classical resistivity
contribution from the trapped electrons as the temperature drop and the
density rise after the injection would prevent the electrons to complete a
bounce period before collision. This is consistent with the previous finding
that the mean-free-path of several hundred $eV$ electrons is comparable with
the size of the cold, dense plasmoid around the fragments after ITER SPIs [9].
We assume the plasma is transparent to the impurity radiation, hence our
result might over-estimate the peak radiation heat flux due to neglecting the
plasma opacity.

The SPIs are carried out from the outer equatorial ports EQ-08 and EQ-17,
which locate at $R=8.4m$, $Z=0.685m$ and opposite toroidal positions, along
the major radial direction. The vertex angle of the injection cone is $20$
degrees. The fragment size distribution follows that of the Statistical
Fragmentation model [34, 35]. Our two sets of SPI parameters are the “$1\%$
neon full pellet SPI” dual-SPIs with $2.5\times 10^{22}$ neon atoms and
$1.8\times 10^{24}$ hydrogen atoms for each injector, as well as the “$5\%$
neon quarter pellet SPI” dual-SPIs with $2.5\times 10^{22}$ neon atoms and
$4.5\times 10^{23}$ hydrogen atoms for each injector. For the full pellet
case, the pellet is shattered into 300 fragments, while for the quarter pellet
case the fragment number is 100. Their characteristic fragment size is
$\kappa_{p}^{-1}\simeq 1.43mm$ and $\kappa_{p}^{-1}\simeq 1.30mm$
respectively. Both cases use $500m/s$ reference injection velocity with $\pm
40\%$ velocity spread. We assume the fragment velocity distributes uniformly
over the spreading range, although in reality the velocity distribution has a
correlation with the fragment size [36]. We do not expect this to strongly
alter the dynamics we investigate here, since one could usually modify the
shattering angle to partly cancel this effect. It should also be noted that,
for our baseline H-mode case, we enabled the velocity spread only in the EQ-08
plume and set all fragments in the EQ-17 plume with the same velocity value,
albeit with correct direction spread so that there exists an asymmetry in the
dual-SPI despite the two plumes being injected at the same time. We
nonetheless don’t expect this to greatly impact our results, since in reality
there would always exist some asymmetry between the plumes even with very good
synchronization between the injectors.

Notation | Equilibria | Neon | Hydrogen | Frag. | Delay
---|---|---|---|---|---
BH-FP-dt0 (case 1) | baseline | $2\times 2.5\times 10^{22}$ | $2\times 1.8\times 10^{24}$ | $300$ |  $0ms$ (Asymm.)
DH-FP-dt0 (case 2) | degraded | $2\times 2.5\times 10^{22}$ | $2\times 1.8\times 10^{24}$ | $300$ | $0ms$
DH-QP-dt0 (case 3) | degraded | $2\times 2.5\times 10^{22}$ | $2\times 4.5\times 10^{23}$ | $100$ | $0ms$
DH-QP-dt1 (case 4) | degraded | $2\times 2.5\times 10^{22}$ | $2\times 4.5\times 10^{23}$ | $100$ | $1ms$
Table 1: The injection parameters for the SPI considered in this study. Note
that for BH-FP-dt0 there exists an asymmetry between the plumes although they
are injected at the same time.

We have summarized the cases we investigate in this paper in Table 1, where
the target equilibria, the neon and hydrogen injection amount, the fragment
number and the delay between the injectors are shown. For all cases, the first
fragment of the EQ-08 plume arrives on the LCFS approximately at $t=0.38ms$
and the last one arrives approximately at $0.89ms$. All the fragments from the
EQ-17 plume in the BH-FP-dt0 case arrive approximately at $t=0.54ms$ since
there is no velocity spread. The arrival time of the fragments from the EQ-17
plume for the rest of the cases is the same with that of the EQ-08 plume, plus
the delay between the injectors shown in Table 1. We will mainly investigate
BH-FP-dt0, DH-FP-dt0 & DH-QP-dt0 in Section 3, and BH-FP-dt0, DH-QP-dt0 & DH-
QP-dt1 in Section 4.

###  2.2 The energy impact

To estimate the wall material temperature rise in an inertially cooled
scenario, it is important to calculate the accumulated heat on the wall
surface considering both the incoming heat flux and the conductive cooling of
the PFCs. For the simple case of 1D heat conduction concerning a semi-infinite
slab and constant heat flux $q_{s}$, the temperature increase $\Delta T$ after
$\Delta t$ time of heat exposure is [11]

$\displaystyle\Delta T(\Delta t)=\frac{2}{\sqrt{\pi}b}\Delta Q(\Delta
t)=\frac{2}{\sqrt{\pi}b}q_{s}\sqrt{\Delta t}.$ (2)

Here $b=\sqrt{\kappa\rho c}$, $\rho$ is the material density, $c$ is its
specific heat capacity, $\kappa$ is its heat conductivity. The temperature
rise is proportional, via the material-dependent factor $b$, to the so-called
energy impact $\Delta Q\equiv q_{s}\sqrt{\Delta t}$ which itself is material-
independent. The maximum tolerable energy impact could then be obtained for
each individual material [11, 14], thus providing insights into the
probability of severe material damage for a given period of heat deposition.

When considering a time varying heat flux which is more relevant during
transient TQ process, one could use the fundamental solution of the heat
equation to obtain an analytical description of the heat accumulation on the
slab surface without numerically solving the 1D heat conduction equation. This
is equivalent to considering the heat pulse $q_{s}(t)$ as a series of
continuously arriving “heat packets”, each of which individually undergoes a
1D random walk towards infinity with step-length determined by the material
heat conductivity after their arrival. The total heat distribution after a
certain time period $\Delta t$ is then the summation of the 1D spatial
probability distribution of all the packets. For each such heat packet, the
fundamental solution on $x\in\left[0,+\infty\right)$ is

$\displaystyle\Phi\left(x,\Delta t\right)=\frac{1}{\sqrt{\pi\alpha\Delta
t}}\exp{\left(-\frac{x^{2}}{4\alpha\Delta t}\right)}.$ (3)

Here $\alpha$ is the thermal diffusivity and $x$ is the distance along the 1D
direction. Since we are only concerned with the temperature rise at the wall
surface, $x$ could be practically set to zero, and we obtain the time decay of
each individual heat packet to be proportional to $\Delta t^{-1/2}$. The
energy impact integrating over all the heat packets could then be obtained by
the following convolution:

$\displaystyle\Delta
Q(t)=\frac{1}{2}\int_{t_{0}}^{t}{\frac{q_{s}(t^{\prime})}{\sqrt{t-t^{\prime}}}dt^{\prime}}.$
(4)

Here $t_{0}$ is the beginning time of the heat pulse. For a constant $q_{s}$,
one naturally recovers the previous result $\Delta Q(\Delta
t)=q_{s}\sqrt{\Delta t}$. Henceforth we will calculate the energy impact using
Eq. (4) and estimate the potential heat damage based on it.

##  3 Plasmoid drift, accompanied heat flux and its mitigation

For the baseline H-mode BH-FP-dt0 studied here, it is found that the high
thermal energy reservoir could result in strong local pressure peak around the
vanguard fragments, which induces vertical polarization, which in turn results
in ${\bf E}\times{\bf B}$ drift along the major radial direction [23]. In our
simulations, the polarization could be implicitly obtained through the
evolution of the vorticity equation in Ref. [9]:

$\displaystyle R\nabla\cdot\left[R^{2}\frac{\partial}{\partial
t}\left(\rho\nabla_{pol}u\right)\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left\\{R^{2}\left|\nabla_{pol}u\right|^{2},R^{2}\rho\right\\}+\left\\{R^{4}\rho\omega,u\right\\}$
(5)
$\displaystyle-R\nabla\cdot\left[R^{2}\nabla_{pol}u\nabla\cdot\left(\rho{\bf
v}\right)\right]+\left\\{\psi,j\right\\}$
$\displaystyle-\frac{F_{0}}{R}\frac{\partial
j}{\partial\phi}+\left\\{P,R^{2}\right\\}+R^{3}\mu_{\bot}\left(T_{e}\right)\nabla_{pol}^{2}\omega.$

Here $\rho$ is the mass density, $u$ is the stream function ${\bf
v}=v_{\|}{\bf B}-R^{2}\nabla u\times\nabla\phi$. To see this, one realizes
that the stream function $u$ is associated with the electric potential, and
the second term on the RHS of the velocity definition is essentially the ${\bf
E}\times{\bf B}$ drift. Thus the Poisson bracket $\left\\{P,R^{2}\right\\}$
term in Eq.,(5), which includes the $\nabla B$ drift contribution, would drive
the aforementioned polarization in the presence of the local plasmoid over-
pressure. Such polarization would then in turn result in the drift motion. The
above discussed strong drift motion is found when the fragments arrive on the
pedestal of the baseline H-mode. It should be noted that due to the asymmetry
in our two plumes as mentioned in Section 2.1, the EQ-08 plume arrives on the
pedestal first and triggers the strong drift. The resulting drift motion
expels significant amounts of injected materials out of the Last Closed Flux
Surface (LCFS), undermining the injection assimilation, as can be seen in Fig.
3, which shows to the same toroidal location as the EQ-08 plume. The local
electron over-pressure could be readily identified at $t=0.41ms$ in Fig. 3(a)
as the fragments penetrate the equilibrium LCFS represented by the white
contour. This over-pressure region essentially is the envelop of several
plasmoid produced by individual vanguard fragments as they enter the pedestal.
This high pressure plasmoid envelop is then seen convected outward in Fig.
3(b) due to the major radial drift and finally relaxes outside of the
separatrix in Fig. 3(c).

|  |
---|---|---
(a) | (b) | (c)
Figure 3: The electron pressure profile at the same toroidal location as the EQ-08 plume for BH-FP-dt0 at (a) $t=0.41ms$, (b) $t=0.42ms$ and (c) $t=0.46ms$. The white contour indicates the position of the equilibrium LCFS. The over-pressure as a result of the vanguard fragments arriving on the pedestal as well as its fast outward drift and expulsion can be identified. |
---|---
(a) | (b)
Figure 4: The mid-plane cut of the $n=0$ electron density for (a) BH-FP-dt0
and (b) DH-QP-dt0 at various times. The high pedestal temperature and the
consequential plasmoid drift in BH-FP-dt0 results in significant density
deposition outside of the pedestal in the early injection phase as is shown by
the black, blue and green lines in (a). The LCFS on the low field side
approximately locates at $R=8.2m$.

Consequentially, the injection assimilation within the Last Close Flux Surface
(LCFS) suffers in the early injection phase as is shown in Fig. 4, where the
comparison between the mid-plane cut of the $n=0$ electron density profile for
BH-FP-dt0 and DH-QP-dt0 at various times is shown. The LCFS on the low field
side approximately locates at $R=8.2m$ for both cases. In Fig. 4(a), comparing
the black, blue and green lines, it can be seen that the density peak
initially occurs near the LCFS, then goes outward into the scrape-off layer.
This drift motion, coupled with the associated outgoing heat flux, causes
significant density deposition in the open field line region, which is then
lost over time. As time goes on however, the fragments ultimately penetrate
deeper into the plasma and result in core density rise in the late injection
phase, as shown by the cyan and red lines. Note that the density rise in the
axis region is due to inward MHD transport, as at that time the majority of
fragments has not reached the axis yet. This behaviour is compared with DH-QP-
dt0 shown in Fig. 4(b), where no apparent density expulsion is seen and the
density peak gradually moves inward.

Such drift motion is accompanied by strong field line stochasticity in the
edge region, resulting in a stochastic parallel heat flux, as well as
perpendicular convective heat flux due to the drift transport of the high
pressure plasmoid. Poincaré plots of the magnetic field lines at $t=0.41ms$
and $t=0.46ms$ for BH-FP-dt0, as well as the simulation boundary and a surface
$10cm$ inside of that boundary which corresponds to the approximate position
of the first wall, are shown in Fig. 5, where the field lines that remain well
confined after $2000$ toroidal turns are marked by red, while those that are
lost are coloured according to the distance they travel before they hit the
simulation boundary. The figure corresponds to the same toroidal location as
the EQ-08 plume. The simulation boundary is shown in purple, and the first
wall approximation is shown in green. In Fig. 5(a), the edge region begins to
become stochastic, and the “connection length” of the field lines initialized
in the edge region to the wall decreases. Another notable feature is the
outward deformation of the closed flux surfaces close to the plasmoid
position, likely a result of the outward plasma motion. In Fig. 5(b), already
significant region around the previous pedestal has become stochastic with
short connection length to the wall, and the H-mode pedestal is collapsed.
Comparing Fig. 4(a) and Fig. 5(b), it can be seen at $t=0.46ms$ there is no
significant density assimilation within the pedestal, while strong edge
stochasticity and thus thermal energy loss already occurs.

|
---|---
(a) | (b)
Figure 5: Poincaré plots of the magnetic field lines at time (a) $t=0.41ms$
and (b) $t=0.46ms$ for BH-FP-dt0. The field lines remain well confined after
$2000$ toroidal turns are marked by red, while those that are lost are
coloured according to their “loss length”. The outward deformation of the flux
surfaces on the mid-plane in (a) is likely due to the outward plasma motion.

As a result of the above-mentioned stochastic and convective transport, the
conductive and convective heat flux onto the first wall could be significant.
In Fig. 6, the normal component of the combined heat flux normal to the first
wall approximation represented by the green line in Fig. 5 is shown for times
$t=0.42ms$ and $t=0.43ms$. Comparing Fig. 3 and Fig. 6, it can be seen that as
the over-pressured plasmoid drifts towards the wall at $t=0.42ms$, a strong
helical boundary heat flux structure begins to appear at the EQ-08 injection
location centred around $\phi=-2.57$. It successively relaxes along field
lines and the peak heat flux locations move away from the plasmoid drift
location at $t=0.43ms$. During this transient heat pulse the perpendicular
heat flux dominates over the parallel one.

|
---|---
(a) | (b)
Figure 6: The normal component of the combined conductive and convective heat flux for BH-FP-dt0 at time (a) $t=0.42ms$ and (b) $t=0.43ms$ across a surface $10cm$ within the simulation boundary, which represents the approximate location of the thin wall equivalent of the blanket module. The horizontal axis is the toroidal angle and the vertical one is the poloidal angle. The EQ-08 port locates approximately at $\phi=-2.57$. During this transient heat pulse the perpendicular heat flux dominates over the parallel one. |
---|---
(a) | (b)
Figure 7: Comparison of the total thermal energy loss with the energy consumed
by radiation and ionization for (a) BH-FP-dt0 and (b) DH-QP-dt0.

Both the aforementioned density expulsion and the accompanying transport event
are undesirable for TQ mitigation. The H-mode pedestal collapses after the
transport event is triggered and over $100MJ$ of thermal energy is prematurely
lost before the fragments penetrate deep into the core, compared with the
$370MJ$ initial thermal energy. However, the plasma usually would already
experience significant confinement degradation due to the growth of precursor
modes before the disruption onset [22] and the triggering of the SPIs. Hence
it is interesting to see if the reduced thermal energy reservoir of such a
“degraded H-mode” could still trigger such strong drift motion, and whether an
increased neon mixture ratio could help to suppress the plasmoid drift by
radiatively reducing the local over-pressure, and thus mitigate the
accompanying outgoing heat flux. Indeed, if we compare the thermal energy loss
and the radiated power of BH-FP-dt0 and DH-QP-dt0 as is shown in Fig. 7, DH-
QP-dt0 enjoys much better radiation fraction behaviour due to the mitigation
of this drift motion. Here the total loss of the thermal energy is represented
by the black solid lines, and the total radiated energy is represented by the
red solid lines. For BH-FP-dt0, the sudden rise of the discrepancy between the
thermal energy loss and the radiated energy beginning at $t=0.45ms$ is due to
the aforementioned heat flux accompanying the drift motion. This heat flux is
absent in DH-QP-dt0 as will be shown below. For both cases, the ohmic heating
is negligible until the end of the simulation when it begins to pick up due to
the thermal energy reservoir getting depleted.

|
---|---
(a) | (b)
Figure 8: The MHD spectrum of (a) DH-FP-dt0 and (b) DH-QP-dt0. The solid lines
represent the perturbed magnetic energy and the dashed lines represent that of
the kinetic energy. The red chained lines corresponds to the MHD peaks at the
time of which we would investigate the boundary heat flux distribution.

The perturbed magnetic and kinetic energies with toroidal harmonics $n=1$ to
$n=6$ for DH-FP-dt0 and DH-QP-dt0 are shown in Fig. 8, where the perturbed
magnetic energies are shown in solid lines and the kinetic ones are shown in
dashed lines. We identified two MHD peaks in both cases, which are marked by
the vertical chained lines in the figure. For both cases, the first peak
corresponds to the fragment arrival on the plasma which results in edge
stochastization, while the second peak corresponds to the core temperature
relaxation. These strong MHD responses correspond to strongest outgoing heat
flux and plasmoid drift motion, thus we would investigate the boundary heat
flux at these times.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 9: The normal component of the combined conductive and convective heat
flux of DH-FP-dt0 at times (a) $t=0.74ms$ and (b) $t=1.17ms$ across a surface
$10cm$ within the simulation boundary, as well as that of DH-QP-dt0 at time
(c) $t=1.09ms$ and (d) $t=1.90ms$. All four time slices coincide with peaks of
MHD amplitude in their respective cases. The EQ-08 port locates approximately
at $\phi=-2.57$ and the EQ-17 port locates at $\phi=0.57$.

The $\log_{10}$ of the normal component of the boundary heat flux at the two
times when the MHD amplitude reaches its peak values for both cases are shown
in Fig. 9. Fig. 9(a) and (b) correspond to two MHD peaks of DH-FP-dt0 and Fig.
9(c) and (d) correspond to that of DH-QP-dt0. For DH-FP-dt0, a small asymmetry
in the MHD response initially results in a weak drift motion and mild out
going heat flux from the EQ-08 plume position at $\phi=-2.57$ and $t=0.74ms$
as is shown in Fig. 9(a). Later on, at the time of the core relaxation, strong
heat flux coming out from the X point of the $1/1$ magnetic island results in
a local high pressure plasmoid close to the EQ-17 plume location at
$\phi=0.57$ and $t=1.17ms$, which is accompanied by a strong outgoing heat
flux shown in Fig. 9(b). The heat flux associated with the second drift is
comparable in amplitude with the heat fluxes in BH-FP-dt0 as is shown in Fig.
6. With higher neon mixture ratio, DH-QP-dt0 found negligible plasmoid drift,
and almost no accompanying outgoing heat flux even during the MHD peak, as is
shown in Fig. 9(c) & (d). This is responsible for the difference in the
radiation fraction between BH-FP-dt0 and DH-QP-dt0 shown in Fig. 7.

To summarize, the arrival of the vanguard fragments on the baseline H-mode
pedestal would trigger strong outward plasmoid drift which is accompanied by
edge stochastization, resulting in significant density expulsion and a large
heat flux onto the wall. This is detrimental to the TQ mitigation. Considering
the confinement degradation in the precursor phase as well as using higher
neon concentration pellets can effectively mitigate such drift motion and the
associated heat flux, thus increasing the radiation fraction.

##  4 Radiation heat flux and its energy impact after neon dual-SPIs

Another important issue that has to be addressed is the asymmetry of the
radiation heat deposition onto the ITER first wall. Concerns have been raised
that if there exists strong asymmetry in the volumetric radiation power
density, local parts of the first wall could still receive significant amount
of heat flux leading to surface melting of PFCs. For the design of the ITER
in-vessel components such as diagnostics or heating systems, the radiation
peaking and the resulting heat fluxes to these components have been specified.
However, validation of these specifications is still outstanding and the
simulations presented here contribute to this validation process. To
investigate such radiation heat deposition, we utilize the IMAS integrated
RaySect/CHERAB code suite to obtain the wall heat flux via ray-tracing the
radiation power density obtained by JOREK onto the realistic ITER first wall
tiles and the stainless steel armour plate of the diagnostic port windows.
Once the heat flux distribution is obtained, the accumulated energy impact is
calculated by the method detailed in Section 2.2, which then is compared to
the maximum tolerable energy impact limit to assess the potential damage to
the PFCs. This serves as a first demonstration of the JOREK capability to
provide input for future more detailed PFC heat damage investigations. It
should be noted that directly using the high time resolution JOREK data to
calculate the energy impact is impractical due to the enormous data size, so
we used lower time resolution data with linear interpolation between time
slices to carry out the time integration needed in Eq. (4).

Case # | Assim. Ne | Assim. H | $f_{rad}$ |  $t_{TQ}$ ($90\%$-$20\%$) | Max. $\Delta Q$
---|---|---|---|---|---
BH-FP-dt0 | $9.91\times 10^{21}$ | $6.40\times 10^{23}$ | $49.4\%$ | $2.0ms$ | $25.4MWs^{1/2}/m^{2}$
DH-QP-dt0 | $2.53\times 10^{22}$ | $4.40\times 10^{23}$ | $89.8\%$ | $4.4ms$ | $5.7MWs^{1/2}/m^{2}$
DH-QP-dt1 | $2.15\times 10^{22}$ | $3.58\times 10^{23}$ | $86.3\%$ | $4.2ms$ | $14.9MWs^{1/2}/m^{2}$
Table 2: The neon and hydrogen assimilation, radiated fraction, TQ time and
maximum local energy impact for three cases considered in this section.

We mainly focus on three cases in this section: BH-FP-dt0, DH-QP-dt0 & DH-QP-
dt1 as shown in Table 1. The overall assimilated neon, assimilated hydrogen
within the LCFS, radiated energy fraction $f_{rad}$, approximate TQ timescale
defined as the timescale of the thermal energy falling from $90\%$ to $20\%$
of the initial thermal energy, as well as the maximum local energy impact for
those three cases are summarized in Table 2. We are interested in DH-QP-dt0 &
DH-QP-dt1 since the quarter pellet dual-SPIs into the degraded H-mode shows
good radiation fraction, as is demonstrated in Fig. 7(b). We are interested in
the comparison between the perfectly synchronized DH-QP-dt0 and the $1ms$
delayed DH-QP-dt1, since the mismatch between the arrival time of the fragment
plumes means stronger asymmetry in the radiation power density and probably a
larger maximum heat flux, and it is of interest to see if the $1ms$ delayed
case would see local radiation energy impact exceeding the tolerable limit of
the PFCs. The baseline H-mode BH-FP-dt0 is included to demonstrate the
expected strongest radiation heat deposition onto the first wall due to the
large radiation power, despite its smaller radiated fraction.

---
(a)
(b)
(c)
Figure 10: The correlation between the radiation power, the neon ablation rate
of both plumes, as well as the $n=1$ and $n=2$ magnetic energy perturbation
for (a) BH-FP-dt0, (b)DH-QP-dt0 and (c)DH-QP-dt1. The approximate arrival time
of the vanguard fragments on the equilibrium LCFS is represented by the green
vertical dashed line for the EQ-08 plume and the magenta vertical dashed line
for the EQ-17 plume. For the DH-QP-dt0 case, the two arrival time overlap with
each other.
---
Figure 11: The ITER first wall (white) and the diagnostic port armour (red). The EQ-08 port locates at $X=-7.05m$ and $Y=4.57m$, it is the third port from right visible in the figure, marked by the green rectangle. The EQ-17 port locates at $X=7.45m$ and $Y=-3.84m$, on the other side of the torus. |
---|---
(a) | (b)
|
(c) | (d)
Figure 12: The first wall heat flux of DH-QP-dt1 at (a) $t=0.52ms$ looking at
port EQ-08, (b) $t=0.90ms$ looking at port EQ-08, (c) $t=5.74ms$ looking at
port EQ-08 and (d) $t=5.74ms$ looking at port EQ-17. By the time of
$t=5.74ms$, the radiation heat flux already becomes relaxed, and despite the
high total radiation power, the radiation heat flux is low.

The correlation between the radiation power, the neon ablation rate of both
plumes, as well as the $n=1$ and $n=2$ magnetic energy perturbation for all
three cases are shown in Fig. 10, with red lines representing the radiation
power, the blue solid lines representing the neon ablation rate of the EQ-08
plume, the blue chained lines representing that of the EQ-17 plume, the black
solid lines representing the $n=1$ perturbed magnetic energy and the black
chained lines representing the $n=2$ perturbed magnetic energy respectively.
The approximate arrival time of the vanguard fragments on the equilibrium LCFS
is represented by the green vertical dashed line for the EQ-08 plume and the
magenta vertical dashed line for the EQ-17 plume. For BH-FP-dt0 shown in Fig.
10(a), the very sharp peak in ablation and radiation at approximately
$t=0.42ms$ coincide with the strong plasmoid drift discussed in Section 3, the
following peak in the MHD amplitude corresponds to the edge stochasticity
shown in Fig. 5. The later radiation peaks coincide with the core temperature
relaxation and cooling, which corresponds a weaker MHD peak just before
$t=1.5ms$. For DH-QP-dt0 shown in Fig. 10(b), there are two radiation peaks in
the early phase of injection at approximately $t=1.1ms$ and $t=2ms$, the later
peak corresponds to the core temperature relaxation. These peaks are closely
correlated with peaks in ablation and MHD amplitude. This correlation is the
result of very nonlinear interaction between the fragments and the MHD mode.
The ablation creates locally peaked lowly charged impurity density profile,
thus result in strong local radiation cooling. This local cooling and local
pressure peak incurs helical current perturbation which destabilize resonant
MHD modes. The growing MHD modes result in enhanced outgoing heat flux which
further strengthen the ablation in the plume. Thus the strong correlation
between the radiation, ablation and MHD amplitude peaks especially in the
early injection phase. DH-QP-dt1 shows similar correlation as is shown in Fig.
10(c).

|
---|---
(a) | (b)
|
(c) | (d)
Figure 13: The first wall energy impact of DH-QP-dt1 calculated using Eq. 4 at
(a) $t=0.56ms$ looking at port EQ-08, (b) $t=0.94ms$ looking at port EQ-08,
(c) $t=5.74ms$ looking at port EQ-08 and (d) $t=5.74ms$ looking at port EQ-17.

For each case, a few time slices are chosen to show both the instantaneous
heat flux and the accumulated energy impact onto the first wall. The time
around the radiation peaks shown in Fig. 10 are especially of interest, so we
choose the time slices to be at the radiation peak for the heat flux analysis
and right after the peak for the energy impact analysis for the early
injection phase. For completion, a few time slices are also chosen for the
late injection phase when the radiation power density distribution is already
relaxed within the plasma. The ITER first wall tiles as well as the diagnostic
window armour plates are shown in Fig. 11, where the tiles are shown in white
and the armour plates are shown in red. The EQ-08 port locates at $X=-7.05m$
and $Y=4.57m$, it is the third port from right visible in the figure. The
EQ-17 port locates at $X=7.45m$ and $Y=-3.84m$, on the other side of the
torus.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 14: The first wall energy impact of DH-QP-dt0 calculated using Eq. 4 at
(a) $t=2.06ms$ looking at port EQ-08, (b) $t=2.06ms$ looking at port EQ-17,
(c) $t=5.64ms$ looking at port EQ-08 and (d) $t=5.64ms$ looking at port EQ-17.

We first take a look at DH-QP-dt1. The instantaneous heat flux onto the first
wall is shown in Fig. 12. Since the EQ-08 plume arrives earlier in this case,
the radiation heat flux concentrates around the EQ-08 port in the early
injection phase as can be seen from Fig. 12(a) and (b), at $t=0.52ms$ and
$t=0.90ms$ respectively. These two times correspond to the first two peaks in
the radiation power shown in Fig. 10(c). The toroidal elongation of the heat
spots in Fig. 12(a) and (b) is due to the toroidally elongated density source
shape as we have discussed in Section 2.1. It can be seen that initially the
radiation heat flux is almost symmetric around the EQ-08 port in Fig. 12(a),
but shifts to the right in Fig. 12(b). This is a result of the interplay
between the impurity density distribution and the outgoing heat flux which
moves the radiation power density peak away from the fragment location that
has been reported previously [8]. Ultimately, the radiation power density
relaxes over the plasma volume both due to the impurity density relaxation and
the arrival of the second plume. As a result, the first wall radiation heat
flux becomes more or less uniform in the late phase of the SPI, as can be seen
in Fig. 12(c) and (d) which show the heat flux around both port EQ-08 and
EQ-17 at $t=5.74ms$. At this time, despite the high total radiation power as
is shown in Fig. 10(c), the maximum radiation heat flux onto the first wall is
one order of magnitude lower compared with the early injection phase shown in
Fig. 12(a) and (b). Another notable feature is that the heat flux on the
armour plate is smaller compared with that of the first wall tiles, which is
due to a geometric effect.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 15: The first wall energy impact of BH-FP-dt0 calculated using Eq. 4 at
(a) $t=0.43ms$ looking at port EQ-08, (b) $t=0.43ms$ looking at port EQ-17,
(c) $t=2.00ms$ looking at port EQ-08 and (d) $t=2.00ms$ looking at port EQ-17.

The energy impact factor at several chosen times for the same case is shown in
Fig. 13. Fig. 13(a) and (b) show the energy impact just after the radiation
spike at $t=0.56ms$ and $t=0.94ms$, corresponding to Fig. 12(a) and (b), and
indeed one can see the energy impact distribution correlates well with that of
the transient heat flux. At the end of the second radiation spike at
$t=0.94ms$, part of the stainless steel armour plate shows the energy impact
approaching the maximum tolerable limit of $13MJs^{-1/2}/m^{2}$ [14]. Another
feature is that the tiles beside the EQ-08 plate show higher energy impact,
consistent with the heat flux behaviour shown in 12(b). Nevertheless, the
energy impact on those tiles is well below the maximum tolerable limit for the
tungsten at $38MJs^{-1/2}/m^{2}$ [11]. In the late injection phase at
$t=5.74ms$, despite the high total radiation power, the energy impact onto the
first wall is much more uniform and stays well below the stainless steel
limit, as is shown in 13(c) and (d) which look at the EQ-08 and EQ-17 port
respectively.

As a comparison, the first wall energy impact in the early and late injection
phase for DH-QP-dt0 is shown in Fig. 14. The energy impact at $t=2.06ms$,
after the second radiation peak shown in Fig. 10(b), around the EQ-08 and
EQ-17 ports is shown in Fig. 14(a) and (b) respectively, while the late time
energy impact at $t=5.64ms$ is shown in Fig. 14(c) and (d) for both locations.
Due to the good synchronization between the injectors, both ports see similar
energy impact accumulation in the early injection phase, and the maximum
energy impact stays well below that of the stainless steel limit, let alone
the tungsten one. In the late injection phase, due to the relaxation of the
radiation power density, the energy impact becomes rather uniform across the
torus apart from some poloidal asymmetry due to the shape of the first wall
cross-section, as can be seen in Fig. 14(c) and (d). Again, despite the higher
radiation power, the maximum energy impact stays well below the tolerable
limit of the stainless steel.

Finally, we take a look at BH-FP-dt0 to demonstrate the highest radiation
energy impact in our simulations. This high energy impact is because of both
the high total radiation power and the asymmetry in the fragment plume we
explained previously. Due to this asymmetry, we expect the radiation to
concentrate around the EQ-08 port since the vanguard of its plume would first
reach the plasma and trigger the plasmoid drift motion discussed in Section 3.
Indeed, as is shown in Fig. 15(a) and (b) at $t=0.43ms$, after the sharp
radiation peak caused by the vanguard fragments arriving on the baseline
H-mode pedestal, significant energy impact is accumulated around the EQ-08
port while the energy impact around the EQ-17 port is negligible. The energy
impact of this radiation pulse already exceeds the tolerable limit of the
stainless steel, thus we expect to see melting of the EQ-08 armour plate even
without considering the conductive/convective heat flux accompanying the
plasmoid drift. The maximum radiation energy impact is still well below the
tungsten limit, however. As time goes on and the fragments fly further inward,
we see the radiation asymmetry turns to the other side of the torus around the
EQ-17 port, as is shown in Fig. 15(c) and (d) at $t=2.00ms$. The radiation
energy impact also concentrates around the injection port at this time, but
shows stronger helical feature, representing the parallel relaxation of the
impurity density. Again, the maximum energy impact exceeds the stainless steel
limit, but stays well below the tungsten one.

##  5 Conclusion and discussion

JOREK simulations with collisional-radiative impurities are carried out for
dual-SPIs into ITER H-mode plasma to investigate the heat deposition during
disruption mitigations. Two sets of equilibria are considered, one baseline
H-mode following the default ITER H-mode scenario, and one degraded H-mode
where an artificial H-L back transition is induced to mimic the confinement
degradation in the disruption precursor phase. Into these two equilibria, two
sets of dual-SPI configuration are simulated, the full pellet case with small
neon fraction and the quarter pellet case where the neon amount is the same
but the hydrogen amount is reduced to one quarter of the full pellet case.

It is found that with baseline H-mode and full pellet, a large local plasmoid
over-pressure occurs near the vanguard fragments, which induces polarization
and consequential drift motion outward along the major radius. This drift
motion causes evident density expulsion and is accompanied by edge
stochasticity and strong outward heat flux, both detrimental to the disruption
mitigation. The transient heat flux is so large that even with expose time on
the order of $1\mu s$, the resulting energy impact could still be on the order
of $10^{8}MJs^{-1/2}/m^{2}$, exceeding the tolerable limit of the tungsten. By
considering the pre-disruption confinement degradation and using the higher
neon fraction quarter pellet, such local over-pressure could be reduced, thus
mitigating the drift motion and the accompanying heat flux. In these
scenarios, the transient heat flux is significantly diminished and the
radiation fraction during the TQ is improved.

The radiation heat flux and its accumulated energy impact onto the first wall
is obtained by ray-tracing from the volumetric radiation power density
distribution provided by JOREK. Three cases are investigated: the baseline
H-mode with full pellet case and two degraded H-mode with quarter pellet cases
with perfectly synchronization and $1ms$ delay between the injectors
respectively. It is found that the radiation heat deposition tends to
concentrate around the injecting port in the early injection phase, while it
spreads over larger areas towards the end of the TQ. With the quarter pellet,
degraded H-mode and $1ms$ injection delay, the maximum radiation energy impact
approaches the tolerable limit of the stainless steel plate while staying far
below the tungsten limit. Good synchronization is found to further relax the
energy impact distribution and reduce its maximum value, such that it never
reaches the stainless steel limit for the perfect synchronized degraded H-mode
case. These results are reassuring since they suggest the degraded H-mode with
quarter pellet case could achieve high radiation fraction without the risk of
melting the first wall. In the baseline H-mode case, despite the strong
radiation power, its energy impact still stays well below the tungsten limit,
although it exceeds that of the stainless steel substantially, hence one would
expect melting of the diagnostic window armour plates. However, even if such
melting occurs, it would not impede continued operation of ITER, as it has
been found in previous studies with $22MJs^{-1/2}/m^{2}$ heat pulses that the
transient surface melting would result in a surface roughening increasing at a
rate of 1-2 micrometers per pulse, but no loss of material [13]. Another
interesting observation is that the maximum energy impact is approaching the
melting limit of the beryllium at $28MJs^{-1/2}/m^{2}$ [14], so that the new
ITER tungsten tiles provide more margin against radiation heat deposition
compared with previous beryllium ones.

A few assumptions in our investigations might impact our results presented
above. First of all is our toroidally elongated ablation density source, which
result in enhanced impurity density spreading along the toroidal direction.
One could expect that a more realistic, and expensive, $\Delta\phi_{NG}$ would
induce stronger impurity density peak and radiation power density peak, thus
more localized radiation heat flux. This would especially be important if the
radiating source is very close to the wall as in the early injection phase, or
if the plasmoid drift causes the radiating plasmoid to move outward towards
the wall. Another consequence of this elongated toroidal density source is
that it tends to result in less concentrated pressure peak since the density
is more spread out, thus reducing the drive for the plasmoid drift. It is
however hard to say if a more realistic toroidal elongation would result in
stronger drift, as the accompanying stronger radiation would tend to reduce
the local pressure, thus counter-act the pressure peaking effect discussed
above. This is left for future studies. Furthermore, we also did not consider
the opacity of the plasma. The high density plasmoid around the fragment plume
could absorb significant portion of the radiation and re-emit them at a later
time [37], reducing the radiation peak. Resolving these two competing factors
is part of our future plan. Last, our choice of viscosities might impact the
TQ duration, as is demonstrated in Ref[10]. The current values we chose are
the minimum value we could use without rampant numerical instabilities. Using
smaller viscosity might results in more rapid MHD growth and quicker TQ.
However, as is shown in Ref[10], changing the viscosity from $2000m^{2}/s$ to
$500m^{2}/s$ only reduces the TQ duration from $2.9ms$ to $2.6ms$ in their
viscosity scan. Hence we don’t expect drastic changes to our conclusion should
a more realistic viscosity is used.

Nevertheless, our study reported above serves as a first look into the
radiation heat deposition on the ITER first wall during its H-mode TQ
mitigation process using non-linear 3D simulation, which is part of ongoing
validation work for the ITER PFC heat load specification. This study also
provides input for more detailed heat transport and material damage studies in
the future.

Acknowledgments

The authors thank Y. Yuan, L. Cheng and A. Matsuyama for fruitful discussion.
ITER is the Nuclear Facility INB no. 174. The views and opinions expressed
herein do not necessarily reflect those of the ITER Organization. This
publication is provided for scientific purposes only. Its contents should not
be considered as commitments from the ITER Organization as a nuclear operator
in the frame of the licensing process. Part of this work is supported by the
National MCF Energy R&D Program of China under Grant No. 2019YFE03010001. This
work has been co-funded by the ITER Organization under the implementing
agreement IO/IA/19/4300002053. This work was supported in part by the Swiss
National Science Foundation. Part of this work has been carried out within the
framework of the EUROfusion Consortium, funded by the European Union via the
Euratom Research and Training Programme (Grant Agreement No. 101052200 -
EUROfusion). The Swiss contribution to this work has been funded by the Swiss
State Secretariat for Education, Research and Innovation (SERI). Views and
opinions expressed are however those of the author(s) only and do not
necessarily reflect those of the European Union, the European Commission or
SERI. Neither the European Union nor the European Commission nor SERI can be
held responsible for them. This work is carried out partly on Tianhe-3
operated by NSCC-TJ and partly on the ITER cluster.

References

## References

  * [1] M. Lehnen, K. Aleynikova, P.B. Aleynikov et al., “Disruptions in ITER and strategies for their control and mitigation” J. Nucl. Mater., 463, 39-48 (2015).
  * [2] Y. Li, Z.Y. Chen, W. Yan et al., “Comparison of disruption mitigation from shattered pellet injection with massive gas injection on J-TEXT”, Nucl. Fusion 61 126025 (2020).
  * [3] U.A. Sheikh, D. Shiraki, R. Sweeney et al., “Disruption thermal load mitigation with shattered pellet injection on the Joint European Torus (JET)”, Nucl. Fusion 61 126043 (2021).
  * [4] S. Jachmich, U. Kruezi, M. Lehnen et al., “Shattered pellet injection experiments at JET in support of the ITER disruption mitigation system design”, Nucl. Fusion 62 026012 (2022).
  * [5] Juhyeok Jang, Jayhyun Kim, Jaewook Kim et al., “Radiation distribution for shattered pellet injection experiment with AXUV array diagnostics in KSTAR”, Fusion Eng. Des. 180 113172 (2022).
  * [6] B. Stein-Lubrano, R. Sweeney, D. Bonfiglio et al., “3D Radiated Power Analysis of JET SPI Discharges Using the Emis3D Forward Modeling Tool”, Nucl. Fusion accepted (2024).
  * [7] C. C. Kim, Y. Liu, P. B. Parks et al., “Shattered pellet injection simulations with NIMROD” Phys. Plasmas, 26, 042510 (2019).
  * [8] D. Hu, E. Nardon, M. Hoelzl et al., “Radiation asymmetry and MHD destabilization during the thermal quench after impurity Shattered Pellet Injection” Nucl. Fusion, 61, 026015 (2021).
  * [9] D. Hu, E. Nardon, F.J. Artola et al., “Collisional-radiative simulation of impurity assimilation, radiative collapse and MHD dynamics after ITER Shattered Pellet Injection” Nucl. Fusion, 63, 066008 (2023).
  * [10] J. McClenaghan, B.C Lyons, C.C. Kim et al., “MHD modeling of shattered pellet injection in JET” Nucl. Fusion, 63, 066029 (2023).
  * [11] A. Herrmann, “Overview on stationary and transient divertor heat loads” Plasma Phys. Control. Fusion 44 (2002) 883–903.
  * [12] Pietro Barabaschi, “Progress on manufacturing, construction, commissioning and an updated baseline”, 29th Fusion Energy Conference, 16–21 Oct 2023, 2023, London, UK.
  * [13] N.S. Klimov, J. Linke, R.A. Pitts et al., “Stainless steel performance under ITER-relevant mitigated disruption photonic heat loads”, J. Nucl. Mater., 438, S241-S245 (2013);
  * [14] R.A. Pitts, B. Bazylev, J. Linke et al., “Final case for a stainless steel diagnostic first wall on ITER” J. Nucl. Mater. 463 (2015) 748–752.
  * [15] J W Coenen, V Philipps, S Brezinsek et al., “Analysis of structural changes and high-heat-flux tests on pre-damaged tungsten from tokamak melt experiments” Phys. Scr. T145 (2011) 014066.
  * [16] G.G. van Eden, T.W. Morgan, H.J. van der Meiden et al., “The effect of high-flux H plasma exposure with simultaneous transient heat loads on tungsten surface damage and power handling” Nucl. Fusion 54 (2014) 123010.
  * [17] Y. Yuan, J. Du, M. Wirtz et al., “Surface damage and structure evolution of recrystallized tungsten exposed to ELM-like transient loads” Nucl. Fusion 56 (2016) 036021.
  * [18] Jun Wang, Chun Li, Yue Yuan et al., “Surface modification of W–V alloy exposed to high heat flux helium neutral beams” Nucl. Fusion 58 (2018) 096001.
  * [19] Yue Yuan, Wangguo Guo, Peng Wang et al., “Influence of surface melting on the deuterium retention in pure and lanthanum oxide doped tungsten” Nucl. Fusion 59 (2019) 016022.
  * [20] M. Hoelzl, G.T.A. Huijsmans, S. Pamela et al., “The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas” Nucl. Fusion 61 065001 (2021).
  * [21] D Hu, G T A Huijsmans, E Nardon et al., “Collisional-radiative non-equilibrium impurity treatment for JOREK simulations” Plasma Phys. Control. Fusion 63 125003 (2021).
  * [22] V. Riccardo, A. Loarte and the JET EFDA Contributors, “Timescale and magnitude of plasma thermal energy loss before and during disruptions in JET”, Nucl. Fusion 45 1427 (2005);
  * [23] V Rozhansky, I Veselova and S Voskoboynikov, “Evolution and stratification of a plasma cloud surrounding a pellet”, Plasma Phys. Control. Fusion 37 399 (1995).
  * [24] A. Matsuyama, R. Sakamoto, R. Yasuhara et al., “Enhanced Material Assimilation in a Toroidal Plasma Using Mixed H2 + Ne Pellet Injection and Implications to ITER” Phys. Rev. Lett. 129, 255001 (2022).
  * [25] M. Carr et al., Raysect Python Raytracing Package (v0.4.0). Zenodo. https://doi.org/10.5281/zenodo.1205064.
  * [26] M. Carr et al., cherab/core : Zenodo. https://doi.org/10.5281/zenodo.3551871
  * [27] S.H. Kim, T.A. Casper and J.A. Snipes, “Investigation of key parameters for the development of reliable ITER baseline operation scenarios using CORSICA”, Nucl. Fusion 58 056013 (2018);
  * [28] F. Villone, Y.Q. Liu, G. Rubinacci et al., “Effects of thick blanket modules on the resistive wall modes stability in ITER”, Nucl. Fusion 50 125011 (2010);
  * [29] J. Zhang and P. B. Parks, “Analytical formula for pellet fuel source density in toroidal plasma configurations based on an areal deposition model”, Nucl. Fusion 60 066027 (2020);
  * [30] S. I. Braginskii, “Transport Processes in a Plasma”, Rev. Plasma Phys. 1, 205, (1965);
  * [31] E.J. Doyle, W.A. Houlberg, Y. Kamada et al., “ITER Progress Chapter 2: Plasma confinement and transport” Nucl. Fusion 47 S18–S127 (2007).
  * [32] L. Spitzer, Physics of Fully Ionized Gases (Interscience, New York, 1956).
  * [33] S.P. Hirshman, “Neoclassical current in a toroidally‐confined multispecies plasma”, Phys. Fluids 21, 1295 (1978).
  * [34] P. B. Parks, “Modeling dynamics Fracture of Cryogenic Pellets”, GA report GA-A28352, (2016).
  * [35] T.E. Gebhart, L.R. Baylor & S.J. Meitner, “Shatter Thresholds and Fragment Size Distributions of Deuterium-Neon Mixture Cryogenic Pellets for Tokamak Thermal Mitigation”, Fusion Sci. Technol., 76:7 831-835 (2020).
  * [36] T. E. Gebhart, L. R. Baylor, and S. J. Meitner, “Analysis of the Shattered Pellet Injection Fragment Plumes Generated by Machine Specific Shatter Tube Designs”, Fusion Sci. Technol., 77:1, 33-41 (2021);
  * [37] Pavel Aleynikov, Alistair M. Arnold, Boris N. Breizman et al., “Thermal quench induced by a composite pellet-produced plasmoid”, Nucl. Fusion 64 016009 (2024);


# A fractal geometry immersed in a hierarchical magnetic flux distribution

Biplab Pal<EMAIL_ADDRESS>Department of Physics, School of
Sciences, Nagaland University, Lumami-798627, Nagaland, India

###### Abstract

Fractal geometry presents us with a self-similarity in their pattern at
various length scales that is prevalent in our natural world. We present
theoretical model of a Sierpinski gasket (SPG) fractal geometry with a
deterministic perturbation in the form of a hierarchical distribution of
magnetic flux. Such flux configuration induces a deterministic disorder in the
Aharonov-Bohm (AB) phases picked up by the electron wavefunction. Using the
tight-binding formalism, we show that by tunning the strength of the hierarchy
parameter of those AB phases, one can systematically engineer quantum states
in a SPG fractal lattice. In addition to this, we have also observed that by
controlling the strength of this hierarchy parameter in the magnetic flux, one
can effectively regulate the persistent current in the SPG fractal structure.
This characteristic is found to be true for various filling factors. Our
results could be useful for designing nanoelectronic devices using molecular
fractal structures fabricated by chemical synthesis technique.

##  I Introduction

The terminology of fractal geometries which was initially emerged purely as a
mathematical concept mandelbrot-book , has generated intense research interest
in various domains of condensed matter physics and material science in recent
times. Although deterministic fractal geometries belongs to the class of
aperiodic systems, they are not completely random; they possess certain self-
similar beautiful pattern and symmetrical structure in their construction.
This often helps theoreticians to analytically investigate various interesting
properties appearing in such fractal structures. Also, deterministic growth
rule allows one to grow such fractal structures in the laboratory using the
state-of-the-art self-assembly or lithography techniques self-assembly-1 ;
self-assembly-2 ; hexagonal-SPG , and perform fascinating experiments with
them.

One of the key reasons why fractal lattice structures have drawn a lot of
research attention of late is due to the fact that they possess non-integer
Hausdorff dimension $D_{f}$. This has ignited a lot of curiosity in the mind
of the researchers and led to a number of intriguing investigations in this
direction, viz., topological character of electronic states on fractals titus
; shreya , flat bands in fractal-like lattices biplab ; atanu ; hahafi-
fractal-FB ; chen-fractal-FB-1 ; chen-fractal-FB-2 ; sougata-fractal-FB , Hall
conductivity askar-Hall ; cm-smith-Hall  and possible existence of anyons in
fractal structures biplab-PRR-2020 , are to name a few of them. The rapid
advancement in experimental tools and fabrication methods has also played a
major role in developing renewed interest in the physics of fractal
geometries. These technological progress has enabled us to synthesize tailor-
made fractal structures in the laboratory using various techniques, such as,
self-assembled molecular growth of fractal materials self-assembly-1 ; self-
assembly-2 , designing quantum fractal using atomic manipulation in a scanning
tunnelling microscope cm-smith-nat-phys , fractal photonic lattice consisting
of waveguides realized using femtosecond-laser-writing technology m-segev
etc.

Inspired by these experimental advancements in creating the fractal structures
in reality, various interesting theoretical studies concerning the
deterministic fractal structures have taken place in recent times. This
includes works on energy level statistics in fractal clusters yorikawa-JPC2019
; askar-PRB2019 , confinement of electrons in fractal plane using electric
field yuan-PRB2020 , electronic fractal pattern in building molecular fractal
systems latge-pccp2022 , role of anisotropy on magneto-transport in fractal
network maiti-PLA2020 , and Josephson effect in a fractal geometry vladimir-
arXiv2024  among others.

Figure 1: Schematic diagram of a Sierpinski gasket fractal lattice with
hierarchical distribution of a magnetic flux following the pattern
$\Phi_{\ell}=\lambda^{\ell}\Phi$, where $\ell$ represents the generation
index.

It is expected that, in future, further advancement in nanotechnology and
fabrication techniques will enable us to do more delicate experiments with
these fractal lattice models. Such possibility inspires us to undertake this
theoretical investigation presented in this article. We consider a Sierpinski
gasket (SPG) fractal lattice structure immersed in a deterministic
hierarchical distribution of magnetic fluxes. The ramification of a SPG
fractal structure is such that, it offers us with down-headed triangular loops
of gradually increasing size (See Fig. 1). Following this ramification
pattern, we introduce a hierarchical distribution of magnetic fluxes
$\Phi_{\ell}=\lambda^{\ell}\Phi$ ($\lambda$ being the hierarchy parameter)
through these down-headed triangular closed loops. Such distribution of the
magnetic fluxes will be manifested as a distribution of Aharonov-Bohm (AB)
phases, whenever an electron will hop from one atomic site to another atomic
site along these down-headed triangular closed loops. This provides us with an
interesting situation, where we have a deterministic disorder in the
distribution of the AB phases in this closed-loop fractal geometry controlled
by a hierarchy parameter $\lambda$.

In this work, we investigate the effect of this usual distribution of magnetic
fluxes on the energy eigenvalue spectrum of the SPG fractal network as a
function of the hierarchy parameter $\lambda$. The presence of the magnetic
flux through such closed-loop fractal lattice geometry induces a circulating
current in these systems known as persistent current maiti-PRB2010 . It is a
pertinent question to ask, may we have some mechanism by which we can control
the persistent current in such closed-loop fractal lattice geometries? In this
study, we address this question and show that, it is indeed possible to
effectively control the magnitude of the persistent current in this SPG
fractal lattice by tunning the hierarchy parameter $\lambda$ in the
distribution of the magnetic flux. Our model may be suitable for mesoscopic or
nanoscopic experiments using fractal structures. In what follows, we present
the schematic and mathematical framework of our model in Sec. II. We discuss
the effect of the hierarchy parameter on the energy eigenvalue spectrum of the
SPG fractal lattice model in Sec. III. The role of the hierarchy parameter in
controlling the persistent current in such complex systems is elaborated in
Sec. IV. Finally, in Sec. V, we draw our conclusion highlighting the summary
of our findings and their utility, and the scope of future study in this
direction.

##  II The model and the mathematical framework

We consider a theoretical model of a SPG fractal lattice which follows certain
deterministic ramification pattern in its growth rule as shown in Fig. 1. The
number of sites in a $\ell$-th generation system is given by the following
formula:

$\mathcal{N}(\ell)=\dfrac{3}{2}\big{(}1+R^{\ell}\big{)},$ (1)

where $R$ is the repetition factor (which is equal to 3 for this particular
fractal model) and $\ell$ denotes the generation index. The single-electron
states in such a lattice structure can be described by the following
Hamiltonian within the tight-binding framework:

$\bm{\mathcal{H}}=\sum_{n}\epsilon_{n}c_{n}^{\dagger}c_{n}+\sum_{\langle
m,n\rangle}\big{(}t_{mn}e^{i\theta_{\ell}}c_{m}^{\dagger}c_{n}+\textrm{h.c.}\big{)},$
(2)

where first term represents the on-site energy and the second term represents
the nearest-neighbor hopping integral.
$\theta_{\ell}=\dfrac{2\pi\Phi_{\ell}}{3(2^{\ell-1})\Phi_{0}}$ is the
Aharonov-Bohm phase A-B-Phase  associated with the nearest-neighbor hopping
parameter $t_{mn}$, whenever an electron hops along a bond in the down-headed
triangular loops pierced by an external magnetic flux
$\Phi_{\ell}=\lambda^{\ell}\Phi$, where $\lambda$ is the hierarchy parameter.
$\Phi$ is measured in units of fundamental flux quantum $\Phi_{0}=hc/e$.
$c_{n}^{\dagger}\,(c_{n})$ denotes the creation (annihilation) operator for an
electron at a site $n$, and $\epsilon_{n}$ is the on-site potential for the
$n$-th atomic site.

Following the Eq. (2), we can construct the Hamiltonian matrix for this
fractal lattice model of different system sizes, which will allow us to
extract the information about the energy spectrum as well as for the
persistent current of this system as a function of the deterministic disorder
distribution in the AB phases. It is to be noted that, because of the self-
similar characteristics of a fractal system, once we get some interesting
properties in such systems for a small system size, it generally retains those
properties even when we go to a large system size. Thus, considering the
possibility of experimental realization of this model with such an unusual
hierarchical distribution of magnetic fluxes in the real-life laboratory, we
consider moderately large system for our theoretical calculations. At first,
we study the role of the hierarchy parameter $\lambda$ which controls the
hierarchical distribution of the magnetic flux $\Phi_{\ell}$ on the energy
eigenvalue spectrum for this model. This is elaborated in the next section.

##  III Effect of the hierarchy parameter on the energy eigenvalue spectrum

In this section, we discuss the role of the hierarchy parameter $\lambda$ on
the energy eigenvalue spectrum of the SPG lattice. We have first explicitly
constructed the Hamiltonian matrix for the SPG lattice using Eq. (2). We take
the on-site potential $\epsilon_{n}=0$ for all the sites, and consider only
nearest-neighbor constant coupling between the sites $t_{mn}=t=1$. The
direction of the hierarchical fluxes $\Phi_{\ell}$ (where $\ell=1,2,3,\ldots$)
in different down-headed triangular closed plaquettes is taken to be in the
out-of-plane direction. Accordingly, the electrons will pick up the AB phases
$\pm\theta_{\ell}$ (‘$+$’ is for forward hopping and ‘$-$’ is for backward
hopping). Once the Hamiltonian matrices are constructed for different
generations of the SPG lattice, one can easily diagonalize them to evaluate
the energy eigenvalue spectrum for the lattice structure. Here, we present the
results for the 3rd and the 4th generation SPG lattice structure in Fig. 2 and
Fig. 3, respectively.

Figure 2: Variation of the energy eigenvalue spectrum for a 3rd generation SPG
lattice structure [with $\mathcal{N}(3)=42$] as a function of the hierarchical
magnetic flux for different values of the hierarchy parameter, viz., (a)
$\lambda=0$, (b) $\lambda=0.1$, (c) $\lambda=0.5$, and (d) $\lambda=0.9$.

We start our discussion with $\lambda=0$, i. e., in absence of the hierarchy
parameter. It is clear that for $\lambda=0$, the hierarchical magnetic flux
distribution in all the plaquettes vanishes. Hence, under this condition, we
get the energy eigenvalue distribution for the SPG lattice geometry in absence
of the magnetic flux. The energy eigenvalue spectrum exhibit highly
degenerate, discrete, fractal characteristics [see Fig. 2(a)]. Then, we
incorporate a small non-zero value of $\lambda=0.1$. With this small non-zero
value of $\lambda$, we instantly observe that, the degeneracy in the energy
spectrum is slightly lifted [Fig. 2(b)]. As we gradually tune the strength of
the hierarchy parameter $\lambda$ more, the energy eigenvalue spectrum become
increasingly non-degenerate [Fig. 2(c) and 2(d)]. Thus, by controlling the
strength of the hierarchy parameter $\lambda$, one can systematically control
the degeneracy of the energy spectrum in a SPG fractal structure.

Figure 3: Energy eigenvalue distribution for a 4th generation SPG lattice
structure [with $\mathcal{N}(4)=123$] as a function of the hierarchical
magnetic flux with values of the hierarchy parameter as (a) $\lambda=0$, (b)
$\lambda=0.1$, (c) $\lambda=0.5$, and (d) $\lambda=0.9$.

It is to be noted that, the self-similar aspect of a fractal structure is, in
general, reflected in all its properties. Therefore, one can expect that, our
observations appearing in Fig. 2 for a 3rd generation SPG fractal structure
will qualitatively persist even when we go to a higher generation fractal
structure. To confirm this fact, we have also evaluated the energy spectrum
for a 4th generation SPG fractal with different values of the hierarchy
parameter $\lambda$. This is depicted in Fig. 3(a)-(d). It is clearly evident
that, the qualitative structure of the energy spectrum remains same as in Fig.
2(a)-(d), apart from appearance of some additional energy eigenvalues which
appear in the spectrum due to increase in the system size. We anticipate the
similar feature to be sustained even for larger SPG fractal structures. It is
worth to point out that, the ramification of the SPG fractal structure is such
that one cannot have such hierarchical distribution of magnetic fluxes in the
upward triangular plaquettes (see Fig. 1). In upward triangular plaquettes,
one can only have a uniform distribution of the magnetic flux as discussed in
Ref. maiti-PRB2010 ; chakrabarti-PRB2005 .

Figure 4: Average density of states (ADOS) for a 3rd generation SPG fractal
structure as a function of the energy and hierarchical magnetic flux. (a) is
for $\lambda=0$ and (b) is for a non-zero value of $\lambda=0.9$.

Figure 5: Variation of the persistent current ($I_{\textrm{max}}$) as a
function of the hierarchy parameter $\lambda$ for a 3rd generation SPG fractal
lattice. (a) is for a filling factor $\nu=1/3$, (b) is for $\nu=1/2$, and (c)
represents $\nu=2/3$.

To corroborate our results on the energy spectrum of the SPG fractal lattice
as a function of the hierarchy parameter, we have calculated the average
density of states (ADOS) for a finite generation SPG fractal structure. To
compute the ADOS, we have used the Green’s function technique which people
uses in general for such fractal lattice models nandy-JPCM2015 ; pal-PRB2012 ;
pal-PRA2013 ; pal-EPJB2012 . Using the Green’s function technique, ADOS can be
defined as,

$\rho(E,\Phi)=-\dfrac{1}{\mathcal{N}\pi}\textrm{Im}\Big{[}\textrm{Tr}\left[\bm{\mathcal{G}}(E,\Phi)\right]\Big{]},$
(3)

where
$\bm{\mathcal{G}}(E,\Phi)=\left[z^{+}{\bm{I}}-\bm{\mathcal{H}}\right]^{-1}$ is
the Green’s function of the system with $z^{+}=E+i\eta\ (\eta\rightarrow
0^{+})$, $\mathcal{N}$ is the number of sites in the lattice, and ‘Tr’
represents the trace of the Green’s function $\bm{\mathcal{G}}$. The results
are displayed in Fig. 4. For $\lambda=0$, we find that, the ADOS spectrum is
not very densely populated [Fig. 4(a)]; whereas for a non-zero value of
$\lambda=0.9$, the ADOS spectrum becomes very densely packed [Fig. 4(b)]. This
indicates that, with a non-zero value of the hierarchy parameter $\lambda$,
one can create more number of allowed states for the electrons to live in such
a deterministic fractal geometry. This suggests towards a possible mechanism
to engineer quantum states in a SPG fractal geometry.

It is important to see the manifestation of these properties in terms of some
experimentally measurable quantities. This brings us to the next section,
where we demonstrate how one can control the persistent current as a function
of the hierarchy parameter in this closed-loop fractal geometry.

##  IV Controlling the persistent current

In isolated mesoscopic metallic rings threaded by a magnetic flux, one can
have a dissipationless circulating current termed as persistent current
buttiker-PLA1983 ; gefen-PRL1989 . It is important to emphasize the fact that,
this current has nothing to do with the persistent current found in a
superconducting material. This idea was initially proposed for a single
mesoscopic conducting loop structure gefen-PRL1989 ; georges-bookchapter2007 ,
but later on people have also investigated the behavior of such persistent
current in multi-loop fractal geometries maiti-PRB2010 . The persistent
current as a function of the magnetic flux $\Phi$ is related to the energy
eigenvalues (energy levels) $E_{\alpha}$ at absolute zero temperature ($T=0$
K) as follows:

$I_{\alpha}(\Phi)=-c\;\dfrac{\partial E_{\alpha}(\Phi)}{\partial\Phi},$ (4)

where $c$ is a constant which is set equal to 1 for our calculation. The total
persistent current in the system is given by,

$I(\Phi)=\sum_{\alpha=1}^{\mathcal{N}_{e}}I_{\alpha}=-\sum_{\alpha=1}^{\mathcal{N}_{e}}\dfrac{\partial
E_{\alpha}}{\partial\Phi},$ (5)

where $\mathcal{N}_{e}$ is the total number of electrons in the system. We
define $\nu=\mathcal{N}_{e}/\mathcal{N}(\ell)$ as the “filling factor” for a
$\ell$-th generation SPG fractal system. Note that, we consider non-
interacting spinless electron model.

For our model, the magnetic flux $\Phi$ depends on the hierarchy parameter
$\lambda$. So, we can study the magnetic response of the SPG fractal geometry
in terms of the persistent current as a function of $\lambda$. We have
systematically investigated the behavior of the persistent current in a 3rd
generation SPG fractal network as a function of the hierarchy parameter
$\lambda$ for different values of $\nu$. This is illustrated in Fig. 5(a)-(c).
It is observed that, for different filling factors $\nu$, the persistent
current in the system gets enhanced as one increases the value of $\lambda$.
This opens up an interesting possibility by which one can flexibly control the
persistent current in such multi-loop fractal geometries. Once again, this
feature is also tested for a 4th generation SPG fractal structure (Fig. 6),
and found to be true as an evidence of the generic self-similar property of a
fractal structure.

Figure 6: Variation of the persistent current ($I_{\textrm{max}}$) as a
function of the hierarchy parameter $\lambda$ for a 4th generation SPG fractal
structure with (a) $\nu=1/3$ and (b) $\nu=2/3$.

It is to be noted that, in principle, one could study these properties for any
other higher generation SPG fractal structures. We predict that, our results
will be qualitatively retained for those cases also. However, from the
experimental point-of-view, one cannot go beyond certain system size in the
real-life laboratory when it comes to fabricating such complex fractal
geometries.

##  V Summary and future outlook

In conclusion, we have studied the role of a hierarchical distribution of
magnetic flux on the energy spectrum and persistent current in a multi-loop
Sierpinski gasket fractal structure. We have shown that, by controlling the
hierarchy parameter of the embedded magnetic flux in this SPG fractal network,
one can effectively engineer quantum states in such a complex system. From the
application point-of-view, we have shown that this directly impacts the
behavior of the persistent current in this SPG fractal structure. We can
suitably control the persistent current in the system by tunning the hierarchy
parameter in a flexible way for different filling factors of the electronic
states. One could in principle, test these results by designing an experiment
using a self-assembled molecular SPG fractal structure. This could be helpful
to devise an unconventional nanoelectronic device using a molecular SPG
fractal structure. In future, we plan to generalize this idea for other
similar kind of closed-loop fractal geometries, such as Vicsek fractal pal-
PRB2012 ; pal-PRA2013 .

###### Acknowledgements.

The author would like to thank Prof. Arunava Chakrabarti from Presidency
University, Kolkata for his useful feedback on the initial stage of this work.
The author also acknowledges helpful discussion with Amrita Mukherjee from
TIFR, Mumbai on the calculation of the persistent current.

## References

  * (1) B. B. Mandelbrot, The Fractal Geometry of Nature (W. H. Freeman, San Francisco, 1982).
  * (2) X. Zhang, N. Li, G.-C. Gu, H Wang, D. Nieckarz, P. Szabelski, Y. He, Y. Wang, C. Xie, Z.-Y. Shen, J.-T. Lü, H. Tang, L.-M. Peng, S.-M. Hou, K. Wu, and Y.-F. Wang, ACS Nano 9, 11909 (2015).
  * (3) J. Shang, Y. Wang, M. Chen, J. Dai, X. Zhou, J. Kuttner, G. Hilt, X. Shao, J. M. Gottfried, and K. Wu, Nat. Chem. 7, 389 (2015).
  * (4) G. R. Newkome, P. Wang, C. N. Moorefield, T. J. Cho, P. P. Mohapatra, S. Li, S.-H. Hwang, O. Lukoyanova, L. Echegoyen, J. A. Palagallo, V. Iancu, and S.-W. Hla, Science 312, 1782 (2006).
  * (5) M. Brzezińska, A. M. Cook, and T. Neupert, Phys. Rev. B 98, 205116 (2018).
  * (6) S. Pai and A. Prem, Phys. Rev. B 100, 155135 (2019).
  * (7) B. Pal and K. Saha, Phys. Rev. B 97, 195101 (2018).
  * (8) A. Nandy, Phys. Scr. 96, 045802 (2021).
  * (9) H. Hanafi, P. Menz, and C. Denz, Adv. Optical Mater. 10, 2102523 (2022).
  * (10) Y. Xie, L. Song, W. Yan, S. Xia, L. Tang, D. Song, J.-W. Rhim, and Z. Chen, APL Photonics 6, 116104 (2021).
  * (11) L. Song, Y. Xie, S. Xia, L. Tang, D. Song, J.-W. Rhim, and Z. Chen, Laser Photonics Rev. 17, 2200315 (2023).
  * (12) S. Biswas and A. Chakrabarti, Physica E 153, 115762 (2023).
  * (13) A. A. Iliasov, M. I. Katsnelson, and S. Yuan, Phys. Rev. B 101, 045413 (2020).
  * (14) M. Fremling, M. van Hooft, C. M. Smith, and L. Fritz, Phys. Rev. Research 2, 013044 (2020).
  * (15) S. Manna, B. Pal, W. Wang, and A. E. B. Nielsen, Phys. Rev. Research 2, 023401 (2020).
  * (16) S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J. M. Zevenhuizen, D. Vanmaekelbergh, I. Swart, and C. M. Smith, Nat. Phys. 15, 127 (2019).
  * (17) Z. Yang, E. Lustig, Y. Lumer, and M. Segev, Light Sci. Appl. 9, 128 (2020).
  * (18) H. Yorikawa, J. Phys. Commun. 3, 085004 (2019).
  * (19) A. A. Iliasov, M. I. Katsnelson, and S. Yuan, Phys. Rev. B 99, 075402 (2019).
  * (20) X. Yang, W. Zhou, P. Zhao, and S. Yuan, Phys. Rev. B 102, 245425 (2020).
  * (21) L. L. Lage and A. Latgé, Phys. Chem. Chem. Phys. 24, 19576 (2022).
  * (22) B. Mal, M. Banerjee, and S. K. Maiti, Phys. Lett. A 384, 126378 (2020).
  * (23) M. Amundsen, V. Juričić, and J. A. Ouassou, arXiv:2404.01373 [cond-mat.supr-con] (2024).
  * (24) S. K. Maiti and A. Chakrabarti, Phys. Rev. B 82, 184201 (2010).
  * (25) Y. Aharonov and D. Bohm, Phys. Rev. 115, 485 (1959).
  * (26) A. Chakrabarti, Phys. Rev. B 72, 134207 (2005).
  * (27) A. Nandy, B. Pal, and A. Chakrabarti, J. Phys.: Condens. Matter 27, 125501 (2015).
  * (28) B. Pal and A. Chakrabarti, Phys. Rev. B 85, 214203 (2012).
  * (29) B. Pal, P. Patra, J. P. Saha, and A. Chakrabarti, Phys. Rev. A 87, 023814 (2013).
  * (30) B. Pal and A. Chakrabarti, Eur. Phys. J. B 85, 307 (2012).
  * (31) M. Büttiker, Y. Imry, and R. Landauer, Phys. Lett. A 96, 365 (1983).
  * (32) H.-F. Cheung, E. K. Riedel, and Y. Gefen, Phys. Rev. Lett. 62, 587 (1989).
  * (33) G. Bouzerar, in Physics of Zero- and One-Dimensional Nanoscopic Systems, Springer Series in Solid-State-Sciences, Vol. 156, edited by S. N. Karmakar, S. K. Maiti, and J. Chowdhury Springer-Verlag, Berlin, 2007, p. 229.