\section{INTRODUCTION}

Planetesimals are $\sim1-1000$ km sized bodies in protoplanetary disks (PPDs) as building blocks of planets.
Following the growth of dust grains to mm-cm sizes in PPDs, planetesimals are believed to form from such grown dust particles, representing the intermediate stage of planet formation. This stage is also considered to be the most difficult, as surface forces can no longer stick particles together upon collisions, nor does self-gravity which dominates among solids of larger sizes \citep{chiang10}. Generally, planetesimal formation is understood as a two-phase process, where particles first concentrate into high-density clumps, followed by their gravitational collapse. What drives the formation of particle clumps thus becomes the key to understanding planetesimal formation. Regardless of the specific mechanism, particles of mm-cm sizes strongly interact with gas via aerodynamic drag, and the gas dynamics in PPDs, especially disk turbulence, must play an important role.

The leading theory of planetesimal formation is the streaming instability \citep[SI,][]{youdin05}, resulting from the two-way drag force between gas and dust in the presence of a background pressure gradient in PPDs. Hybrid gas-dust simulations found that the non-linear evolution of the SI leads to dust clumping \citep[e.g.,][]{johansen09b,bai10b,carrera15,yang17}, provided the dust abundance $Z$ (ratio of dust to gas surface density) exceeds some threshold (typically $\gtrsim0.02$, but see \citealp{li21}). While highly successful, vast majority of SI simulations assumed no external turbulence. However, recent studies found that even modest level of external disk turbulence suppresses the SI \citep{umurhan20,chen20,gole20}, challenging the SI paradigm.

\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\textwidth]{Figures/ppII_intro.pdf}
\end{center}
\caption{
A snapshot of 3D visualization of the dust spatial distribution in our simulation with strong gas bump (run Z2H) at the time of $T=1380\Omega_K^{-1}$ after inserting dust particles. Dust density is mapped to brightness, where brighter region indicates higher dust density. Also shown in the figure are vertically averaged gas density projected to the $x$-$y$ ($R$-$\phi$) plane, and vertically and azimuthally averaged radial profiles of dust and gas density, with amplitudes of the profiles re-scaled for better illustration. \label{fig:intro}}
\end{figure*}

The outer disk regions ($\gtrsim20$ AU) are the most accessible to spatially-resolved disk observations. At sub-mm (where mm-sized dust radiates most efficiently), the prevalence of ring-like substructures has recently been firmly established \citep[e.g.,][]{andrews18,long18}, widely believed to be indicative of dust trapping in pressure bumps \citep{dullemond18}\footnote{Unless otherwise noted, ``pressure bump" refers to dust trapping pressure bumps in this paper, i.e., strong enough to overcome the background pressure gradient and generate local pressure maxima.}. Dust trapping in local pressure maximum is firstly demonstrated by \cite{whipple72}, and the formation mechanism of pressure bumps can be planet-disk interaction \citep[see review of][]{paardekooper22} or non-planet origins such as MHD zonal flows \citep[e.g.,][]{suriano18,riols20,cui21}. While there are direct and indirect evidence pointing to weak disk turbulence \citep{pinte16,flaherty17,flaherty18,teague18}, modeling dust ring widths suggest modest level of turbulence \citep{dullemond18,rosotti20}. Together, these observational constraints lead to the emerging picture that most mm-sized dust particles in outer PPDs are most likely concentrated in weakly-to-modestly turbulent, largely-axisymmetric gas pressure bumps. The strongly enhanced local dust abundance in such pressure bumps yields favorable condition for planetesimal formation. However, given that the SI is unexpected to operate without pressure gradient, and can be suppressed by external turbulence, it is unclear whether it will further lead to dust clumping and hence planetesimal formation.

To address the feasibility of planetesimal formation in more realistic environment of turbulent dust rings, we conduct hybrid gas-particle simulations in the standard shearing-box framework that mimics a local patch of the disk. It is anticipated that the magneto-rotational instability (MRI) is one primary mechanism for driving turbulence in the outer PPDs. Due to very low level of disk ionization, the MRI is weakened by ambipolar diffusion (AD) as the dominant non-ideal magnetohydrodynamic (MHD) effect in the outer disk. Turbulence in our simulations are self-consistently generated by the MRI with AD. We further implemented a novel forcing prescription to drive a pressure bump in the center of our simulation box with essentially no impact on disk turbulence.
Our simulations generalize from our earlier work \citep[][hereafter \citetalias{xu22}]{xu22} by imposing and enforcing the pressure bump, allowing us to achieve a steady-state configuration with dust trapping towards pressure bumps balanced by dust diffusion in realistic disk turbulence in a self-consistent and well-controlled manner. They differ from the recent works of \citep[][]{carrera21,carrera22a,carrera22b}, who conducted simulations with a similar procedure to enforce the pressure bump but without external turbulence, thus cannot achieve such steady-state balance of dust trapping and diffusion as envisioned in recent disk observations \citep{dullemond18,rosotti20}.

This letter is organized as follows.
We first introduce the implementation of the pressure bump and simulation setup in \S \ref{sec:method}, then present our main results in \S \ref{sec:results}. In \S \ref{sec:discu}, we will further analyze the results, focusing on dust clumping conditions and ring dynamics, followed by comparison with recent studies and discussion on observational implications. We summarize and discuss future prospects in \S \ref{sec:conclusions}.

\begin{table}[t]
    \centering
     \caption{List of simulation runs.}
    \label{tab:sumsetup}
    \begin{tabular}{lcccccccc}
   Runs  & $Z$ & $N_x \times N_y \times N_z$  & $A_{bump}$ & $H_d/H$ & $w_{d}/H$ \\
   \hline
   Z0 & 1e-30 & 512 $\times$ 512 $\times$ 64  & 0.5 & 0.067 & 0.55\\
   Z2  & 0.02 & 512 $\times$ 512 $\times$ 64  & 0.5 & 0.014 &0.16\\
   Z0H & 1e-30 & 1024 $\times$ 1024 $\times$ 128  & 0.5 & -&-\\
   Z1H  & 0.01 & 1024 $\times$ 1024 $\times$ 128  & 0.5 & 0.025 &0.25\\
   Z2H  & 0.02 & 1024 $\times$ 1024 $\times$ 128 & 0.5 & 0.017 &0.16\\
      \hline
   Z0w & 1e-30 & 512 $\times$ 512 $\times$ 64 & 0.25 & 0.075 &0.83\\
   Z1w  & 0.01 & 512 $\times$ 512 $\times$ 64 & 0.25 & 0.031 &0.87\\
   Z2w  & 0.02 & 512 $\times$ 512 $\times$ 64 & 0.25 & 0.022 &0.74\\
   \hline
    \end{tabular}
    \begin{tablenotes}
    \item All simulations have a boxsize of $8H \times 16H \times H$, with $\beta_0=12800$, $Am=2$. Gas bump width  $w_{bump} = 2H$, with no background pressure gradient. Particle number is 2 per cell on average, with particle stopping time $\tau_s=0.1$.
    \item All parameters in this table are initial setups, except that $H_d$ and $w_{d}$ are measured from simulation.
    \end{tablenotes}
\end{table}

\section{Method and Simulations} \label{sec:method}

We use the Athena MHD code \citep{stone08} to perform 3D hybrid gas-particle local shearing-box simulations considering non-ideal MHD with ambipolar diffusion.
The formulation and basic simulation setup is similar to \citetalias{xu22}, which we only describe briefly here. The $x, y, z$ coordinates in a local shearing-box correspond to radial, azimuthal and vertical directions in a global disk, as illustrated in Figure \ref{fig:intro}. We use an isothermal equation of state with sound speed $c_s=1$ and local Keplerian frequency $\Omega_K=1$ in code units, which define our length unit of disk scale height $H=c_s/\Omega_K=1$. Gas density, velocity and pressure are denoted by $\rho, {\boldsymbol v}$ and
$P=\rho c_s^2$, respectively.
Dust particles are characterized by their dimensionless stopping time $\tau_s\equiv t_s\Omega_0$, and the height-integrate dust-to-gas mass ratio $Z$ for background gas surface density. Our simulations focus on the midplane region, and hence we neglect vertical gravity in the gas (unstratified), but vertical gravity is included for particles to allow dust settling. We impose a net vertical magnetic flux characterized by $\beta_0=12800$, the ratio of background gas pressure to magnetic pressure of the net vertical field. Ambipolar diffusion is characterized by the Ambipolar Elsasser number $Am$, where we fix $Am=2$ appropriate for typical outer disk conditions.

\subsection{Implementation of the gas pressure bump} \label{subsec:formulation}

On top of the simulation setup in \citetalias{xu22}, we further implement a gas pressure bump by a prescribed forcing procedure.
Let $\rho_0=1$ be the constant background gas density. The forcing will redistribute the density, and we aim to achieve a Gaussian-like density (and hence pressure) bump in our simulation box of the following form
\begin{equation}
\frac{\rho(x)}{\rho_0}= 1+ A_{\rm bump}\exp\bigg[-\frac{x^2}{2w_{\rm bump}^2}-0.2\bigg(\frac{x^2}{2w_{\rm bump}^2}\bigg)^4\bigg],\\
\end{equation}
where
$A_{\rm bump}$ is the dimensionless amplitude of the pressure bump, $w_{\rm bump}$ is the effective width of the bump. The first term in the exponential gives the overall Gaussian profile, while the second term damps the profile towards the domain edges to satisfy the shearing-periodic radial boundary condition.

The radial pressure gradient associated with this bump must be balanced by an azimuthal flow:
\begin{equation}
    \frac{v'_{y,{\rm bump}}(x)}{c_s} =
    \frac{H}{2}\frac{d\ln\rho(x)}{dx}\ ,
\end{equation}
where the prime $'$ in $v_y$ denotes deviation from the Keplerian velocity $v_K=-(3/2)\Omega_Kx$.
Note that in our simulations, as we focus on regions near the pressure maxima, we do not impose any background pressure gradient.

To generate and sustain the bump,
we apply an azimuthal forcing in our simulation, so that the azimuthal velocity relax to the desired profiles within a timescale of $t_{\rm relax}=60\Omega_K^{-1}$:
\begin{equation}
\Delta v'_y(x) = [v_{y,{\rm bump}}(x) - \overline{v'_y}(x)] [1-\exp(-\Delta t/t_{\rm relax})],
\end{equation}
where $\overline{v'_y}(x)$ is azimuthally and vertically averaged. Our forcing procedure is analogous to that of \cite{carrera21}, but with several differences. First, our forcing acts on azimuthal velocity (momentum) only, thus is mass conserving. It directly mimics the torque exerted on the disk responsible for making the bump
\citep[be it of planet origin, or non-planet origin as in][]{cui21},
and the desired density profile is a direct consequence from this forcing. Second, our forcing is homogeneously applied in the $y$ and $z$ directions, as opposed to relaxation {\it at individual grid cells} which introduces local friction. Third, our relaxation time is much longer, as opposed to $\sim1\Omega_K^{-1}$. These facts minimizes the impact of our forcing procedure to disk turbulence, especially avoiding artificial damping at small scales.\footnote{Our choice of relaxation time is generally consistent with the typical timescale for magnetic flux evolution (and thus bump evolution) indicated in \cite{cui21}. We further tested $t_{\rm relax}=30\Omega_K^{-1}$ and $120\Omega_K^{-1}$ respectively (not shown in paper), and the result is not sensitive to the change of $t_{\rm relax}$ in our simulation in this paper. }

\subsection{Simulation Procedures and Parameters} \label{subsec:simsetup}

The list of our simulation runs are shown in Table \ref{tab:sumsetup}.
We explore strengths of the gas pressure bump, presence/absence of dust feedback, different levels of solid abundance, and simulation resolutions. The simulation setup and parameters are described as follows.

We first run simulations without dust particles for a time of $t_0 = 120 \Omega_K^{-1}$ for the MRI to grow and fully saturate into turbulence, following our previous work (\citetalias{xu22}).
The gas pressure bump is initiated at the beginning of the simulation, with bump width $w_{\rm bump} = 2H$.
We explore two cases: strong and weak bumps, with bump amplitude $A_{\rm bump} = 0.5$ and $0.25$ separately, and the  weak bump runs are marked by ``w'' in the run name.\footnote{We have examined that these bump amplitudes are below the threshold for triggering the Rossby wave instability (RWI) for all azimuthal wave numbers under global settings \citep{ono16}, and we anticipate the presence of the MRI turbulence also makes our pressure bump less prone to the RWI.} Our choice of bump width generally agrees with observational constraints \citep{dullemond18,rosotti20}, and both the bump width and amplitude in our simulation are reasonably consistent with the results of \cite{cui21}.

\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\textwidth]{Figures/hist_gbumpw.pdf}
\end{center}
\caption{
Time evolution of dust maximum density (upper panel), dust ring width (middle panel), and dust scale height (lower panel) for all our simulation runs. Dust maximum density for no-feedback runs (Z0, Z0H, and Z0w) are scaled to $Z=0.02$. \label{fig:hist}}
\end{figure*}

Dust particles are then inserted into the simulation box at the time of $t_0 = 120 \Omega_K^{-1}$ without backreaction to the gas, with random horizontal positions and a vertical distribution following a Gaussian profile $\propto \exp(-z^2/2H_{\rm d0}^2)$, with $H_{\rm d0} = 0.02 H$. We define $T=T_{\rm sim}-t_0$, where $T_{\rm sim}$ is the actual simulation time, and use $T$ instead of $T_{\rm sim}$ throughout this paper.
We use a single particle species with fixed stopping time $\tau_s=0.1$,
corresponding to mm-cm size particles in the midplane regions of the outer disk in standard disk models \citep[e.g.,][]{chiang10}.
As dust particles passively follow the gas,
settling balanced by turbulent diffusion in the vertical direction makes them maintain certain scale height (see \citetalias{xu22}), while the balance between trapping toward pressure maxima and turbulent diffusion in the radial direction allow them to concentrate into a ring with certain ring width $w_d$ (see Figure \ref{fig:hist}), approaching a steady state after $T=600\Omega_K^{-1}$.
Next, we gradually introduce dust feedback to the gas by increasing dust abundance gradually to the desired value with over a time interval of $30 \Omega_K^{-1}$.
We consider two {\it global} dust abundances $Z=0.01$ and $Z=0.02$,
and these runs are labeled by ``$Z$" followed by a number indicating the $Z$ values.
In all simulations, there are on average 2 particles per cell, ensuring sufficient particle statistics for studying dust dynamics.

We choose simulation box size to be $8 H\times 16 H \times H$ in $x$, $y$ and $z$ directions respectively.
We choose a fiducial resolution of 64/H (in the $x$ and $z$ directions; the resolution in the $y$ direction is half of that).
We further conducted simulations with a higher resolution of 128/H for comparison, labeled with an ``$H$" at the end of the run names.

\begin{figure*}[t]
\begin{center}
\includegraphics[width=.9\textwidth]{Figures/snapshot_xz_tsim1500.png}
\includegraphics[width=.9\textwidth]{Figures/snapshot_xy_tsim1500.png}
\end{center}
\caption{Snapshots of projected dust density in the $x$-$z$ (upper four panels) and $x$-$y$ (lower four panels) planes for runs Z0, Z2H, Z0w, Z2w, at time $T=1380\Omega_K^{-1}$, when the secondary filament in Z2w is clearly seen. Dust density in the $x$-$z$ plane is azimuthally averaged over the whole domain, while in the $x$-$y$ plane, dust density is vertically averaged around the midplane ($z = \pm 0.05H$) to highlight dust structures.  \label{fig:xyxz}}
\end{figure*}

\section{Results} \label{sec:results}

An overview of our simulations is shown in Figure \ref{fig:intro}, featuring a 3D visualization of the dust spatial distribution for our run Z2H at time $T=1380\Omega_K^{-1}$.
Our forcing procedure generates a wide and steady gas bump, with additional density fluctuations (i.e., density waves) as a result of the MRI turbulence \citep{heinemann09}. Clearly, the dust concentrates to the center of the gas bump, whose width is much narrower than the gap bump width. Besides settling and trapping, and the dust distribution is highly non-smooth and clumpy under background turbulence.

In this section, we present our main simulation results.
We measure the maximum dust density $\rho_{d, {\rm max}}$, normalized by background gas density $\rho_0$, as a proxy for dust clumping. Dust density in no-feedback runs (Z0 and Z0H) are rescaled to Z=0.02 for better comparison. For reference, we consider $\rho_{d, {\rm max}}\gtrsim200\rho_0$ as a threshold condition for triggering gravitational collapse, and hence planetesimal formation, which is applicable for typical outer disk conditions (see Equation 16 and Section 5.1 in \citetalias{xu22}). We further measure the dust ring width $w_d$ and dust scale height $H_d$, defined as the rms value of the radial and vertical coordinates for all dust particles. The time history of $\rho_{d, {\rm max}}/\rho_0$, $w_d$ and $H_d$ from all simulations is shown in Figure \ref{fig:hist}. The time-averaged values of $w_d$ and $H_d$ are presented in Table 1, averaged within $T=800-1500\Omega_K^{-1}$. These are accompanied by Figure \ref{fig:xyxz}, where we show the projected dust density in the $x-y$ and $x-z$ planes for Z0, Z2H, Z02 and Z2w runs at the saturated states.
We find that the results for the strong bump and weak bump cases are distinct, which we discuss separately below.

\subsection[]{The strong bump case: dust clumping in a single thin ring}\label{subsec:strong}

Top panel of Figure \ref{fig:hist} shows that dust density presents modest concentration in the absence of dust feedback, with $\rho_{d, {\rm max}}/\rho_0$ reaching $\sim100$, approaching yet barely reach the threshold for triggering planetesimal formation. This is to due to a combination dust trapping and intrinsic fluctuations. Reading from the second and third panels of Figure \ref{fig:hist}, we find a dust ring width of $w_d\sim0.55H$, and a scale height of $H_d\sim0.067H$ (comparable to $\sim 0.08H$ found in \citetalias{xu22} without enforced pressure bump). This yields a midplane density enhancement factor of $\sim50$ at bump center. Indeed, as shown in Figure \ref{fig:xyxz}, the mean midplane density at bump center for run Z0 (renormalized assuming $Z=0.02$) is about $\overline{\rho_d}\sim1$. We do not observe obvious dust clumping from the averaged dust densities along the vertical and azimuthal directions, and thus attribute an additional factor of $\sim50-100$ enhancement in $\rho_{d, {\rm max}}$ from localized density fluctuations. The results are largely independent of resolution, where $\rho_{d, {\rm max}}$, $w_d$ and $H_d$ are consistent with each other between run Z0 and Z0H within fluctuation levels.

Turning on dust feedback leads to a rapid boost in $\rho_{d, {\rm max}}/\rho_0$, well exceeding $10^3$ (and hence the Roche density) for all runs with Z=0.01 (Z1H) and Z=0.02 (Z2 and Z2H), which is therefore expected to trigger planetesimal formation. There are also visible bright spots seen in the bottom panel of Figure \ref{fig:xyxz} for run Z2H, suggesting that such clumping is physical instead of random Poisson fluctuations. There is a trend of higher resolution leading to higher level of dust clumping, where for runs Z1H and Z2H, $\rho_{d, {\rm max}}/\rho_0\gtrsim5000$, a factor of $\sim50$ enhancement compared to the no feedback case. We have further examined the cumulative probability density distribution (CDF, see Appendix \ref{appsec:cdf}), which reveal that dust density is indeed substantially boosted with a significant fraction (a few tens of percent) residing in over-Roche-density regions.

Somewhat surprisingly, with dust feedback, we observe a dramatic reduction of the ring width to $w_d\approx0.2H$,
$\sim 3$ times narrower than that without feedback. Moreover, the thickness of the dust layer is also reduced substantially by a factor of $3-4$. Examining Figure \ref{fig:xyxz} reveal that this reduction is not homogeneous: there is a strong/spiky dust concentration and more substantial settling at the bump center (see also the top left panel of Figure \ref{fig:vy_x1} and Figure 3 in \citetalias{xu22}), which reduces the measured $w_d$ and $H_d$ towards lower values. We note that the reduction of $H_d$ was also reported in \citetalias{xu22}, which was mainly attributed to reduction of turbulence correlation time.
However, in the presence of pressure bump, $H_d$ is further reduced by a factor of nearly $\sim 2$.
The dramatic reduction of ring width, as well as further promotion of dust settling, are likely the outcome of strong dust mass loading in the bump center (see further discussions in \S \ref{subsec:discu.interplay} and Appendix \ref{appsec:diff}).
The results are convergent with resolution, as the values of $w_d$ and $H_d$ are consistent with each other between runs Z2 and Z2H.

\subsection[]{The weak bump case: diverse ring properties}

The most prominent feature for the weak bump case compared to the strong bump case is the dust ring properties.
The rest of the results, especially dust clumping and settling, are similar to the case with a strong bump.

The middple panel of Figure \ref{fig:hist} shows that without dust feedback, $w_d$ in the weak bump case is generally higher compared to the strong bump case by a factor of up to $\sim2$. This is easily accounted for due to the reduced gas pressure gradient in the two sides of the gas bump, resulting in weaker dust trapping. This is also seen in Figure \ref{fig:xyxz}, showing a smooth and radially extended distribution of dust from run Z0w, together with azimuthally extended stripes akin to the MRI density waves.

When turning on dust feedback, interestingly, the dust ring width $w_d$, measured from the rms radial coordinates, remains similar for all cases.
However, when examining the dust spatial distribution in Figure \ref{fig:xyxz}, we find that
instead of concentrating into a single ring at the bump center, dust dynamics in weak bump is more complicated, potentially forming multiple fine-scale filaments within the gas bump. In the case of run Z2w, we observe two dust rings where the primary is close to bump center, and a secondary ring more than $H$ away from the bump center, in addition to a few minor filaments.
The formation of this more complex ring/filamentary structure is likely caused by the interplay between dust feedback and gas dynamics within the bump. In addition, these structures generally stay while slowly evolve over tens to hundreds of $\Omega_K^{-1}$. Details of the formation and evolution of (multiple) fine-scaled ring structures will be further discussed in \S \ref{subsec:discu.interplay} and Appendix \ref{appsec:longterm}.

\section{Discussion} \label{sec:discu}

\subsection{Clumping conditions}
\label{subsec:discu.condition}

For both the strong and weak bumps, sufficient dust clumping is seen to trigger planetesimal formation. Here we further discuss the clumping conditions to gain insight into the clumping mechanism. In doing so, we show in the top four panels of Figure \ref{fig:vy_x1} the radial profiles of vertically and azimuthally-averaged dust density, gas density and Keplerian-subtracted azimuthal gas velocity $v'_y$, at the very midplane of our simulation ($\pm$ 1 cell). The latter two are useful for quantifying the pressure gradient, and we further show the time evolution of the radial profiles of $\overline{v'_y}$ in four representative runs in the bottom panels.

\paragraph{Dust feedback}

Dust feedback is required for dust clumping in our simulations, as no clumping is seen in our run without feedback (Z0 and Z0w).
As seen in the upper panels of Figure \ref{fig:vy_x1}, the dust midplane density in the no-feedback cases generally exhibits a Gaussian-like radial profile, conforming to the standard expectation of concentration balancing turbulent diffusion \citep{dullemond18}. Including dust feedback, the radial density profiles substantially deviate from Gaussian, showing density cusps at the location of dust rings, where dust density is enhanced by a factor of near $100$ relative to the no-feedback counterparts. The cusps are associated with the locations of dust clumping, and can be the consequences of dust clumping itself.

\paragraph{Solid abundance}
Dust clumping occurs in all simulations for mean solid abundance $Z$ (averaged over the simulation domain) as low as $0.01$. Note that due to dust trapping, the effective solid abundance (i.e., vertically integrated dust-to-gas mass ratio) in the bump region is much higher.
We may estimate the enhancement factor to be $\sim 4-8$ in the weak bump and strong bump cases based on the value of $w_d$, corresponding to $Z\gtrsim0.04$ around bump center in all our simulations.
Therefore, the overall solid abundance $Z$ in the bulk disk appears less relevant, as dust trapping can easily bring sufficient amount of dust to the bump region.

\paragraph{Local pressure maxima}

In the strong bump case (runs Z1H, Z2H), dust clumping primarily occurs at the central cusp, which is marked in a dashed line in Figure \ref{fig:vy_x1}. We further see from the radial profiles of $\rho_g$ and $v'_y$ that this location exactly coincides with the location of the gas pressure maxima, corresponding to $v'_y=0$.
In the weak bump case, the primary ring close to the bump center is also located
at the pressure maxima, which in this case is widened into a plateau with a width of $\sim0.5H$ as seen from the gas density profile as well as an extended region with $v'_y\approx0$. The cusp where strong dust clumping occurs is located within the plateau, again consistent with the requirement of local pressure maxima, or at least a zero-pressure-gradient region being a necessary condition for strong dust clumping.
These findings are consistent with the results reported in \citetalias{xu22}, where the pressure bumps were not enforced but generated by the MRI zonal flows (albeit somewhat unrealistic in local unstratified simulations).

\paragraph{Weak clumping in the secondary ring}
For the weak bump case of Z1w and Z2w, in addition to the primary ring near the bump center, an off-centered ``secondary'' ring is identified at $x \sim -1.2H$.
We show in Appendix \ref{appsec:cdf} that the dust density CDF in the primary ring is similar to that in the strong bump case, with strong clumping. However, the dust CDF of the secondary ring is significantly different. The mean dust density is more than a factor of $\sim10$ lower, and only a very small fraction ($\lesssim3\%$) of particles reside in over-Roche-density regions (which exhibit as individual spots seen in the bottom right panel of Figure \ref{fig:xyxz}), indicating that clumping is much less efficient.
Examining the density (and hence pressure) profiles in the top right panel of Figure \ref{fig:vy_x1}, the location of the secondary ring is clearly not at the pressure maximum.

Overall, for both strong and weak bump cases, dust clumping is tightly related to the presence of pressure maxima as long as dust backreaction is included. While secondary rings/filaments can form for the weak bump case, they generally reside outside of the pressure maxima and do not show prominent/efficient clumping.

\begin{figure*}[t]
\begin{center}
\includegraphics[width=.9\textwidth]{Figures/pdens_vyx1_Hw.pdf}

\includegraphics[width=.9\textwidth]{Figures/vy1d_all.png}
\end{center}
\caption{Upper four panels: radial profiles of vertically and azimuthally averaged dust density ($\rho_d$, solid curves in upper two panels), gas density ($\rho_g$, dashed curves in upper two panels), and Keplerian-subtracted azimuthal gas velocity ($v'_y$, solid curved in lower two panels), at the midplane ($\pm 1$ cell from $z=0$), averaged between $T=800-1500\Omega_K^{-1}$. Vertical dashed lines mark the locations of the dust rings (peak of dust density), with colors corresponding to that used for radial profiles of each run.\label{fig:vy_x1} \\ Lower four panels: time evolution of the radial profile of $v'_y$, between $T=800-1500\Omega_K^{-1}$.}
\end{figure*}

\subsection{Interplay between dust feedback and gas bump}
\label{subsec:discu.interplay}

\paragraph{Bump/ring profile}

In addition to dust clumping, the interaction between dust and gas also influences gas bump properties, leading to diversity in dust structures, as we discuss here.

Dust feedback changes the radial profile of the gas bump.
This can be seen in the yellow region in the time-evolution map of $v'_y$ at the midplane in Figure \ref{fig:vy_x1} by comparing the upper and lower panels. In the strong bump case, as well as around the primary ring in the weak bump case, dust feedback leads to a plateau of Keplerian rotation at the bump center, in comparison to that in run Z0 and Z0w. This is consistent with the previous finding in \citep{huang20}.

The secondary dust ring, despite not in the pressure maxima, can be sustained for at least dozens of orbits without drifting towards the bump center.
Examining the profile of $v'_y$, we find that its value shows a dip at the secondary dust filament, implying near-Keplerian rotation. This is the result of feedback from the dusty filaments modifying the rotational profile of the background gas. It occurs in our weak bump runs likely because the gas pressure gradient is not strong enough to force the heavily loaded dust filaments to drift to bump center, resulting in dust trapping.\footnote{
One might argue that the formation of secondary dust rings seen in run Z1w and Z2w be a consequence of the special timing when we turn on dust feedback due to pre-existing local density enhancements from the MRI turbulence (e.g., zonal flows). We thus further perform a simulation where dust feedback is turned on at a much later time $T_{sim}=1200\Omega_K^{-1}$.
We still find similar secondary dust rings/filaments in this case.}
Such a location shares similarities with a pressure maxima in that there is almost no radial or azimuthal drift between gas and dust, and is likely responsible for the weak clumping discussed in \S \ref{subsec:discu.condition}.
These secondary rings also evolve over the timescale of tens to hundreds of orbits,
as we discuss further in Appendix \ref{appsec:longterm}.

\paragraph{Turbulent diffusion}
The enhanced dust mass loading in turbulent dust rings also affects the properties of background MRI turbulence.
With feedback included, we find that the gas rms velocities in both the radial and vertical directions are largely reduced in the dust ring (see details in Appendix \ref{appsec:diff}).
Such reductions can sufficiently explain the results reported in Section \ref{subsec:strong}, where the dust rings in the strong bump case becomes much narrower when dust feedback is included, and that dust becomes even more settled than that reported in \citetalias{xu22}.

\subsection{Comparison with other recent studies}

Recently, \cite{carrera21,carrera22a,carrera22b} conducted a series of hydrodynamic simulations of the SI, imposing a pressure bump by a similar (albeit more dissipative) forcing procedure as ours but without external turbulence. However, the pressure bump in most of their simulations are not sufficiently strong to overcome the background pressure gradient. As a result, particles get accumulated (but {\it not trapped}) in regions with lower pressure gradient, making such regions the most favorable for clumping and planetesimal formation \citep[e.g.,][]{bai10c}. In other words, their simulations identified the favorable conditions for dust clumping by the conventional SI in the presence of radial-dependent pressure gradient. Our simulations, on the other hand, focus on planetesimal formation in {\it dust-trapping rings}, i.e., the pressure maxima. The problem is only well-defined in the presence of external turbulence, as expected from observations \citep{dullemond18,rosotti20}, otherwise, dust would indefinitely get concentrated towards ring center over the trapping timescale (this issue may also apply to the analytic work of \cite{auffinger18}), forming planetesimals along the way. Overall, the physical picture in our simulations is distinctly different from the works of \cite{carrera21,carrera22a}, and likely better reflect outer disk conditions. Together with \citetalias{xu22}, we showed that for the MRI turbulence, the condition for planetesimal formation requires the presence of pressure maxima, again different from the conventional SI scenario.

The outer disk may also be subject to turbulence driven by the vertical shear instability \citep[VSI,][]{nelson13}. Using 2D global simulations, \cite{schafer20} found that the SI can co-exist with the VSI, and the VSI even facilitates the development of the SI and dust clumping, likely due to the overdensities induced by the VSI. We note that the overdensities from the VSI itself is also insufficient to overcome the background pressure gradient.
This is different from the case with the MRI, where we found that dust clumping requires the presence of pressure maxima, suggesting that dust clumping conditions likely depends on the nature turbulence.
\cite{lehmann22} found that the presence of an initial pressure bump (but without forcing) further facilitates dust concentration, though their resolution is insufficient to capture dust clumping. Moreover, the VSI leads to vortex formation in 3D \citep{manger20}, and hence dust clumping in dust-laden vortices \citep[e.g.,][]{raettig21}, which represents another alternative for planetesimal formation in turbulent disks. It is worth noting that the VSI and the MRI can co-exist (Cui \& Bai 2022, submitted), and turbulence is dominated by the MRI/VSI when the ambipolar Elsasser number $Am$ is large/small compared to unity. It is likely that clumping conditions still depend on the nature of turbulence (be it due to the VSI or the MRI), although more investigations are needed to understand the clumping mechanisms, and how the dust dynamics is affected when both instabilities are present.

\subsection{Observational implications}

\paragraph{Dust clumping and fine-scale substructures}

There can be multiple observational diagnostics to probe the turbulent dust rings. First, as such dust rings are clumpy, there can be both optically thick and optically thin regions within the ring. Observations at sub-mm by ALMA can hardly resolve such regions, but the clumpiness could be reflected in the spectral indices, as suggested by \cite{scardoni21}. Further analysis of our results with radiative transfer and comparison with ALMA observations may provide evidence of dust clumping in the observed dust rings.

Second, we have seen the formation of secondary rings/filaments in our weak-bump simulations. These fine-scale substructures are separated on scales $\lesssim H$, but are all located within a single gas pressure bump. We note that generally we do not expect ring-like substructures in the gas component on scales $\lesssim H$ \citep{dullemond18}, whereas observationally, evidence of such fine-scale ring-like substructures is mounting \cite{jennings22}, with hints of multiple rings within previously identified rings. Our results provide a natural explanation of forming multiple sub-$H$ scaled substructures within a single gas bump. For instance, they may provide an alternative explanation to the compact double rings observed in HD 169142 at 57 and 64 AU separately, previously interpreted as the consequence of a migrating protoplanet \citep{perez19}.

Third, we note that the dust rings from our simulations are not entirely axisymmetric, with fine-scale substructures along the azimuth. Although our simulations are local and do not cover the entire azimuth, we notice that fine-scale, low-contrast asymmetries have been observed in some systems \citep[e.g., DM Tau,][]{hashimoto21} and may be common among gapped protoplanetary disks \citep{vandermarel21}.

\paragraph{Disk turbulence}

Dust ring width has recently been used to constrain the level of turbulence (the $\alpha$ parameter) in disks \citep{dullemond18,rosotti20}.
However, as discussed in \S \ref{subsec:discu.interplay}, disk turbulence in the presence of a dust ring can be both anisotropic and non-uniform due to dust feedback.

More specifically, turbulent diffusion in the radial direction in our simulations is overall stronger than that in the vertical direction, leading to a large aspect ratio of $w_d/H_d \sim 10$ for the single dust ring in run Z2H, which is 5 times of aspect ratio of the gas bump inserted in the simulation. This effect is even more extreme for the weak bump case, with $w_d/H_d \sim 30$ for Z2w, due to the formation of multiple dust filaments.
This flat morphology of dust ring is consistent with recently reported highly settled disk with relatively extended rings in Oph 163131 \citep{villenave22}, likely suggesting different turbulent levels in radial and vertical directions.

In addition, dust feedback weakens turbulent diffusion at the dust ring center (see \S \ref{subsec:discu.interplay}).
It makes the dust ring narrower in our strong bump case, but the presence of additional off-centered filaments increases the overall $w_d$ (if they are unresolved), which largely cancels out the effect of reduced ring width.
Thus, there is a degeneracy between turbulence level, dust feedback, and bump amplitude in determining the dust ring width.
A single $\alpha$ parameter is not necessarily a good indicator of disk turbulence level, and special care is needed when interpreting the level of disk turbulence from observations of dust rings.

\paragraph{Prospects for future observations}
Overall, our simulations suggest that the dynamics of dust-trapping rings can be constrained by the presence or absence of fine-scaled substructures. Simultaneously constraining the gas bump and characterizing the fine-scale dust substructures are crucial for testify such expectations.
High resolution observations searching for kinematic features (e.g. ExoALMA \footnote{\url{https://www.exoalma.com}}) will provide essential information on the properties of the gas pressure bump, while identifying sub-$H$ scaled substructures requires observations with the longest baseline by ALMA for nearby disks, and also potentially by future upgrade of ALMA \citep[e.g.,][]{burrill22}, ngVLA \citep[e.g.,][]{andrews18b} and/or SKA \citep{ilee20}.

\section{Conclusions and Perspectives} \label{sec:conclusions}

In this paper, we conduct the first controlled experiments of dust dynamics and clumping in the MRI-turbulent pressure bumps. Our simulations include realistic level of the MRI turbulence applicable to the ambipolar-diffusion dominated outer disks, and a carefully-designed forcing scheme to achieve self-consistent balance between dust trapping and turbulent diffusion. Our main findings include:

\begin{itemize}

\item At solid abundance compatible to solar, we find robust and efficient dust clumping in radial pressure maxima, with a sizable fraction of dust mass residing in clumps over the Roche density. Secondary dust filaments may form outside of the pressure maxima, but shows only weak or no clumping.

\item Dust feedback substantially affects ring properties. Strong mass loading in dust rings reduces turbulent diffusion, making the ring narrower with stronger dust settling. Weaker pressure bump can result in formation of meta-stable secondary dust filaments near the bump.

\end{itemize}

Our work strongly support that dust-trapping rings are robust sites for planetesimal formation. Our simulations do not include self-gravity, and it is yet to examine how efficient the dust in dense clumps get converted into planetesimals. This could result from a competition between self-gravity and turbulent diffusion \citep{gerbig20}, though the internal turbulent diffusivity is uncertain (but see \cite{klahr20,klahr21}) and requires further study.

Our results also have important observational implications. In particular, depending on the strength of the gas bump, as well as dust abundance, dust feedback leads to a diversity of dust ring morphologies that may evolve over time. Conventional estimate of the level of turbulence based on ring width and dust thickness does not necessarily reflect true level of disk turbulence. Moreover, our results suggest the presence of fine-scale substructures that may be resolvable at the highest resolution by ALMA and the future ngVLA.

As a first study, we consider only a single dust species with fixed stopping time $\tau_s=0.1$. The natural next step is to incorporate a dust size distribution. This will allow us to further examine how particles of different sizes participate in dust clumping, and to yield more realistic observational signatures. Moreover, the pressure bumps in our shearing-box simulations rely on an artificial forcing prescription. Future work should consider more realistic bump formation scenarios such as magnetic flux concentration and planet-induced gas bumps in the global context.

\section*{Acknowledgements}

We thank Daniel Carrera, Kees Dullemond, Guillaume Laibe and Chao-Chin Yang for useful discussions, and Greg Herczeg for overseeing the completion of this work. This work is supported by the National Key R\&D Program of China No. 2019YFA0405100, and the China Manned Space Project with NO. CMS-CSST-2021-B09. ZX Acknowledges the support of the ERC CoG project PODCAST No. 864965 and NSFC project 11773002. Numerical simulations are conducted on the Orion cluster at Department of Astronomy, Tsinghua University, and on TianHe-1 (A) at National Supercomputer Center in Tianjin, China.

\bibliographystyle{apj}


\section{Introduction}
\label{sec:intro}

Scintillation of global navigation satellite system (GNSS) signals is caused by electron density irregularities that are present within the ionosphere \cite{crane1977ionospheric}. Ionospheric scintillations induce significant fluctuations in global navigation satellite systems (GNSS) signals. These fluctuations include variations in the signal-to-noise ratio (SNR), a phenomenon known as amplitude scintillation. Amplitude scintillation may lead to degradation of the quality of service provided by GNSS satellites to users. A typical measure of amplitude scintillation is the S4 index, which is defined as the ratio of the standard deviation of the signal power to the average signal power \cite{yeh2019superposition}. Scintillation monitoring receivers typically provide the S4 index as they continuously monitor the SNR of GNSS signals.\par
As ionospheric scintillations have varying effects based on the geographical location of the user, knowledge of the spatial distribution of scintillation-causing ionospheric irregularities would be very useful. Furthermore, the ability to forecast the severity of these scintillations based on historical data would be advantageous when real-time data is unavailable. Motivated by the points mentioned earlier, in this work, we study the possibility of using historical data from a single GPS scintillation receiver to train a machine learning (ML) model to forecast the severity of amplitude scintillation with respect to temporal and spatial parameters.\par

\subsection{Related Works and Contribution}

Several works have utilized ML to predict or forecast amplitude scintillation, such as \cite{rezende2010survey,linty2018detection,savas2019impact,dey2021automatic,atabati2021ionospheric}. Out of these, the most recent and related to our work are \cite{dey2021automatic,atabati2021ionospheric}. In \cite{dey2021automatic}, the issue of classifying ionospheric amplitude scintillation was addressed by attempting to solve a binary classification problem by considering two levels of scintillation (scintillation or no scintillation). Their best performing classifier was the XGBoost method, which achieved an accuracy of $99.76\%$ for testing. However, they did not consider the spatial allocation of amplitude scintillation even though the variations of amplitude scintillation depend heavily upon spatial parameters. Thus, by not mapping the scintillation values, it would be difficult to localize the predicted scintillation. Furthermore, different scintillation intensities have varying effects on ground users, and thus, unlike \cite{dey2021automatic}, we consider the problem as a 3-class classification problem, by considering weak, moderate, and severe levels of scintillation.\par
In \cite{atabati2021ionospheric}, the use of neural networks for the prediction of S4 values was proposed. The problem was addressed as a regression analysis problem. Furthermore, they included different parameters related to space weather such as F10.7, sunspot number, and peak electron density height in the prediction process. Their designed model was able to predict daily ionospheric scintillations with $81\%$ accuracy. However, like \cite{dey2021automatic} they did not consider the effect of spatial parameters on the predicted scintillation values. Furthermore, their utilization of F2 peak electron density height and F2 vertical drift velocity necessitates the use of ionosonde data, which is less common compared to scintillation monitoring receivers.\par
This work attempts to solve the problem of forecasting the severity of spatially and temporally allocated amplitude scintillation values using bagged trees. This is done by solving a 3-class classification problem using ML methods, to predict weak, moderate, and severe scintillation. The S4 data used in this work was retrieved from a scintillation monitoring receiver at the McMurdo Station in Antarctica, corresponding to a period of four years.\par

\section{Data and Methodology}

In this work, we use 1-minute observations of S4 from a GPS scintillation monitoring receiver at McMurdo station ($-77.83\degree$, $166.66\degree$) in Antarctica, from 18-Jan-2011 to 14-Nov-2014. These observations are freely available at \url{http://cedar.openmadrigal.org/}. To map the S4 values to a geographical location, the thin-shell ionospheric model was assumed at an altitude of 350 km. The point at which the GNSS signal intersects with the thin shell is known as the ionospheric pierce point (IPP). The IPP is mapped into latitude and longitude coordinates. The observation data was provided as Ionospheric Scintillation Monitoring Records (ISMR) files. We extracted the time, satellite elevation, latitude and longitude coordinates of the IPP, and S4 values from these files.\par
One of the most important steps in ML is the preparation and preprocessing of data to improve the learning ability of ML models. For preprocessing, several steps were followed, as shown below:
\begin{enumerate}
  \item Set an elevation cutoff of 20 degrees, to eliminate the effects of multi-path, and removed any $S4<0$ values.
  \item Removed discontinuity in longitude coordinates by adding $360\degree$ to negative longitude values.
  \item At this point the total number of observations was $12\,377\,607$, $67.91\%$ of which were below $0.05$. We assumed S4 values below $0.05$ to represent no scintillation and they were therefore discarded. This reduced the total number of observations to $3\,972\,227$. This served two main purposes. First, removing this large number of S4 observations under $0.05$ would prevent the models from being biased toward very low S4 values. Second, the number of observations would be more manageable, and take considerably less time to train and evaluate.
  \item Added three commonly used solar indices to represent solar activity. These are the daily averaged Planetary Kp-index (KP), Sunspot Number (SSN), and F10.7 index. These indices were obtained from \url{https://omniweb.gsfc.nasa.gov/form/dx1.html}. Next, observations with missing indices, i.e., F10.7 $=999$, were discarded. Furthermore, we did not include parameters such as ionospheric vertical drift velocity and peak electron density height \cite{rezende2010survey,atabati2021ionospheric} to eliminate dependence on ionosonde data.
  \item The data was then divided into three classes. Class 1 represents $S4<0.2$, class 2 represents $0.2\leq S4<0.3$, and class 3 represents $S4\geq 0.3$. These classes correspond to weak, moderate, and severe amplitude scintillation respectively.
  \item For classes 1, 2, and 3, the number of observations were $3\,788\,881$, $157\,163$, and $22\,959$, respectively. Thus, the total imbalanced dataset had a total of $3\,969\,003$ observations.
  \item We also created a balanced dataset by taking $22\,959$ random samples, corresponding to the smallest class size (class 3), from each class of the imbalanced dataset, adding up to a total of $68\,877$ observations. The balanced dataset should provide another perspective of the performance of ML models as imbalanced data was found to substantially compromise the learning process of ML models \cite{he2009learning}.
\end{enumerate}

After preprocessing, the variables in the dataset consisted of: day of year (DOY), hour of day (HOD), IPP latitude coordinate, IPP longitude coordinate, KP, SSN, F10.7, and S4 class (either 1, 2, or 3). The ML models considered in this work are discussed in the following paragraphs.\par\noindent
\underline{Decision Trees (DT)} consist of nodes, branches, and leaves. It is constructed by splitting data iteratively into branches to form a hierarchical structure. The upper edge of the tree starts from the root node, and the lower edge of the tree consists of leaf nodes. By following a path from the root node to the leaf nodes, the user can obtain a classification from the DT model \cite{sutton2005classification}.\par\noindent
\underline{Naive Bayes (NB)} classifiers assume the attribute values are conditionally independent of each other given the target value or class \cite{chen2020novel}. Thus, in our case, each feature contributes to predicting an S4 value without depending on each other.\par\noindent
\underline{Support Vector Machines (SVM)} create the best decision boundary that can segregate data into their respective classes. This decision boundary is known as the hyper-plane. The data which are closest to the hyper-plane and affect its position, are termed as support vectors, hence the name \cite{steinwart2008support}.\par\noindent
\underline{K-Nearest Neighbors (KNN)} are algorithms that assume that similar data points neighbor each other. Accordingly, it classifies the data to the nearest corresponding class \cite{liao2002use}.\par\noindent
\underline{Boosted Trees} are a type of supervised learning where a model combines the weak learners in an iterative manner and progressively improves the learners. Such method uses a weighted average of the results obtained from applying a prediction model to various samples \cite{sutton2005classification}.\par\noindent
\underline{Bagged Trees} is a method where multiple versions of a classifier, in this case a DT, are generated with each of them generating their own predictions, and the final prediction is the average of all of these individual predictions \cite{sutton2005classification}.

\begin{table*}
\centering
\caption{Results of the unoptimized ML models using the balanced and imbalanced datasets.}
\label{tab1}
\begin{tabular}{|l|l|c|c|}
\hline
Model & Parameters & \multicolumn{1}{l|}{Balanced Dataset Accuracy} & \multicolumn{1}{l|}{Imbalanced Dataset Accuracy} \\ \hline
Decision Tree & \begin{tabular}[c]{@{}l@{}}Maximum number of splits = 100\\ Split criterion = Gini's Diversity Index\end{tabular} & $67.8\%$ & $95.5\%$ \\ \hline
Naive Bayes & Kernel type = Gaussian & $56.1\%$ & $95.5\%$ \\ \hline
SVM & \begin{tabular}[c]{@{}l@{}}Kernel function = Gaussian\\ Kernel scale = 0.66\end{tabular} & $64.1\%$ & $95.5\%$ \\ \hline
KNN & \begin{tabular}[c]{@{}l@{}}Preset = Weighted KNN\\ Number of neighbors = 10\\ Distance metric = Euclidean\\ Distance weight = Squared Inverse\end{tabular} & $60.4\%$ & $94.6\%$ \\ \hline
Boosted Trees & \begin{tabular}[c]{@{}l@{}}Method = AdaBoost\\ Maximum number of splits = 20\\ Number of learners = 30\\ Learning rate = 0.1\end{tabular} & $62.1\%$ & $95.5\%$ \\ \hline
Bagged Trees & \begin{tabular}[c]{@{}l@{}}Maximum number of splits = $68\,880$\\ Number of learners = 30\end{tabular} & \boldmath{$80.1\%$} & \boldmath{$97.1\%$} \\ \hline
\end{tabular}
\end{table*}

\section{Results and Discussion}

For training and validation, the balanced dataset was randomly split using 10-fold cross-validation \cite{refaeilzadeh2009cross}, while the imbalanced dataset was randomly split using the holdout technique (training: $90\%$, validation: $10\%$).\par
In this work, the performance of commonly used ML methods was tested in terms of their accuracy, precision, and recall \cite{davis2006relationship} (see Table \ref{tab1}). Accuracy in this case is defined as the number of correct predictions divided by the total number of predictions, and can be represented by the following:
\begin{equation}
\label{eq1}
\text{Accuracy}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}
\end{equation}
where TP, TN, FP, and FN are the number of true positives, true negatives, false positives, and false negatives, respectively. Furthermore, precision and recall can be defined as:
\begin{equation}
\begin{split}
\label{eq2}
\text{Precision}&=\frac{\text{TP}}{\text{TP}+\text{FP}}\\
\text{Recall}&=\frac{\text{TP}}{\text{TP}+\text{FN}}.
\end{split}
\end{equation}

As shown in Table \ref{tab1}, the best performing model was the bagged trees for both balanced and imbalanced datasets. It achieved $80.1\%$ accuracy using a balanced dataset and $97.1\%$ accuracy using the imbalanced dataset. We note that the baseline accuracy of the balanced dataset is $33.3\%$ since all categories of the balanced dataset have an equal number of observations. The baseline accuracy can be achieved by a simple model that randomly allocates S4 values to either one of the three categories. Furthermore, the baseline accuracy of the imbalanced dataset is $95.46\%$, corresponding to the percentage of observations of class 1 in the imbalanced dataset, and is achieved by assuming that all predictions belong to class 1. Therefore, it is clear that the accuracy improvement gained by the model utilizing the balanced dataset is considerably higher than the baseline, by around $46.8\%$.\par

To improve the performance of the bagged trees model, Bayesian optimization was used to tune its hyperparameters, namely, the maximum number of splits and the number of learners. The optimal hyperparameters for each dataset are presented in Table \ref{tab2}.\par

\begin{table}
\centering
\caption{Optimal hyperparameters after 50 iterations of Bayesian optimization.}
\label{tab2}
\begin{tabular}{|l|c|c|}
\hline
Optimized Variables & \multicolumn{1}{l|}{\begin{tabular}[c]{@{}l@{}}Balanced\\ Dataset\end{tabular}} & \multicolumn{1}{l|}{\begin{tabular}[c]{@{}l@{}}Imbalanced\\ Dataset\end{tabular}} \\ \hline
\textbf{Maximum Number of Splits} & $378$ & $485$ \\ \hline
\textbf{Number of Learners} & $68\,828$ & $537\,360$ \\ \hline
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Confusion matrix for balanced dataset, using optimized Bagged Trees.}
\label{tab3}
\begin{tabular}{|cc|ccc|c|}
\hline
\multicolumn{2}{|c|}{\multirow{2}{*}{}} & \multicolumn{3}{c|}{\textbf{Ground Truth}} &  \\ \cline{3-6}
\multicolumn{2}{|c|}{} & \multicolumn{1}{c|}{Class 1} & \multicolumn{1}{c|}{Class 2} & Class 3 & \emph{Precision} \\ \hline
\multicolumn{1}{|c|}{\multirow{3}{*}{\rotatebox[origin=c]{90}{\textbf{Predict.}}}} & Class 1 & \multicolumn{1}{c|}{$19\,225$} & \multicolumn{1}{c|}{$2\,076$} & $276$ & $89.11\%$ \\ \cline{2-6}
\multicolumn{1}{|c|}{} & Class 2 & \multicolumn{1}{c|}{$3\,336$} & \multicolumn{1}{c|}{$17\,063$} & $2\,780$ & $73.61\%$ \\ \cline{2-6}
\multicolumn{1}{|c|}{} & Class 3 & \multicolumn{1}{c|}{$398$} & \multicolumn{1}{c|}{$3\,820$} & $19\,903$ & $82.51\%$ \\ \hline
\multicolumn{1}{|c|}{} & \emph{Recall} & \multicolumn{1}{c|}{$83.74\%$} & \multicolumn{1}{c|}{$74.32\%$} & $86.69\%$ & \begin{tabular}[c]{@{}c@{}}\emph{Accuracy}\\ $81.58\%$\end{tabular} \\ \hline
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Confusion matrix for the imbalanced dataset, using optimized Bagged Trees.}
\label{tab4}
\begin{tabular}{|cc|ccc|c|}
\hline
\multicolumn{2}{|c|}{\multirow{2}{*}{}} & \multicolumn{3}{c|}{\textbf{Ground Truth}} &  \\ \cline{3-6}
\multicolumn{2}{|c|}{} & \multicolumn{1}{c|}{Class 1} & \multicolumn{1}{c|}{Class 2} & Class 3 & \emph{Precision} \\ \hline
\multicolumn{1}{|c|}{\multirow{3}{*}{\rotatebox[origin=c]{90}{\textbf{Predict.}}}} & Class 1 & \multicolumn{1}{c|}{$377\,197$} & \multicolumn{1}{c|}{$7\,969$} & $400$ & $97.83\%$ \\ \cline{2-6}
\multicolumn{1}{|c|}{} & Class 2 & \multicolumn{1}{c|}{$1\,627$} & \multicolumn{1}{c|}{$7\,460$} & $715$ & $76.11\%$ \\ \cline{2-6}
\multicolumn{1}{|c|}{} & Class 3 & \multicolumn{1}{c|}{$64$} & \multicolumn{1}{c|}{$288$} & $1\,180$ & $77.02\%$ \\ \hline
\multicolumn{1}{|c|}{} & \emph{Recall} & \multicolumn{1}{c|}{$99.56\%$} & \multicolumn{1}{c|}{$47.47\%$} & $51.42\%$ & \begin{tabular}[c]{@{}c@{}}\emph{Accuracy}\\ $97.21\%$\end{tabular} \\ \hline
\end{tabular}
\end{table}

In Tables \ref{tab3} and \ref{tab4} we present the confusion matrices of the optimized bagged tree model using the balanced and imbalanced datasets, hereinafter referred to as the balanced and imbalanced models, respectively. The accuracy improvement of the balanced model after optimizing its hyperparameters was higher, by around $1.4\%$, while the accuracy of the imbalanced model improved by about $0.1\%$.\par

The overall accuracy of the imbalanced model was higher by $15.63\%$, compared to the balanced model. However, this is because the imbalanced model is biased towards the class with the higher number of observations \cite{brodersen2010balanced}. This is evident because the class 1 recall and precision of the imbalanced model were higher by $15.82\%$ and $8.72\%$, respectively, compared to the balanced model. The downside of this is the comparatively poor class 3 classification accuracy of the imbalanced model, as the class 3 recall and precision are $35.27\%$ and $5.49\%$ lower, respectively, compared to the balanced model. For class 2, the precision was $2.5\%$ higher for the imbalanced dataset, while the recall was $26.85\%$ higher for the balanced dataset. Therefore, the two models present the user with two different choices. If the user prioritizes the correct classification accuracy of weak scintillation (class 1) then the imbalanced model would be optimal. Otherwise, the balanced model would be preferable.\par

\section{Conclusion}
In this work, we attempted to solve the problem of forecasting the severity of spatially and temporally allocated amplitude scintillation values using machine learning. Six different models were tested, with the best performing of them being the bagged trees model. Next, Bayesian optimization was used to tune the hyperparameters and improve the classification accuracy. The bagged tree model trained using a balanced dataset achieved good accuracy for all three classes. On the other hand, the bagged tree model trained using an imbalanced dataset achieved excellent accuracy that was biased towards the classification of weak scintillation. We aim to expand the study by testing the models with data from different stations in different regions. Furthermore, the usage of deep learning models will be explored in future work as they may benefit from the large quantity of training data.\par

\bibliographystyle{IEEEbib}


\section*{Supplementary Material}

\subsection{A. Renormalization group analysis}
We begin by reiterating that our RG procedure is carried out for the Hamiltonian defined in Eq.~\ref{eq:ham}.
The effective actions
$S$ and $S_{\chi}$ are introduced to describe the next-order fluctuations and calculate the correlations functions.

The relevant correlation functions are:
\ba
&&\langle \sum_\mu [\phi^\mu(\bm{r},\tau=0)]^2 \rangle _{imp}  \, ,
\nonumber \\
&&\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle _{imp} \, .
\label{eq:correlation}
\ea
Here, $\langle \cdot \rangle _{imp}$ means we subtract the free-bulk contributions, and the bulk system has spacetime dimension $d=4-\epsilon$. Before performing calculations, we note that the exact formula of above correlations functions are
\ba
&&\langle  \sum_\mu \phi^\mu(\bm{r},\tau=0)^2\rangle_{imp}=\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }
\langle \sum_\mu S^\mu(\tau) S^\mu(\tau') \rangle_{imp} \nonumber \\
&& -\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle_{imp} = K
\bigg\langle
\bigg[
\delta_{\tau,\tau'}\delta_{k,k'}(\partial_{\tau'}  + \epsilon_{\bm{k}})
+ \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\sigma^\mu\bigg)
\bigg]
\bigg\rangle_{imp} -G_0(\bm{r}_0,\tau),
\label{eq:corr_exact}
\ea
where
\ba
h =-4A_\phi g^2\Gamma(-1+\epsilon)\sin(\frac{\epsilon \pi}{2}) \approx  [2\pi A_\phi +O(\epsilon) ]g^2.
\ea
$A_\phi$ is the normalization factor of boson propagator
and $G_0(\bm{r}_0,\tau)$ is the free fermion propagator. Above equations suggest the fermion and boson correlation functions can be represented as functional of $n$-point correlation functions of $S^\mu$ fields. Then, it's sufficient to evaluate the correlation functions of $S^\mu$. To do so, we first derive the effective action of $S^\mu$ by integrating out fermions and bosons:
\ba
S_{eff}
=-iSS_B
-K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg).
\ea
Using $S_{eff}$ and Eq.~\ref{eq:corr_exact}, we're able to calculate fermion and boson correlation functions. The detailed calculations are shown in Supplementary Material (SM) B and the results are
\ba
\langle \sum_\mu [ \phi^\mu(\bm{r},\tau=0)]^2\rangle_{imp}
&\propto &  \frac{1}{|\bm{r}|^2}
h\bigg[
1-
\frac{1}{S\pi r}\int_0^J\frac{F(j)}{j}dj +\frac{1}{S\pi \epsilon }[\frac{\log(1+h/c(J)}{c(J)} + \frac{\log(1+h/c(J)^*)}{c(J)^*}]
\bigg] \nonumber \\
\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle_{imp} &\propto&
\frac{J^2\text{sgn}(\tau)}{(1+J^2)|\tau|}
\bigg[
1
+ \frac{-1-J^2}{(1+J^2)^2\pi Sr}\int_0^J\frac{F(j)}{j}dj
+ \frac{1}{(1+J^2)\pi Sr}\int_0^J\frac{2jF(j)}{1+j^2}dj \nonumber\\
&&
-\frac{2}{(1+J^2)r\pi S}
\int_0^J \frac{H(j)}{1+j^2}dj -\frac{-1+J^2}{(1+J^2)^2 \pi S\epsilon}
\bigg(
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) \nonumber \\
&&
+\frac{2iJ}{(1+J^2)^2\pi \epsilon S}
\bigg( -
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg)
\bigg] . \nonumber \\
 \label{eq:corr}
\ea

We are now in the position to carry out renormalization group study. We introduce the following relation between renormalized couplings $J,h$ and bare couplings $J_B,h_B$:
\begin{eqnarray}
&&h_B = Z_h\mu^{\epsilon}h \nonumber \\
&&J_B = Z_J\mu^{r}J
\label{eq:bare_coupling_supp}
\end{eqnarray}
where $Z_J,Z_h$ are coupling renormalization factor, and $\mu$ is a renormalization energy scale.
Without loss of generality, we let
\ba
Z_h = 1 + \frac{1}{S}x_h \nonumber \\
Z_J = 1 + \frac{1}{S}x_J \nonumber.
\ea
We determine $x_J$ and $x_h$ by requiring all the poles as a function of $r,\epsilon$  in Eq.~\ref{eq:corr} are minimally removed. This leads to \ba
x_h &=& \frac{1}{\pi r}\int_0^J\frac{F(j)}{j}dj
-\frac{1}{\pi \epsilon }
\bigg[
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg] +O(\frac{1}{S})\nonumber\\
x_J&=& \frac{1}{2\pi r}\int_0^J\frac{F(j)}{j}dj
-\frac{1}{\pi r}\int_0^J\frac{jF(j)}{1+j^2} dj
+\frac{1}{r\pi}\int_0^J\frac{H(j)}{1+j^2}dj +\frac{-1+J^2}{2(1+J^2)\pi \epsilon}\bigg(
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) \nonumber \\
&&
-\frac{iJ}{(1+J^2)\pi \epsilon}
\bigg(
-\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) +O(\frac{1}{S}) \nonumber \\
\ea

Finally, we take $\mu$ derivative of Eq.~\ref{eq:bare_coupling_supp} at fixed $h_B$ and $J_B$ and obtain the following beta functions:
\ba
\beta_J &=&  \frac{J-J^3}{2S\pi (1+J^2)}F(J) + \frac{J^2}{S\pi(1+J^2)}H(J) +\bigg[ \frac{-(1-J^2)Jh+2iJ^2h}{2S(1+J^2)\pi c(J)(c(J)+h)}+\text{h.c.}
\bigg] +o(\frac{1}{S}) \nonumber  \, , \\
\beta_h& =& -\frac{1}{S}\tilde{\epsilon} h + \frac{h}{S\pi}F(J)
-\frac{1}{ S\pi }\frac{(c(J)+c(J)^*)h}{|c(J)|^2}
-\frac{1}{S\pi }\frac{-2h^2-(c(J)+c(J)^*)h}{(c(J)+h)(c(J)^*+h)} +o(\frac{1}{S})  \, .
\ea
where we also set $r=0$ at the last step.

Note that an alternative way to perform renormalization group study is to start from the action in Eq.~\ref{eq:action} and then expand in powers of $\chi$ fields~\cite{nahum2022fixed}. Combining the quadratic order term $S_{\chi}$ and fourth-order term $S_4$, we could find an interacting field theory of $\chi$: $S_{\chi}+S_4$. Then we are able to perform RG analysis based on this model. The resulting beta functions are expected to be consistent with current approach.
{We reserve such an analysis for a future study.}

\subsection{B. Calculations of correlation functions}
Here we present the detailed calculation of Eq.~\ref{eq:corr_exact} using the effective action $S_{eff}
$.
We first introduce the spinor representation of $S^\mu$ with $z^\dag \frac{\sigma^\mu}{2} z$ and let $z=z_0+\delta z$ (see Eq.~3). This gives us
\ba
S^\mu &= & S\delta_{\mu,z} + \delta S^\mu \nonumber \\
\delta S^x &=& \frac{i\sqrt{S}}{\sqrt{2}}(\chi^\dag -\chi) + O(\frac{1}{S^{1/2}})\nonumber \\
\delta S^y &=& \frac{\sqrt{S}}{\sqrt{2}}(\chi^\dag +\chi) + O(\frac{1}{S^{1/2}}) \nonumber  \\
\delta S^z &=& -|\chi|^2+O(\frac{1}{S})
.
\label{eq:s_exp}
\ea
We expand the action in Eq.~\ref{eq:action} to the quadratic order of $\chi$ fields:
\ba
S_{eff}
&=& S_0 +S_\chi \nonumber \\
S_0&=&-\frac{Sg^2}{4}\int_{\tau,\tau'} \frac{\sum_\mu(z_0^\dag \sigma^\mu z_0)^2}{2}G_\phi(\tau-\tau')
-K\text{Tr}\log\bigg(
\delta_{\tau,\tau'}\delta_{k,k'}(\partial_{\tau'}  + \epsilon_{\bm{k}})
+ \delta_{\tau,\tau'}J_K\sum_\mu (z_0^\dag \sigma^\mu z_0)\sigma^\mu\bigg) ,
\label{eq:simp_exp}
\ea
with $S_0$ the saddle point contribution at order of $S$ and $S_\chi$ describing fluctuations at order $S^0$. We now derive the formula of $S_\chi$. It has contribution from three part: the Berry phase, bosonic coupling and Kondo
interaction.
The Berry phase term
gives
\ba
\int i\Omega|\chi(\Omega)|^2 \frac{d\Omega}{2\pi} \, ,
\label{eq:berry_chi}
\ea
{while}
 the bosonic coupling
yields
\ba
&&-\frac{g^2}{2}\int_{\tau,\tau'}G_\phi(\tau-\tau')[\chi^\dag(\tau)\chi(\tau') + \chi(\tau)\chi^\dag(\tau') -|\chi(\tau)|^2 -|\chi(\tau')|^2] \nonumber \\
&=&\int h|\Omega|^{1-\epsilon} |\chi(i\Omega)|^2 \frac{d\Omega}{2\pi} \, ,
\label{eq:bose_chi}
\ea
and, finally, the
 Kondo coupling gives (by expanding the fermion determinant)
\ba
&&\frac{\kappa J_K^2}{4} \int_{\tau,\tau'}\bigg[ G_s(\tau-\tau')G_s(\tau'-\tau) -G_a(\tau-\tau')G_a(\tau'-\tau)\bigg]
\bigg[
\chi^\dag(\tau)\chi(\tau') + \chi^\dag(\tau')\chi(\tau) - |\chi(\tau)|^2 - |\chi(\tau')|^2
\bigg]\nonumber \\
&&+\frac{\kappa J_K^2}{2}\int_{\tau,\tau'}\bigg[
G_s(\tau-\tau')G_a(\tau'-\tau)
\bigg]
\bigg[
-\chi^\dag(\tau)\chi(\tau') + \chi(\tau)\chi^\dag(\tau')
\bigg] \nonumber \\
&=&\int \frac{2\kappa J^2 |\Omega|^{1-2r} }{\pi(1+J^2 |\Omega|^{-2r})^2}|\chi(i\Omega)|^2 \frac{d\Omega}{2\pi}+\int \frac{i\Omega \kappa }{\pi }
\bigg[
\frac{J|\Omega|^{-r} -J^3 |\Omega|^{-3r}}{(1+J^2|\Omega|^{-2r})^2} - \arctan(J|\Omega|^{-r})
\bigg] |\chi(i\Omega)|^2 \frac{d\Omega}{2\pi} ,
\label{eq:fermi_chi}
\ea
where $G_s(\tau) = \sum_\sigma G_\sigma(\tau)$ and $G_a(\tau) = \sum_\sigma \sigma G_\sigma(\tau)$ and $G_\sigma(\tau)$ denotes the conduction electron Green's function at leading order $(S^1)$. The explicit formula of these functions are
\ba
G_s(i\omega) &= \frac{-2i\rho_0\pi |\omega|^{-r}}{1+\rho_0^2\pi^2J_K^2 |\omega|^{-2r}} \text{sgn}(\omega) \nonumber \\
G_a(i\omega) &= \frac{-2\rho_0^2\pi^2J_K |\omega|^{-2r}}{1+\rho_0^2\pi^2J_K^2 |\omega|^{-2r}} .
\ea
In addition, the coefficient before each $|\Omega|^{-r}$ should be a function of $r$. However, we only keep the leading order term ($r^0$) since we would set $r=0$ in the final step. Combining Eqs.~\ref{eq:berry_chi}, \ref{eq:bose_chi} and ~\ref{eq:fermi_chi}, we have the quadratic theory of the $\chi$ fields:
\ba
&&S_{\chi} = \int  G^{-1}_\chi(i\Omega)|\chi(i\Omega)|^2 \frac{d\Omega}{2\pi} \nonumber \\
&&G_\chi^{-1}(i\Omega) = i\Omega \bigg[ 1+p(J|\Omega|^{-r})\bigg] +|\Omega|\bigg[ q(J|\Omega|^{-r}) + h|\Omega|^{-\epsilon}\bigg].
\label{eq:Schisupp}
\ea

For future reference, we introduce the series expansion of $G_\chi(i\Omega)$
\ba
G_\chi(i\Omega) &=& \frac{1}{|\Omega|}
\frac{1}{i(\text{sgn}(\Omega)+p(J|\Omega|^{-r})) + q(J|\Omega|^{-r})
} + \frac{1}{|\Omega|}\sum_{m=1}^\infty \frac{(-h)^m|\Omega|^{-m \epsilon} }
{[i(\text{sgn}(\Omega) +p(J)) + q(J)]^{m+1}}
\ea
Since we set $r=0$ in the last step, we drop the $r$ dependence in the summation by setting $r=0$. Including this contribution won't change the beta functions at the leading order. We can further expand the real and imaginary part of $G_\chi(i\Omega)$ as the following
\ba
&&G_\chi(i\Omega) + G_\chi(-i\Omega) =
\sum_{n=0}^\infty a_n J^n |\Omega|^{-1-nr} + \sum_{m=1}^\infty
(-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)
|\Omega|^{-1-m\epsilon} \nonumber \\
&&G_\chi(i\Omega)-G_\chi(i\Omega) =
\sum_{n=0}^\infty ib_n J^n|\Omega|^{1-nr} +   \sum_{m=1}^\infty
(-h)^m \bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)
|\Omega|^{1-m\epsilon}
\ea
where $a_n, b_n$ are the Taylor coefficients of $F(j)$ and $H(j)$ at $j=0$ and
\ba
F(j) &=& \frac{2q(j)}{(1+p(j))^2+q(j)^2}\nonumber \\
H(j) &= &\frac{-2(1+p(j))}{(1+p(j))^2 + q(j)^2}\nonumber \\
c(J) &= & i(1+p(J)) + q(J) \nonumber .
\ea

Now, we combine Eqs.~\ref{eq:corr_exact}, \ref{eq:s_exp}, \ref{eq:simp_exp} and \ref{eq:Schisupp} and calculate the correlation functions of fermions and bosons.
We first expand Eq.~\ref{eq:corr_exact} using Eq.~\ref{eq:s_exp}. For bosons, we have
\ba
&&\langle  \sum_\mu \phi^\mu(\bm{r},\tau=0)^2\rangle_{imp}\nonumber \\ &=&\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }S^2 \nonumber \\
&&+
\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }\langle -S\delta S^z(\tau) -S\delta S^z(\tau') + \delta S^x(\tau) \delta S^x(\tau') + \delta S^y(\tau) \delta S^y(\tau') \rangle_{S_\chi}  + O(S^0) \nonumber \\
\ea
where the last line denotes the contribution of the next-order fluctuations and are evaluated with respect to $S_{\chi}$.
For the fermionic fields, it gives
\ba
&&-\langle  \sum_{a,\sigma} c_{\bm{r_0},a,\sigma}(\tau)
c_{\bm{r_0},a,\sigma}^\dag(0)\rangle \nonumber \\
&=&K\bigg[G_s(\tau) -G_0(\tau) \nonumber \\
&&
+\int_{\tau_1}\sum_\mu \frac{J_K\text{Tr}[G(\tau-\tau_1) \sigma^\mu G(\tau_1)]}{S}
 \langle \delta S^\mu(\tau_1)\rangle_{S_{\chi}}
+\int_{\tau_1,\tau_2}\sum_{\mu,\nu}\frac{J^2_K\text{Tr}[G(\tau-\tau_1)  \sigma^\mu G(\tau_1-\tau_2) \sigma^\nu G(\tau_2)]}{S^2}
\langle
\delta S^\mu(\tau_1) \delta S^\nu(\tau_2)
\rangle_{S_{\chi}}
\bigg]\nonumber \\
&&+ O(S^0)
\ea
where the third line denotes the contribution of the next-order fluctuations and are evaluated with respect to $S_{\chi}$ and $G(\tau) = \text{diag}[\frac{G_s(\tau)+G_a(\tau)}{2},\frac{G_s(\tau)-G_a(\tau)}{2}]$.

To get the final results, we only need to calculate the correlation functions of $\delta S^\mu$ (or $\chi$) that appearing in the above equations. The boson correlation function becomes
\ba
&&\langle  \sum_\mu \phi^\mu(\bm{r},\tau=0)^2\rangle_{imp}\nonumber \\
&=&\int_{\tau,\tau'}\frac{A_\phi}{2\pi } \frac{h}{(\bm{r}^2+\tau^2)^{(2-\epsilon)/2}(\bm{r}^2+\tau'^2)^{(2-\epsilon)/2} }
S^2 \nonumber\\
&&\bigg[
1 + \frac{1}{S}\langle \chi^\dag(\tau)\chi(\tau') + \chi^\dag (\tau')\chi(\tau') - |\chi(\tau)|^2 -|\chi(\tau')|^2\rangle
\bigg]+O(S^0)\nonumber \\
&=&\frac{A_\phi \pi S^2}{2|\bm{r}|^2}h
+
\frac{A_\phi Sh}{2\pi }
\int_0^\infty
2(G_\chi(\Omega)+G_\chi(-\Omega))
\bigg(
2\pi x^{-(1-\epsilon)}B[\frac{\epsilon-1}{2},x]^2|\Omega|^{2-2\epsilon} - \frac{\pi^2}{|\bm{r}|^2}
\bigg)
\frac{d\Omega}{2\pi} \nonumber +O(S^0)\\
&=&\frac{A_\phi \pi S^2}{2|\bm{r}|^2}h
+
\frac{A_\phi Sh}{2\pi }
\int_0^\infty
2\bigg[
\sum_{n=0}^\infty  a_nJ^n|\Omega|^{1-nr}
+\sum_{m=1}^\infty (-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)
|\Omega|^{-1-m\epsilon}
\bigg] \nonumber \\
&&\frac{1}{|\bm{r}|^{2}}
\bigg(
2\pi x^{1-\epsilon}B[\frac{\epsilon-1}{2},x]^2-
\pi 2^{-\epsilon}\Gamma((1-\epsilon)/2)^2
\bigg)
\frac{d\Omega}{2\pi} +O(S^0)\nonumber \\
&=&
\frac{A_\phi \pi S^2}{2|\bm{r}|^2}h
+
\frac{A_\phi S h}{\pi |\bm{r}|^2} \bigg[ \frac{1}{r}
\sum_{n=0}^\infty \frac{-a_nJ^n\pi}{2n}
+\frac{1}{\epsilon}\sum_{m=1}(-h)^m
\bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) \frac{-\pi}{2n\epsilon}
+O(r^0,\epsilon^0)\bigg] +O(S^0) \nonumber \\
&=&\frac{A_\phi \pi S^2}{2|\bm{r}|^2}
h\bigg\{
1- \frac{1}{S}\bigg[\frac{1}{\pi r}\int_0^J\frac{F(j)}{j}dj
+\frac{1}{\pi \epsilon }[\frac{\log(1+h/c(J)}{c(J)} + \frac{\log(1+h/c(J)^*)}{c(J)^*}]
+O(\epsilon^0,r^0)\bigg]+o(\frac{1}{S})
\bigg\}
\label{eq:chi_prop}
\ea
where $x=|\Omega||\bm{r}|$ and $B(a,x)$ is the modified Bessel functions of the second kind. In the above expression, we only keep track of the poles as a function of $\epsilon$ and $r$. We also drop the $\frac{r}{\epsilon}$ poles which goes to zero as we set $r=0$.

For the fermion correlation function,
we first perform Fourier transformation and separate the next-order contributions into four parts $I_1,I_2,I_3,I_4$:
\ba
G_f(i\omega) = -\int_\tau\langle  \sum_{a,\sigma} c_{\bm{r_0},a,\sigma}(\tau)
c_{\bm{r_0},a,\sigma}^\dag(0)\rangle_{imp}
e^{i\omega\tau} = K\bigg[ \frac{2i\pi \rho_0 J^2}{1+J^2} \text{sgn}(\omega)
+I_1+I_2+I_3+I_4+o(\frac{1}{S})\bigg]
\ea
Each part is shown below:
\begin{eqnarray*}
I_1 &=&
\frac{J_K^2}{4S}(G_s(i\omega)^2 +G_a(i\omega)^2)
\int G_\chi(i\Omega) (G_s(i\omega+i\Omega)+G_s(i\omega-i\Omega)-G_s(i\omega)-G_s(i\omega))\frac{d\Omega}{2\pi} \\
&=&\frac{-2i\rho_0\pi J^2(-1+J^2)}{(1+J^2)^2S}
\int_0^\infty \bigg[ \sum_{n=0}^\infty a_nJ^n|\Omega|^{-1-nr} +\sum_{m=1}^{\infty}(-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)
|\Omega|^{-1-m\epsilon}\bigg]\\
&&\sum_{p=0}^\infty (-J^2)^p
\bigg[ \text{sgn}(\omega+\Omega)|\omega+\Omega|^{-2pr-r} +
 \text{sgn}(\omega-\Omega)|\omega-\Omega|^{-2pr-r}
 -2
  \text{sgn}(\omega)|\omega|^{-2pr-r}
\bigg]\frac{d\Omega}{2\pi} +O(r^0)\\
&=&\frac{-i\rho_0\pi J^2(-1+J^2)}{(1+J^2)^2\pi S}\bigg[ \sum_{n=0}^\infty a_nJ^n\sum_{p=0}^\infty (-J^2)^P\frac{-2}{nr}
+
\sum_{m=1}^\infty (-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg) \frac{-2}{m\epsilon}
\bigg] \text{sgn}(\omega)  \\
&=&\frac{2i\rho_0\pi J^2(-1+J^2)}{(1+J^2)^3\pi S}
\bigg[
\frac{1}{r} \int_0^J \frac{F(j)}{j}dj
-\frac{1}{\epsilon}\bigg(
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
+O(r^0,\epsilon^0)\bigg)
\bigg] \text{sgn}(\omega) .
\end{eqnarray*}

\begin{eqnarray*}
I_2 &=&
-\frac{J_K^2}{2S}G_s(i\omega)G_a(i\omega)
\int G_\chi(i\Omega)
(G_a(i\omega+i\Omega)+G_a(i\omega-i\Omega)-G_a(i\omega)-G_a(i\omega))\frac{d\Omega}{2\pi} \\
&=&-\frac{2iJ^3(-2J\rho_0\pi)}{(1+J^2)^2S}\text{sgn}(\omega)
\int_0^\infty
\bigg[
\sum_{n=0}^\infty a_nJ^n|\Omega|^{-1-nr} +\sum_{m=1}^{\infty}(-h)^m \bigg(\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)
|\Omega|^{-1-m\epsilon}
\bigg] \\
&&
\sum_{p=0}^\infty
(-J^2)^p \bigg[
|\omega+\Omega|^{-2pr-2r}
+|\omega-\Omega|^{-2pr-2r}
-2|\omega|^{-2pr-2r}
\bigg]
\frac{d\Omega}{2\pi}\\
&=&-\frac{2iJ^3(-2J\rho_0\pi)}{2\pi (1+J^2)^2S}\text{sgn}(\omega)
\bigg[
\sum_{n=0}^\infty a_n
\sum_{p=0}^\infty (-J^2)^p
J^n(\frac{2}{2pr+2r+nr} - \frac{2}{nr})
+O(r^0,\epsilon^0)
\bigg]
 \\
&=& \frac{2iJ^4\rho_0\pi}{\pi (1+J^2)^2S}\text{sgn}(\omega)
 \bigg[
 \frac{1}{rJ^2}\int_0^J \frac{2jF(j)}{1+j^2}dj
 -\frac{2}{r(1+J^2)}\int_0^J\frac{F(j)}{j}dj
 +O(r^0,\epsilon^0)\bigg].
\end{eqnarray*}

\begin{eqnarray*}
I_3 &=&\frac{J_K^2}{4S}
\bigg[
G_a(i\omega)^2 +G_s(i\omega)^2
\bigg] \int G_\chi(i\Omega) [
G_a(i\Omega+i\omega) - G_a(i\omega-i\Omega)]
\frac{d\Omega}{2\pi} \\
&=& \frac{J^2(-1+J^2)}{2\pi(1+J^2)^2S}
\int_0^\infty \bigg[
\sum_{n=0}ib_nJ^n|\Omega|^{1-nr} +\sum_{m=1}(-h)^m
\bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)|\Omega|^{-1-m\epsilon} \bigg]
\\
&&\bigg[ -2J\rho_0\pi \sum_{p=0}^\infty (-J^2)^p
\bigg(
|\Omega+\omega|^{-2r-2pr} -
|-\Omega+\omega|^{-2r-2pr}
\bigg)
\bigg]
d\Omega \\
&=&\frac{1}{S}[0+O(r^0,\epsilon^0)].
\end{eqnarray*}

\begin{eqnarray*}
I_4 &=&
\frac{-J_K^2}{2}G_a(i\omega)G_s(i\omega)
\int G_\chi(i\Omega) \bigg(
G_s(i\Omega+i\omega)-G_s(i\omega-i\Omega)
\bigg) \frac{d\Omega}{2\pi} \\
&=&\frac{-2J^3\rho_0\pi}{S(1+J^2)^2}\int_0^\infty
\bigg[
\sum_{n=0}ib_nJ^n|\Omega|^{1-nr} +\sum_{m=1}(-h)^m
\bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)|\Omega|^{-1-m\epsilon} \bigg]\\
&&
\sum_{p=0}(-J^2)^p
\bigg[
\text{sgn}(\omega+\Omega)|\omega+\Omega|^{-2pr-r}
-\text{sgn}(\omega-\Omega)|\omega-\Omega|^{-2pr-r}
\bigg]
\frac{d\Omega}{\pi} \\
&=&\frac{-2J^3\rho_0\pi}{S(1+J^2)^2}
\bigg[
\sum_{n=0}ib_nJ^n\sum_{p=0} (-J^2)^p \frac{2}{\pi(nr+2pr+r)} +\sum_{m=1}(-h)^m\frac{1}{1+J^2}
\bigg(-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg)\frac{2}{\pi m\epsilon} \bigg]\text{sgn}(\omega)  \\
&=&\frac{-2J^3\rho_0\pi}{(1+J^2)^2}
\bigg[ \frac{2}{ r\pi J}\int_0^J\frac{iH(j)}{(1+j^2)}dj
+\frac{2}{\pi \epsilon(1+J^2)}\bigg(
-\frac{\log(1+h/c(J))}{c(J)}
+\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) +O(r^0,\epsilon^0)
\bigg]\text{sgn}(\omega)
.
\\
\end{eqnarray*}
In the above derivation, we only keep terms that diverge as $r\rightarrow 0$ and $\epsilon \rightarrow 0$, with $r$ going to zero first.

Summing over all the contributions, we have
\ba
&&G_f(i\omega) \nonumber \\
&=& 2iK\pi \rho_0\text{sgn}(\omega) \frac{J^2}{1+J^2}
\bigg[
1
+ \frac{-1+J^2}{(1+J^2)^2\pi Sr}\int_0^J\frac{F(j)}{j}dj
+ \frac{1}{(1+J^2)\pi Sr}\int_0^J\frac{2jF(j)}{1+j^2}dj \nonumber\\
&& -\frac{2J^2}{(1+J^2)^2\pi Sr}
\int_0^J\frac{F(j)}{j}dj
-\frac{2}{(1+J^2)r\pi S}
\int_0^J \frac{H(j)}{1+j^2}dj \nonumber \\
&&-\frac{-1+J^2}{(1+J^2)^2 \pi S\epsilon}
\bigg(
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg)
+\frac{2iJ}{(1+J^2)^2\pi \epsilon S}
\bigg( -
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg)
\bigg]
\nonumber \\
\label{eq:g_prop}
\ea

In summary, we have
\ba
\langle \sum_\mu [ \phi^\mu(\bm{r},\tau=0)]^2\rangle_{imp}
&\propto &  \frac{1}{|\bm{r}|^2}
h\bigg[
1-
\frac{1}{S\pi r}\int_0^J\frac{F(j)}{j}dj +\frac{1}{S\pi \epsilon }[\frac{\log(1+h/c(J)}{c(J)} + \frac{\log(1+h/c(J)^*)}{c(J)^*}]
\bigg] \nonumber \\
\sum_{\sigma,a}\langle c_{\sigma,a}(\bm{r_0},\tau) c_{\sigma,a}^\dag(\bm{r_0},\tau)\rangle_{imp} &\propto&
\frac{J^2\text{sgn}(\tau)}{(1+J^2)|\tau|}
\bigg[
1
+ \frac{-1-J^2}{(1+J^2)^2\pi Sr}\int_0^J\frac{F(j)}{j}dj
+ \frac{1}{(1+J^2)\pi Sr}\int_0^J\frac{2jF(j)}{1+j^2}dj \nonumber\\
&&
-\frac{2}{(1+J^2)r\pi S}
\int_0^J \frac{H(j)}{1+j^2}dj -\frac{-1+J^2}{(1+J^2)^2 \pi S\epsilon}
\bigg(
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg) \nonumber \\
&&
+\frac{2iJ}{(1+J^2)^2\pi \epsilon S}
\bigg( -
\frac{\log(1+h/c(J))}{c(J)}
+
\frac{\log(1+h/c(J)^*)}{c(J)^*}
\bigg)
\bigg] . \nonumber
\ea

\subsection{C. Beta function analysis}
We numerically evaluate the beta functions and the resulting RG flows are shown in Fig.~\ref{fig:rg_num}.
Besides numerical evaluations, we analyze the beta equations in certain extreme cases.

At $h=0$, the model becomes a standard Kondo impurity model. The beta function of $J$ is
\ba
\beta_J\bigg|_{h=0} &=& \frac{J-J^3}{2S\pi (1+J^2)}F(J) + \frac{J^2}{S\pi(1+J^2)}H(J)  \nonumber \\
&=&-\frac{2J^2}{S\pi} +\frac{2\kappa J^3}{S\pi^2}  +O(J^4/S).
\ea
Numerically, we find $J$ is always relevant at $\kappa=2$ and the system flows to a strong coupling fixed point $K$ at $J=\infty$. The first two terms in the second line are also consistent with the result obtained from perturbation theory as shown in Ref.~\cite{affleck1995conformal} (note that we've scaled the Kondo coupling).

At $J=0$, the model becomes a Bose-Kondo model with the following beta function of $h$:
\ba
\beta_h\bigg|_{J=0} = -\frac{1}{S}\tilde{\epsilon} h + \frac{2h^2}{S\pi (1+h^2)}.
\ea
This is consistent with previous results obtained in Refs.~\cite{cuomo2022spin,nahum2022fixed}. It gives two stable fixed points at $h = (1-\sqrt{1-\tilde{\epsilon}^2})/ {\tilde{\epsilon}},\infty$ (L and L$'$) and one critical point at $h=(1+\sqrt{1-\tilde{\epsilon}^2})/ {\tilde{\epsilon}}$ (LC). After introducing the Kondo coupling, we find $J$ is irrelevant at both L and LC, but marginal at L$'$.

At $h=\infty$, we have a line of fixed points (L$'$) where $J$ is always marginal. The beta function at the large $h$ is:
\ba
\beta_J = \frac{J-J^3}{\pi S (1+J^2)}\frac{1}{h} + O(\frac{1}{h^2S}).
\ea
Clearly, $\beta_J \rightarrow 0 $ as $h\rightarrow \infty$.

At $J=\infty$, besides the K fixed point, there is an unstable fixed point KD$'$ at $(h_c=\frac{2}{\tilde{\epsilon}\pi},J=\infty)$. To see the existence of such a critical point, we give the beta function of $h$ at large $J$:
\ba
\beta_h= -\frac{\tilde{\epsilon} h}{S}+\frac{2}{S\pi} + O(\frac{1}{SJ}).
\ea
It suggests that $h$ flows to $0$ ($\infty$) when $h<(>) h_c$ at $J\rightarrow \infty$.

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{rg.pdf}
    \caption{Renormalization group flow at $\tilde{\epsilon}<\tilde{\epsilon}_1$ (a), $\tilde{\epsilon}_1<\tilde{\epsilon}<\tilde{\epsilon}_2$ (b) and $\tilde{\epsilon}_2<\tilde{\epsilon}$ (c) obtained from numerical evaluation. For purpose of illustration, we normalize all the vectors to have the same length. }
    \label{fig:rg_num}
\end{figure*}

\subsection{D. Topological contribution from the Kondo coupling}
In this section, we'll show that the Kondo coupling
contributes an additional topological term to the impurity spin at sufficiently large values of $J_K$. For the purpose of discussion, we consider pure Kondo model, in other words, the BFK model at $g=0$ or $h=0$. We then integrate out the fermionic fields, which gives the following action
\ba
S_{K}
=-iSS_B
-K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg).
\ea
We now evaluate the fermion determinant contribution at sufficiently large $J_K$ via  a Trotter decomposition
\ba
\exp\bigg ( K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg) \bigg) = \bigg[\lim_{N\rightarrow \infty }\text{Tr}[\prod^N_i \exp[-H(\tau_i) d\tau]]\bigg]^K
\ea
where
\ba
H(\tau_i) = J_K\sum_\mu c_{\bm{r}_0}\frac{\sigma^\mu}{2} c_{\bm{r}_0}S^\mu(\tau_i) + \sum_{\bm{k} \sigma } \epsilon_{\bm{k}}c_{\bm{k},\sigma}^\dag c_{\bm{k},\sigma}.
\ea
Since the $K$ orbitals are
degenerate, it's sufficient to consider one orbital and multiply the contribution by $K$. When $J_K$ is sufficiently large, the hopping between site $\bm{r}_0$ and the other sites is a high order effect ($\frac{t^2}{J_K}$). Here, we focus on the leading order contribution, and it's sufficient to take $H(\tau_i) = J_K/2\sum_\mu c_{\bm{r}_0}\sigma^\mu c_{\bm{r}_0}S^\mu(\tau_i)$. The eigenstates and eigenvalues of $H(\tau_i)$ are
\ba
&&|v_{+}(\tau_i)\rangle = \cos(\frac{\theta(\tau_i)}{2})c_{\bm{r_0},\up}^\dag
|0\rangle +\sin(\frac{\theta(\tau_i)}{2})e^{i\phi(\tau_i)}c_{\bm{r_0},\dn}^\dag |0\rangle \nonumber \\
&&|v_{-}(\tau_i)\rangle = -\sin(\frac{\theta(\tau_i)}{2})c_{\bm{r_0},\up}^\dag
|0\rangle +\cos(\frac{\theta(\tau_i)}{2})e^{i\phi(\tau_i)}c_{\bm{r_0},\dn}^\dag |0\rangle \nonumber \\
&&|v_0(\tau_i)\rangle = |0\rangle \nonumber  \\
&&|v_d(\tau_i)\rangle = |\up\dn\rangle \nonumber \\
&&E_{\pm}(\tau_i) = \pm \frac{J_KS}{2},\quad E_0(\tau_i) = 0 , \quad E_d(\tau_i) =0
\ea
where we let $\bm{S}(\tau_i) =S (\sin(\theta(\tau_i))\cos(\phi(\tau_i)), \sin(\theta(\tau_i))\sin(\phi(\tau_i)), \cos(\theta(\tau_i)))$. We have the following decomposition via the eigenvalues and eigenstates
\ba
e^{-H(\tau_i)d\tau} &=& e^{-E_+(\tau_i)d\tau}|v_{+}(\tau_i)\rangle\langle v_{+}(\tau_i)|+ e^{-E_-(\tau_i)d\tau}|v_{-}(\tau_i)\rangle\langle v_{-}(\tau_i)|
\nonumber \\
&&
+e^{-E_0(\tau_i)d\tau} |v_0(\tau_i)\rangle \langle v_0(\tau_i)|
+e^{-E_d(\tau_i)d\tau} |v_d(\tau_i)\rangle \langle v_d(\tau_i)|
.
\ea
Combining the above identity with a Trotter decomposition of the fermion determinant, we obtain the following leading order contribution
\ba
&&\exp(K\text{Tr}\log\bigg( \delta_{\tau,\tau'}\delta_{k,k'}(\partial_\tau  + \epsilon_{\bm{k}})- \delta_{\tau,\tau'}J_K\sum_\mu S^\mu(\tau)\frac{\sigma^\mu}{2}\bigg)) \nonumber \\
&\propto& \exp\bigg(K\int_\tau \langle \partial_\tau v_-(\tau)|v_-(\tau)\rangle d\tau +   O(\frac{1}{J_K}) \bigg)\nonumber \\
&\propto& \exp\bigg(i\frac{K}{2}\int_\tau (1-\cos(\theta))\partial_\tau \phi d\tau +   O(\frac{1}{J_K}) \bigg)\nonumber \\
&=& \exp\bigg( -i\frac{K}{2}S_B+
O(\frac{1}{J_K})\bigg).
\ea
The effective action  then
becomes
\ba
S_{K}
= -i(S-\frac{K}{2})S_B + O(\frac{1}{J_K})
\ea
which indicates the non-trivial topological contribution from the coupling to the gapless fermions.
Qualitatively,
 the spin has been Kondo-screened, with the effective spin becoming $S-\frac{K}{2}$.
 {For the}
 perfect screening case
 {that we consider}, $S-\frac{K}{2}=0$, and we obtain a Kondo singlet.

\end{document}