# Millimeter Observations of the Type II SN 2023ixf: Constraints on the
Proximate Circumstellar Medium

Edo Berger Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA Edo Berger<EMAIL_ADDRESS>Garrett K. Keating Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  Raffaella Margutti Department of Astronomy, University of California, Berkeley, CA, 94720  Department of Physics, University of California, Berkeley, CA, 94720  Keiichi Maeda Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan  Kate D. Alexander Department of Astronomy/Steward Observatory, 933 North Cherry Avenue, Rm. N204, Tucson, AZ 85721-0065, USA  Yvette Cendes Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  Tarraneh Eftekhari NASA Einstein Fellow  Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA  Mark Gurwell Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  Daichi Hiramatsu Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  Anna Y. Q. Ho Department of Astronomy, Cornell University, Ithaca, NY 14853, USA  Tanmoy Laskar Department of Physics & Astronomy, University of Utah, Salt Lake City, UT 84112, USA  Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL, Nijmegen, The Netherlands  Ramprasad Rao Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  Peter K. G. Williams Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA

(Received TBD; Revised TBD; Accepted TBD)

###### Abstract

We present 1.3 mm (230 GHz) observations of the recent and nearby Type II
supernova, SN 2023ixf, obtained with the Submillimeter Array (SMA) at
$2.6-18.6$ days after explosion. The observations were obtained as part the
SMA Large Program POETS (Pursuit of Extragalactic Transients with the SMA). We
do not detect any emission at the location of SN 2023ixf, with the deepest
limits of $L_{\nu}(230\,{\rm GHz})\lesssim 8.6\times 10^{25}$ erg s-1 Hz-1 at
2.7 and 7.7 days, and $L_{\nu}(230\,{\rm GHz})\lesssim 3.4\times 10^{25}$ erg
s-1 Hz-1 at 18.6 days. These limits are about a factor of 2 times dimmer than
the mm emission from SN 2011dh (IIb), about an order of magnitude dimmer
compared to SN 1993J (IIb) and SN 2018ivc (IIL), and about 30 times dimmer
than the most luminous non-relativistic SNe in the mm-band (Type IIb/Ib/Ic).
Using these limits in the context of analytical models that include
synchrotron self-absorption and free-free absorption we place constraints on
the proximate circumstellar medium around the progenitor star, to a scale of
$\sim 2\times 10^{15}$ cm, excluding the range $\dot{M}\sim{\rm few}\times
10^{-6}-10^{-2}$ M⊙ yr-1 (for a wind velocity, $v_{w}=115$ km s-1, and ejecta
velocity, $v_{\rm eje}\sim(1-2)\times 10^{4}$ km s-1). These results are
consistent with an inference of the mass loss rate based on optical
spectroscopy ($\sim 2\times 10^{-2}$ M⊙ yr-1 for $v_{w}=115$ km s-1), but are
in tension with the inference from hard X-rays ($\sim 7\times 10^{-4}$ M⊙ yr-1
for $v_{w}=115$ km s-1). This tension may be alleviated by a non-homogeneous
and confined CSM, consistent with results from high-resolution optical
spectroscopy.

Supernovae (1668); Core-collapse supernovae (304); Type II supernovae (1731);
Massive stars (732); Stellar mass loss (1613) Circumstellar matter (241)

††journal: ApJ††facilities: SMA

TBD

##  1 Introduction

Over the past few decades observations of core-collapse supernovae (CCSNe)
ranging from radio to X-rays have revealed that massive stars may undergo
enhanced or eruptive mass loss in the final years and decades prior to core-
collapse (e.g., Smith 2014). The resulting circumstellar medium (CSM) around
the progenitor stars can be probed using optical spectroscopy shortly after
explosion, which can reveal emission from the ionized CSM (so-called flash
spectroscopy; e.g., Gal-Yam et al. 2014; Yaron et al. 2017); optical/UV
observations of the rising light curves, which can reveal energy deposition
from CSM-SN ejecta shock interaction (e.g., Nakar & Piro 2014; Moriya et al.
2018; Hiramatsu et al. 2021); X-ray observations, which track thermal
bremsstrahlung emission from CSM-SN ejecta shock interaction and absorption by
neutral CSM; and radio/millimeter observations which trace synchrotron
emission due to the acceleration of relativistic electrons in CSM-SN shock
interaction (e.g., Chevalier & Fransson 2017).

SN 2023ixf was discovered on 2023 May 19.727 UT shortly after explosion in
M101 (Itagaki, 2023), at $d\approx 6.9$ Mpc (Riess et al., 2022), and was
classified as a Type II SN shortly thereafter (Perley et al., 2023). Thanks to
its early discovery and proximity, extensive data across the electromagnetic
spectrum is being collected, including early optical spectra that reveal flash
ionization features (Jacobson-Galan et al., 2023; Smith et al., 2023) and
narrow absorption due to the pre-shocked CSM (Smith et al., 2023); optical
light curves, which point to early excess emission (Hosseinzadeh et al.,
2023); and hard X-ray observations, indicating a large neutral hydrogen column
density in the vicinity of the SN (Grefenstette et al., 2023). Pre-explosion
Hubble Space Telescope, Spitzer Space Telescope, and ground-based observations
point to a variable and dusty red supergiant progenitor (Kilpatrick et al.,
2023; Neustadt et al., 2023; Pledger & Shara, 2023). Here, we report on early
observations obtained with the SMA spanning about 2.7 to 18.6 days after the
estimated time of first light, and use these data to constrain the CSM density
around the progenitor star. We present the data in §2 and compare to previous
mm-band data for CCSNe in §3 and to numerical and analytical models in §4. We
discuss and summarize the results in §5.

##  2 Observation and Data Analysis

Following the discovery and classification of SN 2023ixf, we used the SMA on
several occasions to observe the SN, starting on 2023 May 21.17, about 2.4
days after the estimated time of first light (2023 May 18.74 UT; e.g.,
Hosseinzadeh et al. 2023; Hiramatsu et al. in prep.) The observations were
obtained as part of the new Large Project POETS (Pursuit of Extragalactic
Transients with the SMA; project 2022B-S046, PI: Berger). The observations are
summarized in Table 1. During these observations, the SMA was tuned to an LO
frequency of 225.5 GHz, providing spectral coverage between 209.5–221.5 and
229.5–241.5 GHz. Across all nights, 3C454.3 was observed as a bandpass
calibrator, Ceres was observed as a flux calibrator, and J1419+543 and
J1506+426 were observed as gain calibrators, with a 15-minute cycle time
cadence.

Figure 1: Map from the aggregate SMA data from all 6 epochs of observations,
of the region around SN 2023ixf (white cross) overlaid on a false-color Hubble
Space Telescope image, with the synthesized beam shown in the lower-left
(white ellipse). We show contours of the 230 GHz (1.3 mm) continuum emission
in red, and the CO(2-1) line emission in cyan. Continuum contours correspond
to 0.6, 0.9, 1.2, and 1.5 mJy; CO(2-1) contours correspond to 2, 4, 8 Jy km
s-1 (uncorrected for primary beam effects). No continuum emission is detected
at the position of SN 2023ixf, which is located near NGC 5461, an H II region
at the southeastern edge of M101 (white box in inset). Table 1: SMA
Observations of SN 2023ixf

UT Date | $\delta t$ | $F_{\nu}\,^{a}$ | $L_{\nu}\,^{a}$
---|---|---|---
| (d) | (mJy) | ($10^{26}$ erg s-1 Hz-1)
May 21.17–21.63 | 2.66 | $<1.5$ | $<0.86$
May 22.17–22.29 | 3.49 | $<3.3$ | $<1.89$
May 24.51–24.62 | 5.83 | $<4.2$ | $<2.40$
May 26.34–26.62 | 7.74 | $<1.5$ | $<0.86$
May 28.44–28.62 | 9.79 | $<2.4$ | $<1.37$
June 6.09–6.60 | 18.60 | $<0.6$ | $<0.34$

Note. — Phases ($\delta t$) are given relative to an estimated time of first
light of 2023 May 18.74 UT and are at the mid-point of each observation.

a Limits are 3 times the image rms.

Analysis of the data was performed using the SMA COMPASS pipeline (Keating et
al., in prep.), which flags spectral data based on outliers in amplitude when
coherently averaging over increasing time intervals for each channel within
each baseline, as well as baselines where little-to-no coherence is seen on
calibrator targets. Flux calibration was performed using the Butler-JPL-
Horizons 2012 (Butler, 2012) model for Ceres. The data were imaged, and
deconvolution was performed via the CLEAN algorithm (Högbom, 1974).

We do not detect emission from the location of SN 2023ixf in any of our
observations, with root-mean-square noise levels spanning about $0.2-1.4$ mJy.
The corresponding luminosity limits are $\lesssim 3.4\times 10^{25}-2.4\times
10^{26}$ erg s-1 Hz-1 (Table 1). An SMA continuum and CO(2-1) line map are
shown in Figure 1 overlaid on a Hubble Space Telescope image of the field.

Figure 2: SMA 230 GHz upper limits for SN 2023ixf (black triangles; $3\sigma$)
compared to the existing sample of CCSNe with early mm-band detections
(squares: Type II; circles: Type Ib/c; Weiler et al. 2007; Soderberg et al.
2008; Gorosabel et al. 2010; Horesh et al. 2013a, b; Chakraborti et al. 2013;
Zauderer et al. 2013, 2014; Maeda et al. 2021).

##  3 Comparison to Other Core-Collapse Supernovae

In Figure 2 we show the SMA spectral luminosity upper limits for SN 2023ixf in
comparison to the mm-band ($100-250$ GHz) light curves of several nearby CCSNe
(Type II and Ib/c) on timescales of $\sim 3-23$ days. Previously detected SNe
span luminosities from only a factor of about 2 times higher than our limits
(e.g., SN 2011dh; Horesh et al. 2013c) and up to $\sim 3\times 10^{27}$ erg
s-1 Hz-1, more than an order of magnitude more luminous than our limits for SN
2023ixf (e.g., SN 2018ivc; Maeda et al. 2023). In recent work Maeda et al.
(2023) inferred mass loss rates (scaled to a wind velocity of 100 km s-1) of
$\sim 10^{-3}$ M⊙ yr-1 for SN 2018ivc, $\sim 10^{-4}$ M⊙ yr-1 for SN 1993J,
and $\sim 10^{-5}$ M⊙ yr-1 for SN 2011dh based on mm-band data. Comparable
mass loss rates, $\sim{\rm few}\times 10^{-6}-10^{-4}$ M⊙ yr-1 have been
inferred for other Type IIb SNe based on cm-band data (e.g., Nayana et al.
2022 and references therein). We note that all these comparison SNe are of
Type IIb or IIL withmore stripped and compact progenitors than that of SN
2023ixf (Kilpatrick et al., 2023; Neustadt et al., 2023) and likely have
higher velocity ejecta that boost their mm-band (and radio) synchrotron
emission.

##  4 CSM Interaction Models

We use the non-detections at 230 GHz to place constraints on the CSM around
the progenitor of SN 2023ixf on a scale of $\sim 2\times 10^{14}-2\times
10^{15}$ cm (for an assumed ejecta velocity of $\sim 10^{4}$ km s-1). In the
mm band the expected emission is due to synchrotron radiation arising from
shock interaction of the SN ejecta with the CSM. We consider the effects of
both synchrotron self-absorption and free-free absorption by an external CSM,
following the well established models of Weiler et al. (1986) and Chevalier
(1998), which have been used to interpret the radio/mm emission from previous
CCSNe.

In the high-density limit, free-free absorption is expected to dominate, with
the optical depth given by Weiler et al. (1986): $\tau_{\rm
ff}=K_{2}(\nu/5\,{\rm GHz})^{-2.1}t_{d}^{-3}$, where $t_{d}$ is time in days
since explosion, and we have assumed that the ejecta are roughly in free
expansion (leading to the temporal power law index of $-3$). Setting
$\tau_{\rm ff}(230\,{\rm GHz})\gtrsim 1$ at 2.7 and 18.6 days we find
$K_{2}\gtrsim 6.1\times 10^{4}$ and $\gtrsim 2\times 10^{7}$, respectively.
The mass loss rate is given by (Weiler et al., 1986):

$\dot{M}\approx 10^{-8}\,M_{\odot}\,{\rm yr}^{-1}\,\tau_{\rm ff,5\,{\rm
GHz}}^{0.5}\,t_{d}^{1.5}\,(v_{\rm ej}/10^{4}\,{\rm km\,s^{-1}})^{1.5}\,\\\
\times(T_{e}/10^{4}\,{\rm K})^{0.68}\,(v_{w}/10\,{\rm km\,s^{-1}}),$ (1)

where $T_{e}$ is the electron temperature, and $v_{w}$ is the CSM wind
velocity. Using our SMA limits at 2.7 and 18.6 days, with $v_{w}=115$ km s-1
(Smith et al., 2023), $v_{\rm ej}=10^{4}$ km s-1, and $T_{e}=10^{4}$ K we find
$\dot{M}\gtrsim 10^{-4}$ and $\gtrsim 4\times 10^{-2}$ M⊙ yr-1, respectively.

In the case of synchrotron self-absorption only, the inferred mass loss rate
is much lower, and is given by (Chevalier, 1998):

$\dot{M}\approx 8.2\times 10^{-7}\,M_{\odot}\,{\rm
yr}^{-1}\,\epsilon_{B,-1}^{-1}\,(\epsilon_{e}/\epsilon_{B})^{-8/19}\,L_{26}^{-4/19}\,\\\
\times\nu_{230\,{\rm GHz}}^{2}\,t_{d}^{2}\,(v_{w}/10\,{\rm km\,s^{-1}}),$ (2)

where $\epsilon_{e}$ and $\epsilon_{B}$ are the post-shock energy fractions in
the relativistic electrons and magnetic fields, respectively. Using our limits
at 2.7 and 18.6 days, and assuming $\epsilon_{e}=\epsilon_{B}=0.1$, we find
$\dot{M}\lesssim 7\times 10^{-5}$ and $\lesssim 4\times 10^{-3}$ M⊙ yr-1,
respectively.

A full calculation of the effects of free-free absorption and synchrotron
self-absorption (following Chevalier 1998) is shown in Figure 3, where we
constrain the joint phase-space of $\dot{M}$ and $v_{\rm ej}$ ruled out by our
SMA data. In this calculation we determine the mm-band luminosity for a range
of ejecta velocities and mass loss rates using an input kinetic energy,
$v_{w}=115$ km s-1, and assuming two sets of equipartition values
($\epsilon_{e}=\epsilon_{B}=0.1$ and $\epsilon_{e}=0.1,\epsilon_{B}=0.01$) and
two values for $T_{e}$ ($10^{4}$ and $10^{5}$ K). We find that for a
reasonable range of $v_{\rm ej}=(1-2)\times 10^{4}$ km s-1 we can exclude a
mass loss rate range of $\dot{M}\sim{\rm few}\times 10^{-6}-10^{-2}$ M⊙ yr-1
(for $\epsilon_{e}=\epsilon_{B}=0.1$ and $T_{e}=10^{5}$ K).

Figure 3: The phase-space of progenitor mass-loss rate versus ejecta velocity
that is excluded by our most constraining SMA limits at 2.7, 7.7, and 18.6
days. The models (Chevalier, 1998) assume $v_{w}=115$ km s-1 (Smith et al.,
2023) and are shown for two sets of equipartition parameters ($\epsilon_{e}$
and $\epsilon_{B}$), an electron power law distribution with
$N_{e}(\gamma_{e})\propto\gamma_{e}^{-p}$ with $p=3$, a volume filling factor
of 0.5, and two values for the electron temperature ($T_{e}$). The vertical
grey shaded region indicates an ejecta velocity range of $v_{\rm
ej}=(1-2)\times 10^{4}$ km s-1. The blue dot marks the value of $\dot{M}$
inferred from X-ray observations acquired at $\lesssim 11$ days (Grefenstette
et al., 2023), while the horizontal blue dot-dashed line marks the value of
$\dot{M}$ inferred from optical spectra at $\lesssim 14$ days (Jacobson-Galan
et al., 2023). Both values have been updated to $v_{w}=115$ km s-1.

Figure 4: Model light curves at 230 GHz (lines) using the formalism of Maeda
et al. (2021) and Maeda et al. (2023) for several choices of the mass loss
rate and for two sets of equipartition parameters:
$\epsilon_{e}=\epsilon_{B}=0.1$ (Left) and $\epsilon_{\rm e}=0.04$ and
$\epsilon_{\rm B}=0.02$ (Right; calibrated with models for SN 2020oi: Maeda et
al. 2023). The triangles are the upper limits for SN 2023ixf.

The constraints shown in Figure 3 are demonstrated by a comparison between
model light curves and the observed upper limits in Figure 4. Here, we use the
model prescriptions described in Maeda et al. (2021) and Maeda et al. (2023).
We adopt a CSM with a wind profile truncated at $2\times 10^{15}$ cm to
represent a “confined” CSM. The ejecta mass and energy are set to be
$11M_{\odot}$ and $10^{51}$ erg, and free-free absorption is included assuming
$T_{\rm e}=10^{5}$ K. Inverse Compton cooling is taken into account, using the
$r/R$-band light curve (Jacobson-Galan et al., 2023) as a proxy for the
bolometric light curve of SN 2023ixf. While the model here adopts the self-
similar solution for the shock-wave dynamics (Chevalier, 1998), it reproduces
the results of the more detailed calculation by Matsuoka et al. (2019).

Adopting our fiducial values of $\epsilon_{\rm e}=\epsilon_{\rm B}=0.1$, we
find that the 230 GHz upper limits for SN 2023ixf rule out $\dot{M}\sim
5\times 10^{-6}-2\times 10^{-3}$ M⊙ yr-1, comparable to the model exclusion
region shown in Figure 3. A key point is that the model accounts for the
deceleration of the outer ejecta, which is important in the case of the high
mass loss rate regime (this means that we need to evaluate the constraints on
mass loss rate at somewhat different velocities for the lower and upper bounds
in Figure 3). Additionally, we stress that the choice of (unconstrained)
microphysical parameters ($\epsilon_{\rm e}$ and $\epsilon_{\rm B}$) are
important, as demonstrated in the right panel of Figure 4 where we adopt the
smaller values inferred for SN 2020oi (Maeda et al., 2021). In this case, the
allowed range of the mass-loss rate is $\dot{M}\lesssim 5\times 10^{-5}$ M⊙
yr-1 or $\gtrsim 2\times 10^{-3}M_{\odot}$ yr-1; interestingly, an
intermediate mass-loss rate of $\sim 10^{-4}$ M⊙ yr-1 is marginally allowed.
In any case, it is important to use additional information such as the inverse
Compton cooling or multi-wavelength information to robustly constrain the
mass-loss rate (e.g., Maeda et al., 2023).

In recent work, Jacobson-Galan et al. (2023) estimated a mass loss rate of
$\dot{M}\sim 10^{-2}$ M⊙ yr-1 to a scale of $\sim 10^{15}$ cm by comparing the
ionization features in early optical spectra (extending to about 14 days) with
spectral models of CSM interaction. This result is marked in Figure 3 and is
in agreement with the results presented here. On the other hand, Grefenstette
et al. (2023) used hard X-ray data at 4 and 11 days to estimate $\dot{M}\sim
3\times 10^{-4}$ M⊙ yr-1 on a similar radial scale (or $\sim 7\times 10^{-4}$
M⊙ yr-1 if using $v_{w}=115$ km s-1; see Figure 3). This value is in tension
with our results for $T_{e}=10^{5}$ K and marginally for $T_{e}=10^{4}$ K.
However, it is essential to note that each of these approaches to determining
the CSM density makes key simplifying assumptions that may lead to disparate
estimates of the CSM properties. This includes numerical factors (e.g.,
$\epsilon_{e}$, $\epsilon_{B}$, and $T_{e}$ in our case) and underlying
geometrical assumptions such as spherical symmetry and a homogeneous CSM. For
example, it is possible that a non-homogeneous and radially-confined CSM,
which may also affect the ejecta velocity angular profile, would lead to X-ray
and mm emission regions of different densities. Indeed, an asymmetric and
confined CSM on a radial scale of $\sim 10^{15}$ cm has been inferred from
high resolution optical spectroscopy (Smith et al., 2023).

##  5 Conclusions

We have presented mm-band observations of the recently-discovered SN 2023ixf,
covering about 2.7 to 18.6 days post-explosion. The non-detections place an
upper bound on the luminosity at 230 GHz of $\lesssim 8.6\times 10^{25}$ erg
s-1 Hz-1 at 2.7 and 7.7 days, and $\lesssim 3.4\times 10^{25}$ erg s-1 Hz-1 at
18.6 days. These limits are about a factor of 2 times lower than the mm-band
emission detected from SN 2011dh (IIb) and about an order of magnitude lower
than the emission in SN 1993J (IIb) and SN 2018ivc (IIL). Using these limits
we place constraints on the proximate CSM of the progenitor of SN 2023ixf, out
to a scale of $\sim{\rm few}\times 10^{15}$ cm of $\dot{M}\gtrsim 10^{-2}$ M⊙
yr-1 (free-free absorption) or $\lesssim 10^{-6}$ M⊙ yr-1 (synchrotron self-
absorption). We note that these limits are in agreement with the inferences
from early optical spectroscopy (Jacobson-Galan et al., 2023; Smith et al.,
2023) but are in tension with the inferred mass loss based on early X-ray
observations (Grefenstette et al., 2023). A likely scenario that may alleviate
this tension is a non-homogeneous and confined dense CSM (Smith et al., 2023),
which can lead to the X-ray emission emerging from lower density regions while
the mm-band emission is absorbed.

We anticipate that continued observations of SN 2023ixf in the radio/mm and
X-ray regimes over the coming weeks, months, and years, as well as joint
modeling with the available optical photometry and spectroscopy, will
eventually better delineate the complex and likely non-spherical and non-
homogeneous CSM environment of the progenitor. More broadly, as demonstrated
here, early mm-band observations can provide constraints on the density and
geometry of the CSM around SN progenitors on the same spatial scales as rapid
optical spectroscopy and X-ray observations, providing an independent
constraint with different modeling assumptions. We plan to continue
undertaking such observations for nearby CCSNe as part of POETS.

The Berger Time-Domain research group at Harvard is supported by the NSF and
NASA. T.E. is supported by NASA through the NASA Hubble Fellowship grant HST-
HF2-51504.001-A awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. The Submillimeter Array is a joint
project between the Smithsonian Astrophysical Observatory and the Academia
Sinica Institute of Astronomy and Astrophysics and is funded by the
Smithsonian Institution and the Academia Sinica. .

## References

  * Butler (2012) Butler, B. 2012, Atacama Large Millimeter/Submillimeter Array Memo 594, 594, 1
  * Chakraborti et al. (2013) Chakraborti, S., Petitpas, G., Zauderer, A., et al. 2013, The Astronomer’s Telegram, 4947, 1
  * Chevalier (1998) Chevalier, R. A. 1998, ApJ, 499, 810
  * Chevalier & Fransson (2017) Chevalier, R. A., & Fransson, C. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin
  * Gal-Yam et al. (2014) Gal-Yam, A., Arcavi, I., Ofek, E. O., et al. 2014, Nature, 509, 471
  * Gorosabel et al. (2010) Gorosabel, J., de Ugarte Postigo, A., Castro-Tirado, A. J., et al. 2010, A&A, 522, A14
  * Grefenstette et al. (2023) Grefenstette, B. W., Brightman, M., Earnshaw, H. P., Harrison, F. A., & Margutti, R. 2023, arXiv e-prints, arXiv:2306.04827
  * Hiramatsu et al. (2021) Hiramatsu, D., Howell, D. A., Van Dyk, S. D., et al. 2021, Nature Astronomy, 5, 903
  * Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
  * Horesh et al. (2013a) Horesh, A., Cao, Y., Mooley, K., & Carpenter, J. 2013a, The Astronomer’s Telegram, 5198, 1
  * Horesh et al. (2013b) Horesh, A., Stockdale, C., Fox, D. B., et al. 2013b, MNRAS, 436, 1258
  * Horesh et al. (2013c) Horesh, A., Stockdale, C., Fox, D. B., et al. 2013c, MNRAS, 436, 1258
  * Hosseinzadeh et al. (2023) Hosseinzadeh, G., Farah, J., Shrestha, M., et al. 2023, arXiv e-prints, arXiv:2306.06097
  * Itagaki (2023) Itagaki, K. 2023, Transient Name Server Discovery Report, 2023-1158, 1
  * Jacobson-Galan et al. (2023) Jacobson-Galan, W. V., Dessart, L., Margutti, R., et al. 2023, arXiv e-prints, arXiv:2306.04721
  * Kilpatrick et al. (2023) Kilpatrick, C. D., Foley, R. J., Jacobson-Galán, W. V., et al. 2023, arXiv e-prints, arXiv:2306.04722
  * Maeda et al. (2021) Maeda, K., Chandra, P., Matsuoka, T., et al. 2021, ApJ, 918, 34
  * Maeda et al. (2023) Maeda, K., Chandra, P., Moriya, T. J., et al. 2023, ApJ, 942, 17
  * Matsuoka et al. (2019) Matsuoka, T., Maeda, K., Lee, S.-H., & Yasuda, H. 2019, ApJ, 885, 41
  * Moriya et al. (2018) Moriya, T. J., Förster, F., Yoon, S.-C., Gräfener, G., & Blinnikov, S. I. 2018, MNRAS, 476, 2840
  * Nakar & Piro (2014) Nakar, E., & Piro, A. L. 2014, ApJ, 788, 193
  * Nayana et al. (2022) Nayana, A. J., Chandra, P., Krishna, A., & Anupama, G. C. 2022, ApJ, 934, 186
  * Neustadt et al. (2023) Neustadt, J. M. M., Kochanek, C. S., & Rizzo Smith, M. 2023, arXiv e-prints, arXiv:2306.06162
  * Perley et al. (2023) Perley, D. A., Gal-Yam, A., Irani, I., & Zimmerman, E. 2023, Transient Name Server AstroNote, 119, 1
  * Pledger & Shara (2023) Pledger, J. L., & Shara, M. M. 2023, arXiv e-prints, arXiv:2305.14447
  * Riess et al. (2022) Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7
  * Smith (2014) Smith, N. 2014, ARA&A, 52, 487
  * Smith et al. (2023) Smith, N., Pearson, J., Sand, D. J., et al. 2023, arXiv e-prints, arXiv:2306.07964
  * Soderberg et al. (2008) Soderberg, A. M., Berger, E., Page, K. L., et al. 2008, Nature, 453, 469
  * Weiler et al. (1986) Weiler, K. W., Sramek, R. A., Panagia, N., van der Hulst, J. M., & Salvati, M. 1986, ApJ, 301, 790
  * Weiler et al. (2007) Weiler, K. W., Williams, C. L., Panagia, N., et al. 2007, ApJ, 671, 1959
  * Yaron et al. (2017) Yaron, O., Perley, D. A., Gal-Yam, A., et al. 2017, Nature Physics, 13, 510
  * Zauderer et al. (2013) Zauderer, A., Soderberg, A., Chakraborti, S., et al. 2013, The Astronomer’s Telegram, 4946, 1
  * Zauderer et al. (2014) Zauderer, B. A., Kamble, A., Chakraborti, S., & Soderberg, A. 2014, The Astronomer’s Telegram, 5764, 1


# A New Paradigm for Fast and Repetitive Chirping of Alfvén Eigenmodes

Junyi Cheng  Beijing National Laboratory for Condensed Matter Physics and CAS
Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy
of Sciences, Beijing 100190, China  School of Physical Sciences, University of
Science and Technology of China, Hefei, Anhui 230026, China  School of
Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049,
China  Wenlu Zhang<EMAIL_ADDRESS>Beijing National Laboratory for
Condensed Matter Physics and CAS Key Laboratory of Soft Matter Physics,
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China  School
of Physical Sciences, University of Chinese Academy of Sciences, Beijing
100049, China  CAS Center for Excellence in Ultra-intense Laser Science,
Shanghai 201800, China  School of Physical Sciences, University of Science and
Technology of China, Hefei, Anhui 230026, China  Zhihong Lin  Department of
Physics and Astronomy, University of California, Irvine, California 92697, USA
Jian Bao  Beijing National Laboratory for Condensed Matter Physics and CAS Key
Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  School of Physical Sciences, University of
Chinese Academy of Sciences, Beijing 100049, China  Chao Dong  Beijing
National Laboratory for Condensed Matter Physics and CAS Key Laboratory of
Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences,
Beijing 100190, China  School of Physical Sciences, University of Chinese
Academy of Sciences, Beijing 100049, China  Jintao Cao  Beijing National
Laboratory for Condensed Matter Physics and CAS Key Laboratory of Soft Matter
Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Ding Li  Beijing National Laboratory for Condensed
Matter Physics and CAS Key Laboratory of Soft Matter Physics, Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  Songshan Lake
Materials Laboratory, Dongguan, Guangdong 523808, China  School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
CAS Center for Excellence in Ultra-intense Laser Science, Shanghai 201800,
China

###### Abstract

A novel 2D nonlinear dynamical paradigm is constructed to interpret the fast
and repetitive frequency chirping and amplitude oscillation of Alfvén
eigenmodes excited by energetic particles in fusion plasmas as observed in
global gyrokinetic simulations. In this non-perturbative paradigm of the
collisionless phase-space dynamics, wave-particle resonant interactions cause
the phase-space structure to continuously twist and fold, leading to the
repetitive excitation and decay of the Alfvén eigenmode. The radial
(perpendicular to the dominant wave-particle interaction) dependence of the
mode amplitude and toroidal precessional drifts of the energetic particles
leads to the 2D dynamics of wave-particle interactions, which is found to be
responsible for the repetitive process of formation and destruction of the
mode structure.

1. Introduction

Energetic particles, including energetic ions and electrons produced by the
fusion reaction and auxiliary heating, can excite various Alfvén eigenmodes in
magnetic confinement fusion plasmas such as ITER [1], which may induce
significant transport of energetic particles and degrade the overall plasma
confinement[2, 3, 4]. Increased energetic particle transport due to Alfvén
eigenmodes has been correlated [3] with a fast frequency oscillation
(chirping) with a sub-millisecond period (non-adiabatic regime[5]), which has
been observed in many tokamak experiments [6, 7, 8]. Previous theoretical
studies focus on the adiabatic regime (chirping time much longer than the wave
period) based on the 1D (in toroidal direction) dynamical model for the wave-
particle interactions[9], which has been investigated by numerical
simulations[10, 6, 11, 12, 13]. In this 1D dynamical model[9], the
relaxation/dissipation process[10] are essential to inducing the frequency
chirping.

However, the repetitive fast chirping of frequency and amplitude oscillation
of Alfvén eigenmodes excited by energetic particles in the toroidal geometry
have been observed by both first-principle kinetic [14] and hybrid
magnetohydrodynamic[15] simulations without source and sink, which is in
agreement with experiments and cannot be explained by the 1D dynamical model.
Here, we construct a novel 2D (in toroidal and radial directions)
collisionless nonlinear dynamical paradigm based on the global gyrokinetic
simulation result, which shows that the dynamics of fast and repetitive
frequency chirping and amplitude oscillation is a non-adiabatic[16], non-
perturbative and essentially intrinsic 2D nonlinear wave-particle resonant
interactions, where the evolution of 2D trajectories of resonant particles
restores unstable distribution function. In particular, the radial
(perpendicular to the dominant wave-particle interaction) dependence of the
mode amplitude and toroidal precessional drifts of the energetic particles
leads to the 2D dynamics of wave-particle interactions, which is found to be
responsible for the repetitive fast frequency chirping.

Specifically, the current simulations using the gyrokinetic toroidal code
(GTC) [17] find that beta-induced Alfvén eigenmodes (BAE) excited by energetic
electrons saturate and chirp, exhibiting various stages of nonlinear 2D
dynamics due to the radial dependence of the mode amplitudes and the toroidal
precessional drifts of energetic electrons. Our paradigm of the 2D nonlinear
wave-particle interactions can be applied to other Alfvén eigenmodes.

Figure 1: Time evolution of electrostatic potential $|e\phi/T_{e}|$ (red) and
the frequency $\omega$ (black) of dominant mode $(n,m)=(3,6)$ in upper panel
(a) and frequency power spectrum in lower panel (b). The color in panel (b)
represents the power intensity and the $95\%$ lines indicate the significance
level. The unit of the power intensity in the panel (b) is arbitrary. The mode
stage at different times is marked with letters A-K.

Figure 2: Poloidal contour plots of electrostatic potential $e\phi/T_{e}$. The
dotted circle is the $q=2$ rational surface. The $x$ axis is the major radius
$R/R_{0}$ and the $y$ axis is the vertical distance from the midplane
$Z/R_{0}$. Time steps $b-h$ correspond to $B-H$ in Fig. 1.

Figure 3: Radial profiles of electrostatic potential $e\phi/T_{e}$. The black,
blue and green lines are corresponding to the poloidal mode number $m=5,6,7$,
respectively. The dotted line is the $\psi$ value at the $q=2$ rational
surface. The $x$ axis is $\psi/\psi_{w}$ and the $y$ axis is $e\phi/T_{e}$.
Time steps $(b)-(h)$ are corresponding to $(B)-(H)$ in Fig. 1.

2. Simulation parameters

This electron BAE is chosen to simplify our 2D nonlinear dynamical model,
since only the toroidal precessional resonance of deeply trapped particles is
responsible[18] for the excitation. In our simulations of a tokamak with a
concentric circular cross-section, the safety factor profile is
$q=1.797+0.8(\psi/\psi_{w})-0.2(\psi/\psi_{w})^{2}$ so that the $q=2$ rational
surface is located at a minor radius $r_{s}=0.5a$, where $\psi$ is the
poloidal flux with $\psi=0$ on the axis and $\psi=\psi_{w}$ at the plasma
boundary, and $a$ is the tokamak minor radius at the wall. The inverse aspect
ratio is $\epsilon\equiv a/R_{0}=0.333$ in terms of on-axis major radius
$R_{0}$. The thermal plasma temperature is uniform with $T_{i}=T_{e}$, and the
on-axis beta is $\beta=4\pi n_{0}(T_{e}+T_{i})/B_{0}^{2}=0.00718$ with $B_{0}$
being the on-axis magnetic field. The thermal electron density $n_{0}$ is
uniform and the energetic electron density profile is
$n_{f}=0.05n_{0}(1.0+0.20(\tanh((0.24-\psi/\psi_{w})/0.06)-1.0))$, so that the
energetic electron density gradient peaks near the $q=2$ flux surface with
$R_{0}/L_{n_{f}}=15$, where $L_{n_{f}}\equiv-(d\ln n_{f}/dr)^{-1}$ is the
density gradient scale length of energetic electrons. The thermal ion density
$n_{0i}$ is obtained through the neutrality condition $n_{0i}=n_{0}+n_{f}$.
The energetic electrons are loaded with a local Maxwellian distribution with
the uniform temperature $T_{f}=25T_{e}$.

3. Fast and repetitive frequency chirping of BAE.

In the linear simulations[18], the most unstable toroidal mode number is $n=3$
with $k_{\theta}\rho_{i}=0.115$ at the $q=2$ rational surface, where the
device size is $a=104\rho_{i}$, $k_{\theta}=nq/r$ is the poloidal wave-vector,
$\rho_{i}=\sqrt{m_{i}T_{i}}/eB_{0}$ is the thermal ion gyro-radius, and
$m_{i}$ is the ion mass. Given the thermal plasma profiles, the frequency at
the BAE accumulation point is
$\omega_{0}=\sqrt{11T_{i}/2m_{i}R_{0}^{2}}\approx 2.34v_{i}/R_{0}$, where
$v_{i}=\sqrt{T_{i}/m_{i}}$ is the ion thermal velocity. In our collisionless
simulations, particle sources or sinks were turned off to allow free evolution
of temperature and density profiles.

The nonlinear time evolution of the electrostatic potential of the
$(n,m)=(3,6)$ mode is shown in Fig. 1, where the amplitude (red line in panel
a) and spectrum (panel b) of the electrostatic potential are measured at the
mode rational surface $r_{s}$. The time evolution of the power spectrum in the
lower panel is obtained through the wavelet analysis of the real part of the
electrostatic potential $\phi$. Fig. 1b shows a fast oscillation of the mode
intensity. For each time step in Fig. 1b, the frequency of the maximum power
intensity is plotted in Fig. 1a with a black curve, which shows that the
frequency starts with a linear frequency of $\omega_{BAE}$ and chirps in the
nonlinear stage roughly in phase with the intensity of the electrostatic
potential $\phi$ (red curve). In the linear stage, the measured real frequency
is $\omega_{BAE}=1.215\omega_{0}$ (the straight black line before point C) and
the linear growth rate $\gamma=0.0341\omega_{0}$, The mode saturates at point
D with the amplitude of the perturbed electrostatic potential being
$|e\phi/T_{e}|=6.84\times 10^{-3}$. After saturation (point D), the mode
evolves to a nonlinear chirping state. The first chirping period (points F-I)
is about seven wave periods, which is in a non-adiabatic regime[16]. There is
a small phase difference between the minima of amplitude and frequency. During
such chirping, the range of the frequency change is $0.157\omega_{0}$, which
is $12.9\%$ of the linear eigenfrequency $\omega_{BAE}$.

Fig. 2 and Fig. 3 shows the time evolution of the BAE mode structure of the
electrostatic potential on a poloidal plane and in the radial direction. The
mode structure is coherent and nearly symmetric with respect to the $q=2$
rational surface. It changes little from the linear stage (panel b) to the
nonlinear saturation with a maximal mode amplitude (panel d). It then moves
slightly outward as the mode amplitude decreases and the frequency starts
chirping down (panel e). At the lowest amplitude and frequency (panel f), the
mode structure is destroyed across the $q=2$ rational surface. Later, a
coherent mode structure forms but moves slightly inward as the mode amplitude
increases and the frequency chirps up (panel g). Finally, at the highest mode
amplitude and frequency (panel h), the mode structure becomes nearly symmetric
with respect to the $q=2$ rational surface, similar to that in the initial
nonlinear saturation (panel d). This time oscillation (from panel d to panel
h) for the mode amplitude, frequency, and mode structure repeats itself in the
next cycle (from panel h to panel j) but with a slightly smaller time period.

In controlled simulations, we find no frequency chirping when the nonlinearity
of the energetic electrons is removed. This indicates that energetic electrons
are exclusively responsible for such a frequency chirping process. The
saturation amplitude significantly increases when the nonlinearity of thermal
ions is removed, which suggests that thermal ions contribute to the nonlinear
saturation.

Figure 4: Time evolution of passive energetic-particle trajectories in
$(\delta\zeta,P_{\zeta})$ Lagrangian space. Different color represents initial
canonical angular momentum $P_{\zeta}$ as shown in panel (a). Dashed line is
the $q=2$ rational surface. Horizontal axis is the toroidal angle
$\delta\zeta/2\pi$ in the wave frame, and vertical axis is
$-P_{\zeta}/\psi_{w}$. The diagram of predicted particle trajectories
according to Eq. 1 are labeled in panel (b) with black dot line, respectively.
Labels $a-k$ correspond to $A-K$ in Fig. 1.

Figure 5: Evolution of passive energetic-particle perturbed distribution
function $\delta f(\delta\zeta,P_{\zeta})$ in Eulerian phase space. Different
color represents amplitude of perturbed distribution function $\delta f/f$.
Dashed line is the $q=2$ rational surface. Horizontal axis is the toroidal
angle $\delta\zeta/2\pi$ in the wave frame, and vertical axis is
$-P_{\zeta}/\psi_{w}$. The maximum of the absolute perturbed distribution
$|\delta f/f|$ for each time point is indicated. Labels $a-k$ correspond to
$A-K$ in Fig. 1.

4. Intrinsically 2D particle dynamics.

Given an energy $E$ and a magnetic moment $\mu$, the guiding-center motions of
the deeply trapped particles can be described by a pair of action-angle
variables $(\zeta,P_{\zeta})$ with $\zeta$ being the toroidal angle,
$P_{\zeta}=g\rho_{\|}-\psi$ being the normalized canonical angular momentum
[19], $2\pi g$ being the poloidal current and $\rho_{\|}=v_{\|}/B$ being the
normalized parallel velocity. Here, we examine the structure of the energetic
particle distribution function in the $(\zeta,P_{\zeta})$ phase space. Instead
of directly analyzing the data of energetic particles in the self-consistent
simulations, tracking orbits of passive particles will give a better
resolution and a lower particle noise for diagnosis by loading a large number
of the passive deeply-trapped particles in the resonant region of the phase
space.

These passive particles are uniformly initialized in the $(\psi,\zeta)$ plane
with $\theta=0$, $v_{\|}=-0.01\sqrt{T_{e}/m_{e}}$ and $\mu=32T_{e}/B_{0}$
within the resonant region. $P_{\zeta}/\psi_{w}$ are initialized in the range
of $(0.15,0.42)$. The particle positions in Fig. 4 are plotted in the
$(\zeta,P_{\zeta})$ space, where $P_{\zeta}$ values are approximately the
radial position since $P_{\zeta}\approx-\psi$ for the energetic electrons.

Using the nonlinear gyrokinetic equations [20, 21, 22, 23], we can calculate
the trajectory of the deeply trapped particles using the perturbed
electrostatic field $\phi=\hat{\phi}(\psi,\epsilon t)e^{i\alpha}$, where
$\epsilon$ indicates the slow evolution of the mode amplitude and
$\alpha=n\zeta-m\theta-\omega t$, $\alpha\in[0,2\pi)$, with a single toroidal
mode number $n$. Since the parallel wave wavelength is much longer than the
distance that a deeply trapped particle travels along the magnetic field line,
i.e., $k_{\|}\approx 0$, the guiding center equation of motion is

$\dot{\zeta}\propto\frac{\mu}{e}\frac{\partial
B_{0}}{\partial\psi}-\frac{\partial\hat{\phi}}{\partial\psi}e^{i\alpha},\qquad\dot{\psi}\propto
n\hat{\phi}e^{i(\alpha+\pi/2)},$ (1)

where $e$ is the elementary charge. Since in Eq. 1, the second term of
$\dot{\zeta}$ and $\dot{\psi}$ have a $\pi/2$ phase difference, these two
equations describe a Lissajous curve [24] rotating clockwise with a distortion
described by the first term of $\dot{\zeta}$. A 2D dynamical model can be
constructed using the reduced evolution equation for the perturbed
distribution function,

$\dot{\delta f}\propto nf_{0}\left(\frac{\mu}{e}\frac{\partial
B_{0}}{\partial\psi}+\frac{1}{f_{0}}\left.\frac{\partial
f_{0}}{\partial\psi}\right|_{v_{\perp}}\right)\hat{\phi}e^{i(\alpha+\pi/2)},$
(2)

which is closed with the field equations.

By using Taylor expansion around the $q=2$ rational surface and keeping the
first and second order terms of $B_{0}$ and eliminating the phase dependence
for $\dot{\zeta}$, and keeping the zero order terms of $\hat{\phi}$ for
$\dot{\psi}$, the equations are reduced to

$\dot{\delta\zeta}\propto\frac{\mu}{e}\frac{\partial^{2}B_{0}}{\partial\psi^{2}}\bigg{|}_{\psi_{0}}\delta\psi,\qquad\dot{\delta\psi}\propto
n\hat{\phi}(\psi_{0})e^{i(\alpha+\pi/2)},$ (3)

where $\delta\zeta=\zeta-\omega_{d0}t$, $\delta\psi=\psi-\psi_{0}$,
$\alpha=n\delta\zeta-m\theta-(\omega-n\omega_{d0})t$, the toroidal
precessional drift $\omega_{d}\propto(\mu/e)\partial
B_{0}/\partial\psi\propto\mu B_{0}/r$, $\omega_{d0}=\omega_{d}(\psi_{0})$ and
$\psi_{0}$ is the poloidal flux at the $q=2$ rational surface. Regarding
$(\delta\zeta,\ \delta\psi)$ as the generalized canonical coordinates $(x,\
v)$, Eq. 3 recovers the conventional 1D nonlinear dynamical model[25].

Compared with Eq. 1, Eq. 3 neglects the asymmetric radial dependence of
$B_{0}$ and $\hat{\phi}$ and the phase dependence for $\dot{\zeta}$. Thus, Eq.
1 generalizes the wave-particle resonant interactions to a non-perturbative
and 2D dynamical process, which incorporates the radial (perpendicular to the
dominant wave-particle interaction) dependence of the toroidal precessional
drifts as well as the mode structure and the phase dependence (the second term
of $\dot{\zeta}$ in Eq. 1). In particular, the asymmetric radial dependence is
essential to generate the nonlinear chirping phenomenon, which will be
discussed later.

Using Eq. 1, the energetic particle trajectories (Fig. 4) can be easily
understood. At time B, the trajectory starts twisting. The peak points ① and ②
in the panel (b) are the $\dot{\delta\psi}=0$ turning points, corresponding to
the $\alpha=0$ region. At time D, the symmetric phase space island is formed.
The points {③,④,⑥,⑧} in the panel (d) are the $\dot{\delta\psi}=0$ turning
points with the particle velocity only along the negative and positive
$\boldsymbol{\delta\zeta}$ directions, corresponding to the $\alpha=0$ region.
The points of ⑤ and ⑦ are $\dot{\delta\zeta}=0$ turning points with the
particle velocity only along the negative and positive
$\boldsymbol{\delta\psi}$ directions, corresponding to the $\alpha=3\pi/2+O$
and $\alpha=\pi/2-O$, where the corrections $O$ comes from the radial
dependence of the toroidal precessional drift.

This paradigm is first verified with a simple numerical test, where the
particle trajectories are modeled by Eq. 1 with all the parameters set to the
actual data in the simulations. The curves were shown in Fig. 4 (panel b, d,
black dot line), which are the instantaneous interface of particles initially
residing at the two sides of rational surface. Such a good agreement further
supports the validity of our paradigm.

5. Mode radial envelope plays an essential role.

The model of Eq. 1 is used to understand the nonlinear evolution of the
particle trajectories (Fig. 4). According to Eqs. 1 and 2, $\dot{\delta
f}\propto\dot{\psi}$. Thus, the perturbed distribution function (Fig. 5) of
the resonant particles is similar to the particle trajectories (Fig. 4). As
shown in Fig. 4b, particle trajectories are determined by the oscillations
along the $\boldsymbol{\delta\psi}$ and $\boldsymbol{\delta\zeta}$ directions
plus the motion in the negative and positive $\boldsymbol{\delta\zeta}$
directions for phases between $\alpha=(0,\pi)$ and $\alpha=(\pi,2\pi)$,
respectively. As the mode grows, the structure of particle trajectories
gradually starts twisting and generating the $\dot{\delta\psi}=0$ turning
points. This persistent twisting eventually leads to the generation of the
phase space island (Fig. 4 a→d), when the perturbed density (Fig. 5d) and the
electric field (Fig. 2d) both reach their maxima. Due to the asymmetric radial
dependence of the toroidal precessional drift, the trajectories in the phase
space island are asymmetric. At this moment, almost all the resonant particles
driving the instabilities are trapped in the phase space island.

Due to the asymmetric radial dependence of the toroidal precessional drift
along with the effect from the oscillating electric field, the island
continues its twisting (Fig. 4 e→f), which causes the outer part of the phase
island and the inner part of the neighboring island to twist and form a new
island. During this stage, the clumps and holes gradually line up vertically
in the radial direction so that the mode structure is destroyed and the
density perturbation becomes small at any toroidal angle. This leads to the
perturbed density (Fig. 5 e→f) and electric field (Fig. 2 e→f) decaying
gradually. Since $\dot{\delta\zeta}$ inside the $q=2$ rational surface is
larger than that outside (Fig. 4 e→f), the weight evolution inside is faster
than that outside (Fig. 5 e→f). Therefore, when the mode amplitude decreases,
the perturbed distribution of the outer part decreases more slowly than the
inner part, which leads to the mode structure moving outward radially (Fig.
2e). When the distribution function in the resonance region outside the
rational surface dominates the perturbed distribution function, the
precessional frequency $\omega_{d}$ $(\propto\mu B/r)$ of the resonant
particles decreases, corresponding to the reduction of the mode frequency.
This stage is characterized by a weakening BAE with the broken mode structure
and the down-chirping frequency.
6. Particle dynamics determines mode evolution.

When the BAE amplitude reaches its minimal level, the drive from the energetic
particles is restored and begins to excite the BAE instability again, and the
process of (f→i) proceeds as opposite to the process of (d→f). This leads to
the repetitive process of mode excitation and suppression accompanied by the
down-chirping in the frequency (f→i, i→k). This repetitive process corresponds
to the folding process of the phase space islands, which form the stretched-
nested-folding-layer structure (Fig. 4 f→i, i→k; Fig. 5 f→i, i→k). During this
folding process, the resonance region of the particle phase space and the mode
structure extend radially beyond the rational surface. Comparing with the 1D
Landau damping model[25], the radial envelope plays an important role.
Therefore, our paradigm reveals the intrinsically 2D dynamics in the
collisionless nonlinear Landau damping.

Our new paradigm of 2D wave-particle interaction can be extended to other
Alfvén eigenmodes in toroidal geometry with radial variations of mode
amplitude and radially asymmetric particle dynamics. For example, the reversed
magnetic field has the shear reversed point, so the evolution of the radial
dependence of the toroidal precessional drift is different from current
analysis, which can be the explanation for the RSAE’s [26, 27] chirping
behavior. And within our theoretical framework, considering additional physics
such as the sideband effect, ballooning structure, thermal particle kinetic
effect and so on, we can foresee abundant physics. Finally, our paradigm
suggestions, for the first time, that the perpendicular dependence of the the
wave-particle interaction plays an essential role for the nonlinear dynamics,
which could have implications beyond the fusion application.

Acknowledgments

Authors gratefully acknowledge useful discussions with Dr. W. Chen and Dr. H.
S. Zhang. This work was supported by the National MCF Energy R&D Program under
Grant Nos. 2018YFE0304100, 2018YEF0311300 and 2017YFE0301300; the the National
Natural Science Foundation of China under Grant Nos. 12025508, 11835016,
11875067, 11675256, 11675257, the Strategic Priority Research Program of
Chinese Academy of Sciences under Grant No. XDB16010300; the Key Research
Program of Frontier Science of Chinese Academy of Sciences under Grant No.
QYZDJ-SSW-SYS016; and SciDAC ISEP. Resources used in this research were
provided by the National Supercomputer Center in Tianjin (NSCC-TJ), the Oak
Ridge Leadership Computing Facility at Oak Ridge National Laboratory (OLCF)
and the National Energy Research Scientific Computing Center (NERSC).

## References

  * [1] K. Ikeda, Nucl. Fusion 47, E01 (2007).
  * [2] E. J. Strait, W. W. Heidbrink, A. D. Turnbull, M. S. Chu, and H. H. Duong, Nucl. Fusion 33, 1849 (1993).
  * [3] M. Podestà, R. E. Bell, N. A. Crocker, E. D. Fredrickson, N. N. Gorelenkov, W. W. Heidbrink, S. Kubota, B. P. LeBlanc, and H. Yuh, Nucl. Fusion 51, 063035 (2011).
  * [4] M. Garcia-Munoz, I. G. J. Classen, B. Geiger, W. W. Heidbrink, M. A. Van Zeeland, S. Äkäslompolo, R. Bilato, V. Bobkov, M. Brambilla, G. D. Conway, S. da Graça, V. Igochine, Ph. Lauber, N. Luhmann, M. Maraschek, F. Meo, H. Park, M. Schneller, G. Tardini, and the ASDEX Upgrade Team, Nucl. Fusion 51, 103013 (2011).
  * [5] L. Chen and F. Zonca, Rev. Mod. Phys. 88, 015008 (2016).
  * [6] S. D. Pinches, H. L. Berk, M. P. Gryaznevich, S. E. Sharapov, and JET-EFDA Contributors, Plasma Phys. Control. Fusion 46, S47 (2004).
  * [7] M. P. Gryanzevich and S. E. Sharapov, Nucl. Fusion 46, S942 (2006).
  * [8] W. W. Heidbrink, E. Ruskov, E. D. Fredrickson, N. Gorelenkov, S. S. Medley, H. L. Berk, and R. W. Harvey, Plasma Phys. Control. Fusion 48, 1347 (2006).
  * [9] H. L. Berk, B. N. Breizman, and M. Pekker, Phys. Rev. Lett. 76, 1256 (1996).
  * [10] M. K. Lilley, B. N. Breizman, and S. E. Sharapov, Phys. Rev. Lett. 102, 195003 (2009).
  * [11] J. Candy, H. L. Berk, B. N. Breizman, and F. Porcelli, Phys. Plasmas 6, 1822 (1999).
  * [12] R. G. L. Vann, H.L. Berk, and A.R. Soto-Chavez, Phys. Rev. Lett. 99, 025003 (2007).
  * [13] H. Hezaveh, Z. S. Qu, B. N. Breizman, and M. J. Hole, Nucl. Fusion 60, 056014 (2020).
  * [14] H. S. Zhang, Z. Lin, and I. Holod, Phys. Rev. Lett. 109, 025001 (2012).
  * [15] X. Wang, S. Briguglio, L. Chen, C. Di Troia, G. Fogaccia, G. Vlad, and F. Zonca, Phys. Rev. E 86, 045401(R) (2012).
  * [16] F. Zonca, L. Chen, S. Briguglio, G. Fogaccia, G. Vlad, and X. Wang, New J. Phys. 17, 013052 (2015).
  * [17] Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang, and R. B. White, Science 281, 1835 (1998).
  * [18] J. Cheng, W. Zhang, Z. Lin, I. Holod, D. Li, Y. Chen, and J. Cao, Phys. Plasmas 23, 052504 (2016).
  * [19] R. B. White, The theory of toroidally confined plasmas (World Scientific Pusblishing Co Inc. 2013).
  * [20] A. J. Brizard and T. S. Hahm, Rev. Mod. Phys. 79, 421 (2007).
  * [21] P. Liu, W. Zhang, C. Dong, J. Lin, Z. Lin, and J. Cao, Nucl. Fusion 57, 126011 (2017).
  * [22] P. Liu, W. Zhang, C. Dong, J. Lin, Z. Lin, J. Cao, and D. Li, Phys. Plasmas 24, 112114 (2017).
  * [23] I. Holod, W. L. Zhang, Y. Xiao, and Z. Lin, Phys. Plasmas 16, 122307 (2009).
  * [24] T. B. Greenslade Jr., Phys. Teach. 31, 364 (1993).
  * [25] T. M. O’Neil, J. H. Winfrey, and J. H. Malmberg, Phys. Fluids 14, 1204 (1971).
  * [26] H. L. Berk, D. N. Borba, B. N. Breizman, S. D. Pinches, and S. E. Sharapov, Phys. Rev. Lett. 87, 185002 (2001).
  * [27] S. D. Pinches, H. L. Berk, D. N. Borba, B. N. Breizman, S. Brguglio, A. Fasoli, G. Fogaccia, M. P. Gryanzevich, V. Kiptily, M. J. Mantsinen, S. E. Sharapov, D. Testa, R. G. L. Vann, G. Vlad, F. Zonca, and JET-EFDA Contributors, Plasma Phys. Control. Fusion 46, B187 (2004).


* Schmidt and Enss [2011] R. Schmidt and T. Enss, Excitation spectra and rf response near the polaron-to-molecule transition from the functional renormalization group, Phys. Rev. A 83, 063620 (2011).
  * Schmidt  _et al._ [2012] R. Schmidt, T. Enss, V. Pietilä, and E. Demler, Fermi polarons in two dimensions, Phys. Rev. A 85, 021602 (2012).
  * Kohstall  _et al._ [2012] C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. M. Bruun, F. Schreck, and R. Grimm, Metastability and coherence of repulsive polarons in a strongly interacting Fermi mixture, Nature 485, 615 (2012).
  * Ngampruetikorn  _et al._ [2012] V. Ngampruetikorn, J. Levinsen, and M. M. Parish, Repulsive polarons in two-dimensional Fermi gases, Eur. Phys. Let. 98, 30005 (2012).
  * Zhang  _et al._ [2012] Y. Zhang, W. Ong, I. Arakelyan, and J. E. Thomas, Polaron-to-Polaron transitions in the radio-frequency spectrum of a quasi-two-dimensional Fermi gas, Phys. Rev. Lett. 108, 235302 (2012).
  * Catani  _et al._ [2012] J. Catani, G. Lamporesi, D. Naik, M. Gring, M. Inguscio, F. Minardi, A. Kantian, and T. Giamarchi, Quantum dynamics of impurities in a one-dimensional bose gas, Phys. Rev. A 85, 023623 (2012).
  * Fukuhara  _et al._ [2013] T. Fukuhara, A. Kantian, M. Endres, M. Cheneau, P. Schauß, S. Hild, D. Bellem, U. Schollwöck, T. Giamarchi, C. Gross, I. Bloch, and S. Kuhr, Quantum dynamics of a mobile spin impurity, Nat. Phys. 9, 235 (2013).
  * Burovski  _et al._ [2014] E. Burovski, V. Cheianov, O. Gamayun, and O. Lychkovskiy, Momentum relaxation of a mobile impurity in a one-dimensional quantum gas, Phys. Rev. A 89, 041601 (2014).
  * Ardila and Giorgini [2016] L. A. P. n. Ardila and S. Giorgini, Bose polaron problem: Effect of mass imbalance on binding energy, Phys. Rev. A 94, 063640 (2016).
  * Schmidt  _et al._ [2016] R. Schmidt, H. R. Sadeghpour, and E. Demler, Mesoscopic Rydberg impurity in an atomic quantum gas, Phys. Rev. Lett. 116, 105302 (2016).
  * Grusdt  _et al._ [2017] F. Grusdt, R. Schmidt, Y. E. Shchadilova, and E. Demler, Strong-coupling Bose polarons in a Bose-Einstein condensate, Phys. Rev. A 96, 013607 (2017).
  * Volosniev and Hammer [2017] A. G. Volosniev and H.-W. Hammer, Analytical approach to the bose-polaron problem in one dimension, Phys. Rev. A 96, 031601 (2017).
  * Scazza  _et al._ [2017] F. Scazza, G. Valtolina, P. Massignan, A. Recati, A. Amico, A. Burchianti, C. Fort, M. Inguscio, M. Zaccanti, and G. Roati, Repulsive Fermi polarons in a resonant mixture of ultracold ${}^{6}\mathrm{Li}$ atoms, Phys. Rev. Lett. 118, 083602 (2017).
  * Gamayun  _et al._ [2018] O. Gamayun, O. Lychkovskiy, E. Burovski, M. Malcomson, V. V. Cheianov, and M. B. Zvonarev, Impact of the injection protocol on an impurity’s stationary state, Phys. Rev. Lett. 120, 220605 (2018).
  * Guenther  _et al._ [2018] N.-E. Guenther, P. Massignan, M. Lewenstein, and G. M. Bruun, Bose polarons at finite temperature and strong coupling, Phys. Rev. Lett. 120, 050405 (2018).
  * Schmidt  _et al._ [2018] R. Schmidt, M. Knap, D. A. Ivanov, J.-S. You, M. Cetina, and E. Demler, Universal many-body response of heavy impurities coupled to a Fermi sea: a review of recent progress, Rep. Prog. Phys. 81, 024401 (2018).
  * Mistakidis  _et al._ [2019a] S. I. Mistakidis, G. C. Katsimiga, G. M. Koutentakis, and P. Schmelcher, Repulsive Fermi polarons and their induced interactions in binary mixtures of ultracold atoms, New J. Phys. 21, 043032 (2019a).
  * Mistakidis  _et al._ [2019b] S. I. Mistakidis, G. C. Katsimiga, G. M. Koutentakis, T. Busch, and P. Schmelcher, Quench dynamics and orthogonality catastrophe of Bose polarons, Phys. Rev. Lett. 122, 183001 (2019b).
  * Mistakidis  _et al._ [2021] S. I. Mistakidis, G. M. Koutentakis, F. Grusdt, H. R. Sadeghpour, and P. Schmelcher, Radiofrequency spectroscopy of one-dimensional trapped Bose polarons: crossover from the adiabatic to the diabatic regime, New J. Phys. 23, 043051 (2021).
  * Taglieber  _et al._ [2008] M. Taglieber, A.-C. Voigt, T. Aoki, T. W. Hänsch, and K. Dieckmann, Quantum degenerate two-species Fermi-Fermi mixture coexisting with a Bose-Einstein condensate, Phys. Rev. Lett. 100, 010401 (2008).
  * Tiecke  _et al._ [2010] T. G. Tiecke, M. R. Goosen, A. Ludewig, S. D. Gensemer, S. Kraft, S. Kokkelmans, and J. T. M. Walraven, Broad feshbach resonance in the ${}^{6}\mathrm{Li}\mathrm{\text{$-$}}^{40}\mathrm{K}$ mixture, Phys. Rev. Lett. 104, 053202 (2010).
  * Naik  _et al._ [2011] D. Naik, A. Trenkwalder, C. Kohstall, F. M. Spiegelhalder, M. Zaccanti, G. Hendl, F. Schreck, R. Grimm, T. M. Hanna, and P. S. Julienne, Feshbach resonances in the 6Li-40K Fermi-Fermi mixture: elastic versus inelastic interactions, Europhys. J. D 65, 55 (2011).
  * Cetina  _et al._ [2016] M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. M. Walraven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt, M. Knap, and E. Demler, Ultrafast many-body interferometry of impurities coupled to a Fermi sea, Science 354, 96 (2016).
  * Voigt  _et al._ [2009] A.-C. Voigt, M. Taglieber, L. Costa, T. Aoki, W. Wieser, T. W. Hänsch, and K. Dieckmann, Ultracold heteronuclear Fermi-Fermi molecules, Phys. Rev. Lett. 102, 020405 (2009).
  * Serwane  _et al._ [2011] F. Serwane, G. Zürn, T. Lompe, T. B. Ottenstein, A. N. Wenz, and S. Jochim, Deterministic preparation of a tunable few-fermion system, Science 332, 336 (2011).
  * Wenz  _et al._ [2013] A. Wenz, G. Zürn, S. Murmann, I. Brouzos, T. Lompe, and S. Jochim, From few to many: Observing the formation of a Fermi sea one atom at a time, Science 342, 457 (2013).
  * Krönke  _et al._ [2013] S. Krönke, L. Cao, O. Vendrell, and P. Schmelcher, Non-equilibrium quantum dynamics of ultra-cold atomic mixtures: the multi-layer multi-configuration time-dependent Hartree method for bosons, New J. Phys. 15, 063018 (2013).
  * Cao  _et al._ [2013] L. Cao, S. Krönke, O. Vendrell, and P. Schmelcher, The multi-layer multi-configuration time-dependent Hartree method for bosons: Theory, implementation, and applications, Chem. Phys. 139, 134103 (2013).
  * Cao  _et al._ [2017] L. Cao, V. Bolsinger, S. I. Mistakidis, G. M. Koutentakis, S. Krönke, J. M. Schurer, and P. Schmelcher, A unified ab initio approach to the correlated quantum dynamics of ultracold fermionic and bosonic mixtures, J. Chem. Phys. 147, 044106 (2017).
  * Hummel  _et al._ [2021] F. Hummel, M. T. Eiles, and P. Schmelcher, Synthetic dimension-induced conical intersections in Rydberg molecules, Phys. Rev. Lett. 127, 023003 (2021).
  * Levi  _et al._ [2011] L. Levi, M. Rechtsman, B. Freedman, T. Schwartz, O. Manela, and M. Segev, Disorder-enhanced transport in photonic quasicrystals, Science 332, 1541 (2011).
  * Kolkowitz  _et al._ [2016] S. Kolkowitz, S. L. Bromley, T. Bothwell, M. L. Wall, G. E. Marti, A. P. Koller, X. Zhang, A. M. Rey, and J. Ye, Spin–orbit-coupled fermions in an optical lattice clock, Nature 542, 66 (2016).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * Celi  _et al._ [2014] A. Celi, P. Massignan, J. Ruseckas, N. Goldman, I. B. Spielman, G. Juzeliūnas, and M. Lewenstein, Synthetic gauge fields in synthetic dimensions, Phys. Rev. Lett. 112, 043001 (2014).
  * Boada  _et al._ [2012] O. Boada, A. Celi, J. I. Latorre, and M. Lewenstein, Quantum simulation of an extra dimension, Phys. Rev. Lett. 108, 133001 (2012).
  * Olshanii [1998] M. Olshanii, Atomic scattering in the presence of an external confinement and a gas of impenetrable bosons, Phys. Rev. Lett. 81, 938 (1998).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010).
  * Born and Oppenheimer [1927] M. Born and R. Oppenheimer, Zur Quantentheorie der Moleküle, Ann. Phys. 389, 457 (1927).
  * Schuster  _et al._ [2012] T. Schuster, R. Scelle, A. Trautmann, S. Knoop, M. K. Oberthaler, M. M. Haverhals, M. R. Goosen, S. Kokkelmans, and E. Tiemann, Feshbach spectroscopy and scattering properties of ultracold Li + Na mixtures, Phys. Rev. A 85, 042721 (2012).
  * Dirac [1930] P. A. M. Dirac, On the annihilation of electrons and protons, Proc. Cambridge Phil. Soc. 26, 361–375 (1930).
  * Frenkel [1934] J. Frenkel,  _Wave mechanics, advanced general theory_ , Vol. 1 (Oxford University Press, 1934).
  * Köppel  _et al._ [1984] H. Köppel, W. Domcke, and L. S. Cederbaum, Multimode molecular dynamics beyond the Born-Oppenheimer approximation, in  _Advances in Chemical Physics_ (John Wiley & Sons, Ltd, 1984) pp. 59–246.
  * Pacher  _et al._ [1989] T. Pacher, C. A. Mead, L. S. Cederbaum, and H. Köppel, Gauge theory and quasidiabatic states in molecular physics, Chem. Phys. 91, 7057 (1989).
  * Epstein [1966] S. T. Epstein, Ground-state energy of a molecule in the adiabatic approximation, Chem. Phys. 44, 836 (1966).
  * Baer [2006] M. Baer,  _Beyond Born-Oppenheimer_ , 1st ed. (John Wiley & Sons, Hoboken, New Jersey, 2006).
  * Gurari [1953] M. Gurari, XXXVI. Self energy of slow electrons in polar materials, London, Edinburgh, Dublin Philos. Mag. J. Sci. 44, 329 (1953).
  * Lee and Pines [1952] T.-D. Lee and D. Pines, The motion of slow electrons in polar crystals, Phys. Rev. 88, 960 (1952).
  * Lee  _et al._ [1953] T. D. Lee, F. E. Low, and D. Pines, The motion of slow electrons in a polar crystal, Phys. Rev. 90, 297 (1953).
  * Busch  _et al._ [1998] T. Busch, B.-G. Englert, K. Rzażewski, and M. Wilkens, Two cold atoms in a harmonic trap, Found. Phys. 28, 549 (1998).
  * Farrell and van Zyl [2009] A. Farrell and B. P. van Zyl, Universality of the energy spectrum for two interacting harmonically trapped ultra-cold atoms in one and two dimensions, J. Phys. A: Math. Theor. 43, 015302 (2009).
  * Budewig  _et al._ [2019] L. Budewig, S. I. Mistakidis, and P. Schmelcher, Quench dynamics of two one-dimensional harmonically trapped bosons bridging attraction and repulsion, Mol. Phys. 117, 2043 (2019).
  * Domcke  _et al._ [2004] W. Domcke, D. Yarkony, and H. Köppel,  _Conical Intersection: Electronic Structure, Dynamics & Spectroscopy_ (World Scientific, 2004).
  * Kemper  _et al._ [1977] M. Kemper, J. Van Dijk, and H. Buck, On the comparison between crude and adiabatic Born-Oppenheimer coupling elements, Chem. Phys. Lett. 48, 590 (1977).
  * Anderson [1967a] P. W. Anderson, Infrared catastrophe in Fermi gases with local scattering potentials, Phys. Rev. Lett. 18, 1049 (1967a).
  * Anderson [1967b] P. W. Anderson, Ground state of a magnetic impurity in a metal, Phys. Rev. 164, 352 (1967b).
  * Thouless [1961] D. J. Thouless,  _The Quantum Mechanics of Many-Body Systems_, Pure and Applied Physics, Vol. 11 (Elsevier, London, 1961) pp. 153 – 162.
  * Messiah [1981] A. Messiah,  _Quantum Mechanics_ (Elsevier Science, Amsterdam, 1981).
  * Bouvrie  _et al._ [2014] P. A. Bouvrie, A. P. Majtey, M. C. Tichy, J. S. Dehesa, and A. R. Plastino, Entanglement and the Born-Oppenheimer approximation in an exactly solvable quantum many-body system, Eur. Phys. J. D 68, 346 (2014).
  * Izmaylov and Franco [2017] A. F. Izmaylov and I. Franco, Entanglement in the Born–Oppenheimer approximation, J. Chem. Theory Comput. 13, 20 (2017).
  * Streltsov  _et al._ [2010] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, General mapping for bosonic and fermionic operators in Fock space, Phys. Rev. A 81, 022124 (2010).
  * Girardeau [1960] M. Girardeau, Relationship between systems of impenetrable bosons and fermions in one dimension, J. Math. Phys. 1, 516 (1960).