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

Among the populations of compact astrophysical objects, strongly magnetized isolated neutron stars, or magnetars \citep{DT92}, are by far the most vibrant. Their extreme magnetic fields (10$^{14}$-10$^{15}$ G) are likely powering the underlying mechanisms of the large variety of their observed characteristics \citep{TD95,TD96}. Since the discovery of magnetar-like field strengths with the Rossi X-ray Timing Explorer \citep{CK1998}, long-term monitoring of their bright, persistent, soft X-ray emission revealed new characteristics, such as excessive timing noise \citep{Woods2000} and timing anomalies \citep[glitches,][ and anti-glitches, \citealt{ach13,younes20d}]{kaspi00}. However, their persistent emission is not limited to only soft X-rays. Roughly a dozen magnetars have exhibited emission in hard X-rays (above 15 keV), and a handful of them have been observed in optical and/or radio wavelengths. See \cite{kasbel17} for a recent review covering their broadband properties, and \cite{olausen} for a dynamic web portal listing the properties of currently known magnetars.

What remains the most intriguing property of some magnetars is the repeated emission of multiple, highly energetic bursts. Magnetar bursts appear in broad morphological and energetic varieties; these range from sub-second-long typical short bursts at super-Eddington luminosities (10$^{39}$--10$^{41}$\,erg s$^{-1}$) to giant flares that last up to several minutes, with energies in excess of 10$^{45}$\,erg. In between these classes is a group of bursts, each typically with a duration of a few to a few tens of seconds and energies of around 10$^{42}$\,erg, which are commonly referred to as intermediate events \citep{Kouveliotou2001}. Such events have been seen from SGR 1627$-$41 \citep{mazets99}, SGR 1900+14 \citep{ibrahim,lenters,israel08}, SGR 1806$-$20 \citep{gogus11}, and SGR J1550$-$5418 \citep{mereg09,kaneko10}. These events usually occur at the onset, or in the early phases, of their burst active episodes. Spectrally, the intermediate events exhibit thermal characteristics, represented by a blackbody temperature of a few keV, which decays over a timescale of seconds to minutes. Pulsations are commonly observed during the tail of the decaying emission, with a higher pulsed intensity compared to the level prior to the intermediate event. The case of SGR J1550$-$5418 was particularly interesting, as that burst ``storm" occurred several hours before the intermediate events \citep{mereg09}. At the onset of this storm, an enhancement of the underlying persistent emission in hard X-rays lasted $\sim2$ minutes. The spectrum of this enhanced emission was well described with a blackbody temperature of $kT=17$\,keV, accompanied with pulsed emission detected up to $\sim$70~keV \citep{kaneko10}. The enhanced emission was energetically comparable to that of emission from tails of typical intermediate events, and was interpreted as emission originating from a short-lived hot spot on the surface of the neutron star.

\mbox{SGR\,J1935+2154~} became active again in 2020 April when it emitted hundreds of bursts \citep{lin20}. The onset of this episode was marked with exceptional activity; a storm of highly energetic bursts detected with the Swift/Burst Alert Telescope  \citep[BAT:][]{palmer20} and the Fermi/Gamma-ray Burst Monitor \citep[GBM;][]{fletcher}. Follow-up observations of \mbox{SGR\,J1935+2154~} with Neutron star Interior Composition ExploreR (NICER) captured another cluster of bursts $\sim$6 hours later \citep{younn20}, while INTEGRAL, KONUS-WIND and HXMT captured the X-ray counterpart of the first observed Galactic Fast Radio Burst \citep{mereghetti20,ridnaia2021,Li2021}. A careful inspection of the Fermi/GBM data showed the onset of the event resembled the 2006 ``burst forest" observed from SGR\,$1900+14$, during which more than 40 bursts were detected with Swift/BAT over $\sim30$\,s \citep{israel08}. Nearly one-third of those bursts exhibited highly unusual morphology with durations of $\sim$1\,s, longer than typical magnetar bursts, some of which were flat-topped \citep{israel08}. We also note that a similar burst forest from SGR\,1806$-$20 was observed with INTEGRAL, during which $>$100 bursts were observed within $\sim$600\,s \citep{gotz06}.

\citet{lin20} analyzed time-integrated spectra of 125 Fermi/GBM bursts from the 2020 April activity of \mbox{SGR\,J1935+2154~}, excluding the 130\,s segment, where a cluster of such peculiar events were observed, hereafter referred to as the ``burst forest." Here we present the results of our time-resolved spectral analyses of this 130\,s long burst forest, observed with Fermi/GBM. The observation and the data analysis methods are described in \S{\ref{sec:data}}, followed by the presentation of the analysis results in \S{\ref{sec:res}}.  We interpret and discuss the results in \S{\ref{sec:disc}}. We note that the estimated distance to \mbox{SGR\,J1935+2154~} ranges from $\sim$2 to 9\,kpc \citep[see][and the references therein]{zhong20} including the most recent estimations from the dust-scattering halo of 4.5$^{+2.8}_{-1.3}$\,kpc \citep{mereghetti20} and from the host supernova remnant of 6.6$\pm$0.7\,kpc \citep{zhou20}; here we use the distance of 9\,kpc, consistent with \citet{lin20}.


\section{Observation and Data Analysis} \label{sec:data}
GBM consists of 12 NaI detectors (covering 8\,keV-1\,MeV) and 2 BGO detectors ($\sim$200\,keV-30\,MeV). On 2020 April 27, three of the NaI detectors (\#6, 7, \& 9) were triggered at 18:26:20 UT (609704785.155 MET = $T_0$) by a burst originating from the direction of \mbox{SGR\,J1935+2154~} (trigger\# bn200427768).  Upon inspection of the event data, we identified the burst forest from $T_0 + 310$\,s to $T_0 + 430$\,s (see Figure~\ref{fig:lc_2segs}), which included several broad peaks with a few-seconds duration. The GBM triggered data are accumulated in several different formats.  Of those, the time-tagged event (TTE) datatype provides the finest resolution, both in time  (2\,$\mu$s) and in energy (128 energy channels).  Therefore, we used the TTE data with the time-dependent detector response matrices (rsp2) provided by the GBM team. For this work, we did not include the BGO data as adding the BGO data did not affect the spectral parameters at all due to the low count rates above a few hundred keV.

\begin{figure*}
\begin{center}
\includegraphics[scale=1, trim=45 220 30 320, clip]{lightcurve_bothseg_bkg.pdf}
\end{center}
\caption{
Lightcurve of the entire ``burst forest" episode ($T_0 + 310$\,s to $T_0 + 430$\,s, 8 ms resolution). The count rates are summed over all triggered detectors. The blue dashed line shows the separation between Segments 1 and 2. The yellow-hatched regions are the background time intervals used to determine the background rate (red dashed line).  
\label{fig:lc_2segs}}
\end{figure*}

\subsection{Detector selection and background modeling}
Because the detector-to-source angles change considerably during 130\,s episode, we divided this data into two segments separated at $T_0 +392$\,s, taking into account the detector angles and the source lightcurve. We also calculated the field of view and spacecraft blockage (mainly by the Large Area Telescope) for each NaI detector during this time period. We used detectors having detector-to-source angles $< 55^{\rm o}$ that were unblocked for each of the two segments (1.~$T_0+310$\,s to $T_0+392$\,s and 2.~$T_0+392$\,s to $T_0+433$\,s). These requirements constrained the detectors used for the analysis to NaI\,6, 7, and 9 for \mbox{Segment\,1~} and NaI\,0, 1, 6, and 9 for Segment\,2.

To model the background count rates, we first generated a Bayesian-block representation of the lightcurve \citep{scargle13} using the TTE data of the triggered detectors (binned to 8\,ms, 10--100\,keV) and defined the background time intervals as the blocks longer than 4\,s (hatched regions in Figure\,\ref{fig:lc_2segs}). The background count rates were then interpolated using the 8 ms lightcurve of the background time intervals in the energy range used in the spectral analysis (8--500\,keV), by fitting a second-order polynomial function. We also cross checked our choice of the background time intervals by running the untriggered event search, which uses Poisson statistics to identify count rate increases, over the entire burst-forest period; none of the untriggered bursts identified were within our background time intervals.

\subsection{Evaluation of Deadtime and Pile-up}
Before proceeding with the spectral analysis, we evaluated the possibility of count saturation in the data due to the very high count rates recorded for some of the burst forest peaks. The TTE data suffers from a deadtime of 2.6\,$\mu$s due to the fixed data-packet processing speed of 375~kHz on the spacecraft~\citep{meegan09}. Consequently, bandwidth saturation of the data occurs when the collective count rates of all 14 GBM detectors (12 NaI + 2 BGO), exceed this limit. For the burst forest episode, the maximum count rate summed for all 14 detectors over all energy bands is $\sim250$\, kHz, which is well below this saturation limit.
Additionally, we also evaluated the data for pulse pile-up effects, a phenomenon that inaccurately processes normally discrete spectral and temporal values when the ``recovered" count rates exceed $\sim$10$^5$ counts per second per detector \citep{chaplin13, bhat14}. The deadtime-corrected TTE count rates for the detectors used in this study do not exceed this limit, and thus we conclude that our data is unaffected by this behavior and our analysis of the burst forest is accurate.

\subsection{Spectral Analysis}
We performed time-resolved spectral analysis  ($8-$500\,keV) using the RMFIT spectral analysis tool (version 4.3.2).\footnote{https://fermi.gsfc.nasa.gov/ssc/data/analysis/rmfit} The energy range of 30-40\,keV was excluded in our analysis to avoid possible contribution to the fit statistics due to the iodine $K$-edge at 33\,keV. This exclusion, however, did not affect the analysis outcome.
We fit each of the two time segments with four models commonly used to describe spectra of magnetar bursts: a power law with exponential cutoff (COMPT), a single blackbody (BB), a double blackbody (BB+BB), and a blackbody plus power law (BB+PL). The fit parameters were obtained by minimizing the log-likelihood (C-stat).  The spectral fits were performed in various time resolutions: 4, 16, and 128 ms.  The 4 ms and 16 ms spectral bins were created with a signal-to-noise (S/N) minimum of 25 in the brightest detector. We present in Table\,\ref{tab:param} the spectral parameters obtained with either COMPT, BB, or BB+BB fits for each time interval. We note here that based on the spectral fit results, we identified two distinct kinds of intervals, which we define as peaks and troughs. The former are invariably associated with the higher photon fluxes and can be best fit with either COMPT or BB+BB. The latter are associated with flux minima and are best fit with either COMPT or a single BB. Because the COMPT model provides reasonable fits for all spectra of the burst forest, we use the BB+BB and BB fits to define the peak/trough definition. 

Throughout this Letter we have consistently used equal time bins of 128\,ms for the peaks in all the figures, unless otherwise stated. However, for troughs, the requirement of a 25 S/N results in bins with variable time resolution. All these parameters are presented in Table\,\ref{tab:param}. Finally, we find that the fluence during the entire duration of the burst forest was $F = (2.305 \pm 0.001)\times 10^{-4}$\,erg\,cm$^{-2}$ in the $8–500$ keV energy range; Segment 1 accounted for $(0.562 \pm 0.002) \times 10^{-4}$\,erg\,cm$^{-2}$ and Segment 2 for $(1.743\pm 0.002) \times 10^{-4}$ \,erg\,cm$^{-2}$. 

\section{Results}\label{sec:res}
We find that of the four models mentioned in Section 2.3, the parameters for three of these (COMPT, BB, and BB+PL) are reasonably constrained for {\it all} time-resolved spectra.  When we fit the BB+BB model to the troughs, it often fails to constrain the high-energy BB parameters. Finally, the BB+PL model provides the least acceptable statistics and has not been used in further analysis (see also Section\,\ref{sec:bic}). When fitting the COMPT model, we find that the $\nu F_{\nu}$ spectral peak ($E_{\rm peak}$) varies significantly within each individual temporal peak, while the low-energy spectral index remains constant around 0.5 (see Figure\,\ref{fig:compt_param}). The $E_{\rm peak}$ generally follows the lightcurve, and we find a significant positive correlation between $E_{\rm peak}$ and photon flux; the Spearman rank-order correlation coefficients are $r_s = 0.4$ and 0.9, with $P < 10^{-6}$ and $10^{-17}$ for the \mbox{Segment\,1~} and Segment 2, respectively.  We note that the COMPT model is defined as $f(E) \propto (E/E_{\rm peak})^{\Gamma}$; therefore, $\Gamma > 0$ indicates that the spectral curvature below $E_{\rm peak}$ resembles that of thermal spectra.  The thermal nature of these spectra is also supported by the fact that BB+PL fits provide well-constrained parameters: the PL index remains constant (weighted average of $-2.35 \pm$0.07 for both segments combined), representing the spectra above $\sim$100\,keV, and the BB $kT$ follows the $E_{\rm peak}$ evolution.

\begin{figure*}
\begin{center}
\epsscale{1.16}
\plottwo{seg1part1_compt_param_mod}{seg1part2_compt_param_mod}
\includegraphics[scale=0.5, trim=0 150 0 50, clip]{seg2peaks_compt_param_128ms}
\end{center}
\caption{COMPT-fit parameters ($E_{\rm peak}$ and index $\Gamma$) for \mbox{Segment\,1~} (top two panels) and for the triple peaks of \mbox{Segment\,2~} (bottom panel; 4 ms resolution binned to S/N $\geq$ 25 for troughs and 128 ms resolution for peaks).  The count rates (right y-axis) are plotted with gray dotted lines. The red histograms at the top of each panel show the extrapolated NICER phase cycles; the cycle numbering starts at the GBM trigger time (see \S\ref{sec:disc}).
\label{fig:compt_param}}
\end{figure*}

When we fit the BB+BB model, we find that within the Segment 1 peaks and troughs, the low and high $kT$ do not vary drastically; however, for Segment 2 peaks the $kT$s vary significantly. The $kT$s are also strongly correlated with each other in all peaks in Segment\,2, with $r_s \sim 0.8$ and $P < 10^{-23}$.  In Figures\,\ref{fig:bbbb_kt_1} and \ref{fig:bbbb_kt_2}, we present the $kT$ evolution during the two segments.

\subsection{Spectral Fit Statistics}\label{sec:bic}

The fit statistics that we employed for our spectral fitting (C-stat) does not provide the goodness of fit and cannot be directly used to identify the model that better describes each of the spectra.  Therefore, we calculate for each spectral fit the Bayesian Information Criterion  \citep[BIC;][]{liddle07} as follows:

\begin{equation}
  {\rm BIC} = -2 \ln \mathcal{L}_{\rm max} + k \ln N = {\rm C\mhyphen stat} + k \ln N , 
\end{equation}

\noindent where $\mathcal{L}_{\rm max}$ is the maximum likelihood, $k$ is the number of parameters in the photon model and $N$ is the number of data points.  We then compare the BIC values associated with each of the four model fits in pairs and the model with significantly lower BIC, i.e. with $\Delta {\rm BIC} > 10$, was considered preferred.
According to our $\Delta {\rm BIC}$ comparisons, we find that BB+BB for the peaks and BB for most of the troughs are consistently identified as preferred models.  In Figure~\ref{fig:bic_comparison}, we show the preferred models for all spectra in  Segment\,2.  We see a clear transition of the preferred models from BB+BB to BB in the decaying parts of the broadest peaks.

Finally, within the trough spectra where BB is preferred, the $kT$ does not vary significantly. The weighted means of the BB $kT$ in the Segment\,1 and 2 trough spectra are 8.0$\pm$0.5\,keV and 7.8$\pm$0.3\,keV, respectively.  
We assume, therefore, that this component is manifested as the low-energy $kT$ in the BB+BB spectra, and fix the low-energy $kT$ to the weighted mean values mentioned above. We then
proceed to study the evolution of the  high-energy $kT$.  Although in some parts of the peaks the high-energy $kT$ is then found to be slightly lower ($\lesssim 2$\,keV) than those found with the varying low-energy $kT$, we find that the overall high-$kT$ evolution does not change significantly from what is seen in Figure\,\ref{fig:bbbb_kt_2}.

\begin{figure*}
\begin{center}
\epsscale{1.16}
\plottwo{sgr1935_storm1_pulse_MOD2.pdf}{seg1part2_bbkt_param_mod.pdf}
\end{center}
\vspace{-2cm}
\caption{Blackbody $kT$ evolution during Segment 1 (4 ms resolution binned to S/N $\geq$ 25 for troughs and 128 ms resolution for peaks). $kT_{\rm high}$ is shown in blue for the intervals where BB+BB better represents the spectra over a single BB.  $kT_{\rm low}$ and single BB $kT$ are shown in black. The count rates (left y-axis) and the temperatures (right y-axis) are shown vs. barycentered times along with the phases of the persistent emission pulse profiles extrapolated using NICER data (top panels).
\label{fig:bbbb_kt_1}}
\end{figure*}

\begin{figure*}
\begin{center}
\plotone{sgr1935_storm2_pulse_128ms.pdf}
\end{center}
\vspace{-4cm}
\caption{Blackbody $kT$ evolution during \mbox{Segment\,2~} (4 ms resolution binned to S/N $\geq$ 25 for troughs and 128 ms resolution for peaks). $kT_{\rm high}$ is shown in blue for the intervals where BB+BB better represents the spectra over a single BB.  $kT_{\rm low}$ and single BB $kT$ are shown in black.  The count rates (left y-axis) and the temperatures (right y-axis) are shown vs. barycentered times, along with the phases of the persistent emission pulse profiles extrapolated using NICER data (top panels). The yellow shaded intervals indicate the intervals where BB+BB is consistently preferred over BB (see Figure \ref{fig:bic_comparison} ).
\label{fig:bbbb_kt_2}}
\end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[scale=0.6, trim=0 220 0 0, clip]{seg2_preferredmodel.pdf}
\end{center}
\vspace{-3.5cm}
\caption{Segment\,2 Photon flux history showing the preferred model between BB and BB+BB for each of the time intervals (4-ms resolution binned to S/N $\geq$ 25) determined based on the differences in BIC (i.e., $\Delta {\rm BIC = BIC_{BB} - BIC_{BB+BB}} >$ 10 if BB+BB is preferred).  The y-axis values for BB+BB or single BB points were arbitrarily selected. The hatched regions define the ``peaks" in our analysis, as discussed above. The photon flux values are estimated with the fitted models.
\label{fig:bic_comparison}}
\end{figure*}

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

\citet{younes20} presented the pulse profile of the  persistent emission of \mbox{SGR\,J1935+2154~} as observed with NICER $\sim$6\,hours after the GBM burst forest episode discussed here.  We converted the GBM photon arrival times during the burst forest to barycentric times and matched the episode lightcurves and their spectral parameter evolution with the extrapolated NICER pulse profile.  The extrapolated pulse profile is shown in the top panels of Figures\,\ref{fig:compt_param}, \ref{fig:bbbb_kt_1}, \ref{fig:bbbb_kt_2}. In particular, we observe from the one-component COMPT model fits (Figure\,\ref{fig:compt_param}) that the $E_{\rm peak}$ values of the spectra appear to anti-correlate with the pulse profile, with the spectra peaking at higher energies at or close to the minima of the pulse profile.  Furthermore, we find that in Segment\,2, the transition from peaks to troughs occurs roughly at the spin period ($P_{\rm spin}=3.24$\,s) of the source (see hatched regions in Figure\,\ref{fig:bbbb_kt_2}). Motivated by this trend, we superimpose the COMPT $E_{\rm peak}$ evolution of each pulse phase within \mbox{Segment\,2~} to investigate the phase-resolved spectral evolution, which is shown in Figure\,\ref{fig:compt_ep_phase}.  It is evident that for all three peaks, the $E_{\rm peak}$ stays high ($\sim$45\,keV) during the phase minimum and transitions to lower energies ($\sim$25--30\,keV) at the phase maximum.  It is also clear from the figure that right before and right after the triple peaks in the lightcurve (i.e., beginning of Phase~125 and the end of Phase~128), the spectra do not follow this trend. We also note that there are some time intervals in the trough parts of the lightcurve (i.e., non-hatched regions) for which the spectra are still better described with the BB+BB model, based on the BIC comparisons.

We present in Figure\,\ref{fig:flux_rad_segs} the the energy flux versus the emitting radius of the BB+BB and BB components. We note that especially for Segment 2, the single BB $kT$ (i.e., the trough spectra) remains constant around 7\,keV, while the flux varies around two orders of magnitude.  The $kT_{\rm high}$ also remains roughly constant at around 14\,keV throughout both segments. The high-energy BB corresponds to smaller emitting radii than the low-energy component, while the ranges of the energy fluxes of the four peaks are comparable. In Figure\,\ref{fig:r2_kt}, we present the estimated BB emitting area ($R^2$) as a function of temperature, for all spectra during the 130 s episode.  The $kT_{\rm low}$ and $kT_{\rm high}$ of the peaks represented by BB+BB models are negatively correlated with the emitting area ($r_s = -0.79$, $P < 10^{-20}$). In the same figure, we also show as a reference the luminosity contours calculated with the Stefan-Boltzmann law, assuming a distance to the source of 9 kpc. It is clear that the luminosity values remain below $\sim10^{41}$\,erg\,s$^{-1}$, which corresponds to the upper limit of magnetar short-burst luminosity.
We also note that the distribution of $kT_{\rm low}$ in the burst forest episode is narrower than the that of other short bursts from the same active episode presented in \citet{lin20}.

\begin{figure*}
\begin{center}
\includegraphics[scale=0.7]{seg2_compt_ep_phase_128ms.pdf}
\end{center}
\vspace{-6cm}
\caption{COMPT $E_{\rm peak}$ values of \mbox{Segment\,2~} peaks (Figure\,\ref{fig:compt_param}, bottom panel) as a function of phase. The NICER pulse phase is extrapolated to the Segment 2 data, shown with black dashed lines (right axis; NICER pulse count rate, also shown in Figure\,\ref{fig:bbbb_kt_2}) The colored histograms correspond to the evolution of $E_{\rm peak}$ in different phases.
\label{fig:compt_ep_phase}}
\end{figure*}

\begin{figure*}
\begin{center}
\epsscale{1.16}
\plottwo{seg1_bbflux_vs_radsq_LOG_gray.pdf}{seg2_bbflux_vs_radsq_LOG_gray.pdf}
\end{center}
\vspace{-3cm}
\caption{BB \& BB+BB flux vs.~emitting areas ($R^2$) for \mbox{Segment\,1~} and \mbox{Segment\,2~} separately. Lines corresponding to constant temperatures, $kT = 7$\,keV \& 14\,keV are shown as a reference in both panels.
\label{fig:flux_rad_segs}}
\end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[scale=0.8, trim=0 100 0 100, clip]{bbrad_kt_corr_lumall_gray.pdf}
\end{center}
\vspace{-4cm}
\caption{BB \& BB+BB emitting areas ($R^2$) vs.~$kT$. The $kT_{\rm Low}$ are shown in red, while $kT_{\rm High}$ are shown in blue.  Large squares and diamonds indicate peak parameters in Segments 1 and 2, respectively.  Gray dashed lines show the correlations reported in Figure\,3 of \citet{lin20} from the analysis of all bursts within the same activation period of the source; the blue shaded regions show the mean values of low and high BB kT in the \citet{lin20} burst sample. The green dashed-dotted line shows the $R^2 \propto L(kT)^{-4}$ line corresponding to a luminosity upper limit of $2.6 \times 10^{41}$\,erg\,s$^{-1}$.  The isoluminosity contours, corresponding to $L=10^{38}-10^{41}$\,erg\,s$^{-1}$, are shown as light green dotted lines.
\label{fig:r2_kt}}
\end{figure*}


\subsection{\mbox{Segment\,2~} Peaks}

We expand below on the intriguing properties of Segment 2. Although the overall lightcurve of the burst forest by itself does not clearly exhibit periodic behavior consistent with the source's spin, $P_{\rm spin} =3.24$\,s, the source spectral properties are correlated with the spin period and the extrapolated pulse profile of the persistent emission observed with NICER at energies below 10 keV.  The high-energy spectral component is clearly present in both COMPT and BB+BB fits  within the burst peaks. We note that all transition times from BB+BB to the single BB thermal spectra correspond to a certain spin phase. In addition, the spectral hardness of the peaks is highest at the minima of the pulse phase, while the spectral parameters remain approximately constant during the trough intervals.  

This strong correlation between the flux and the burst spectra, as well as the phases of soft X-ray emission, provides clues to the emission and viewing geometry for this magnetar.  Because of the need to sustain bursts of duration 0.1\,s for many light-crossing times ($< 0.3$ms), it is widely believed that the bursts emanate from closed field line regions in the magnetosphere.  The extrapolated NICER lightcurve represents the persistent emission from the stellar surface that likely emanates from locations relatively near the magnetic poles \citep{younn20}. This is then more visible to an observer whose instantaneous perspective is more or less over the pole.  As the high X-ray luminosity of the forest of bursts likely indicates highly optically thick conditions (e.g., \citealt{lin11}), their zones act to occult the surface emission at various rotational phases. Thus, it is highly likely that the soft X-ray pulse peak corresponds to a phase where the observer line of sight lies roughly over the poles; this is when the GBM burst spectra are at their softest in both COMPT and BB+BB fits.  Yet, because the $<10\;$keV pulse fraction is modest \citep{younn20}, the quasi-polar surface locales are probably at least partially visible for all burst phases.  The key stellar parameters are the angle $\alpha$ between the magnetic and rotation axes, and the angle $\zeta$ between the observer direction and the spin axis.  The soft X-ray pulse profile is most easily interpreted if $\alpha$ and $\zeta$ are both not too large, likely both less than around $50^{\circ}$; more precise determination of these angles requires detailed surface modeling \citep{Younes-2020-ApJ-J1708} and treatment of magnetospheric occultation, which is beyond the scope of this Letter.

The phase dependence of the burst spectra can then be interpreted, and constraints on the magnetar geometry can be discerned.  Presume first that the various zones of burst activation span a full toroidal closed field line volume, with little or no bias in magnetic azimuth for the overall burst ensemble.  Further, suppose that the hotter regions constitute a toroidal ``belt'' at low altitudes close to where the footpoints of the field lines (that define the toroidal surface) touch the stellar surface. According to Figure\,\ref{fig:r2_kt}, this belt is roughly of effective surface area 3--30 km$^2$, a factor of 10--80 smaller than the cooler majority of the burst activation zone. For the rotational phases during which the polar surface is approximately visible to an observer, the hot zones are also readily visible along with large portions of the cooler outer toroid.  There would also be other rotational phases where the observer's view to the hotter burst locales is occulted by the optically thick cooler outer toroidal zones; these putatively correspond to the low flux phases evident in Figure\,\ref{fig:bbbb_kt_2} where the burst emission is cooler.  Given that the hotter BB component is visible for around 60--70\% of phases, the observer likely does not ever view the magnetar very close to the magnetic equator, when the hotter belt would be occulted by the cooler outer zones for larger portions of the rotation period. This translates to $\zeta$ not being very close to $90^{\circ}$, again probably less than $60^{\circ}$ or so.

A contrasting scenario would be to presume that the hotter BB regions are not adjacent to the surface but rather at higher altitudes, perhaps near the outermost portion of the toroid.  In this circumstance, because Figures~\ref{fig:flux_rad_segs} and \ref{fig:r2_kt} indicate that the areas of the hotter zones for various bursts are on average a factor 20 or so smaller than those for the cooler BB spectral components, the hotter blackbody regions would generally be obscured by the larger opaque toroid for most rotational phases.  This is contrary to the phase-dependent spectral landscape depicted in Figure\,\ref{fig:bbbb_kt_2}. Yet, it is still possible to construct pathological morphology using disparate hotspot locales -- for instance, placing the hotter BB zones on opposite faces (``upper'' and ``lower'') of the toroid and well removed from the magnetic equator.  These would have to be activated in a geometrically non-random manner at preferred magnetic longitudes in order to maintain the presence of hotter BB components for the majority of rotational phases, which would considerably constrain physical mechanisms for burst activation.

This sketch of a geometrical interpretation suggests that locales for the hotter burst zones that are well removed from the surface are somewhat disfavored. Such locales are the preserve of the magnetic reconnection scenario for magnetar burst activation \citep{Lyutikov-2003-MNRAS} that serves as an analog of the picture for solar coronal mass ejections.  They could still be accommodated, however, if the reconnection regions putatively triggering bursts do not sample magnetic longitudes in a uniform manner.  In contrast, a hot near-surface belt for the active toroids is easier to render consistent with the data in Figures~\ref{fig:bbbb_kt_2}, \ref{fig:flux_rad_segs} and \ref{fig:r2_kt}.  This is the expected locale for ignition of bursts by magnetar crustal fractures \citep{TD95,TD96}, likely involving hydromagnetic reconfigurations, and adjacent magnetospheric radiative dissipation. To discern between these two pictures requires more precise modeling of the emission geometry and its phase-dependent spectral signatures. The high quality of this burst storm dataset provides the basis for such a study, eliciting the prospect of discrimination between these burst activation scenarios.

Finally, compared with the COMPT-fit parameters presented in \citet{lin20}, the spectra of the broad peaks (at the end of \mbox{Segment\,1~} and the triple peaks in Segment 2) have harder $E_{\rm peak}$ and index values (Figure\,\ref{fig:compt_param}) than the time-integrated spectra of the bursts within the same active episode, observed before and after the burst forest.  However, when the time-resolved spectra of the peaks are represented with thermal components (i.e., single BB or BB+BB), the luminosity during these peaks does not exceed the upper end of the typical magnetar burst luminosity range ($\lesssim 2.6 \times 10^{41}$\,erg\,s$^{-1}$; see Figure\,\ref{fig:r2_kt}).  A further study to compare the burst forest and the time-resolved spectral properties of bursts within the same activation period is planned.


\acknowledgments
We thank Michael S. Briggs for helpful discussions regarding instrumental deadtime and pile-up issues.
Y.K and \"O.K acknowledge the support from the Scientific and Technological Research Council of Turkey (T\"{U}B\.{I}TAK grant No. 118F344).
M.G.B. acknowledges the generous support of the National Science Foundation through grant AST-1813649.



\clearpage

\begin{longrotatetable}
\begin{deluxetable*}{llcllllll}
\tablecolumns{9}
\tablenum{1}
\tablecaption{Time-resolved spectral analysis results of the Burst-forest Episode.  The spectral time resolution is 4\,ms (binned to S/N $\geq$ 25) for trough spectra and 128\,ms for peak spectra. The uncertainties are for the 1$\sigma$ confidence level.  The spectra without kT$_{\rm high}$ values are better described with a single BB either based on the BIC comparison or because the BB+BB fit failed to constrain the kT$_{\rm high}$ value.  This table is available in its entirety in the machine-readable version. \label{tab:param}}
\tablehead{
\colhead{} & \colhead{} & & \multicolumn{3}{c}{COMPT} & \multicolumn{3}{c}{BB or BB+BB} \\
\cline{4-6}
\cline{7-9}
\colhead{Start Time} & \colhead{End time} & \colhead{Flux $^{a}$} &  \colhead{E$_{\rm peak}$} & \colhead{$\Gamma$} & \colhead{C-Stat/DoF~$^{b}$}  & \colhead{kT$_{\rm low}$} & \colhead{kT$_{\rm high}$} & \colhead{C-Stat/DoF~$^{c}$}\\
\colhead{(s $- T_0^{~\dagger}$)} & \colhead{(s $- T_0^{~\dagger}$)} & \colhead{$( 10^{-6}~\rm{erg~cm}^{-2}~s^{-1})$} & \colhead{(keV)} & \colhead{} & \colhead{} & \colhead{(keV)} & \colhead{(keV)} & \colhead{}
}
\input{table1v2.tex}
\tablecomments{$^\dagger$ $T_0$ = 2020 April 27, 18:26:20 UT (609704785.155 Fermi MET) \\
$^a$ Fluence in $8-200$~keV. \\
$^b$ C-Stat for the COMPT model fit per degree of freedom (DoF). \\
$^c$ C-Stat for the BB+BB model fit or BB fit.\\
$^d$ Non-peak time bin with a low total S/N; used only for constraining the BB $kT$ and the flux.}
\end{deluxetable*}
\end{longrotatetable}



