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

Quarkonia have long been considered as golden probes to study the strongly interacting matter consisting of deconfined quarks and gluons, the quark-gluon plasma (QGP), produced in high-energy heavy-ion collisions~\cite{Busza:2018rrf,Karsch:2000ps,Shuryak:1977ut,Matsui:1986dk,Digal:2001ue}.
Quarkonium states are produced at the early stages of the collision via hard parton scatterings, thus experiencing the full space-time evolution of the medium.
Also, their spectral functions are modified due to color screening~\cite{Matsui:1986dk,Digal:2001ue} and interactions with medium constituents such as gluo-dissociation or Landau damping~\cite{laine:2007,Brambilla:2008cx,Brambilla:2010vq}.
Consequently, the quarkonium yields are expected to be suppressed in heavy ion collisions with respect to expectations from proton-proton (\mbox{$p$+$p$}\xspace) data, following the order of their binding energies.
On the other hand, the yields of quarkonia can be enhanced in the presence of the QGP by recombination processes of uncorrelated as well as correlated quarks~\cite{Gorenstein:2000ck,Andronic:2007bi,Ravagli:2007xx,blaizot:2016jp}.

The modification of the quarkonium yields have been studied by various experiments at RHIC and LHC using the nuclear modification factor quantified as the yield ratio in nucleus-nucleus collisions ($A$+$A$) to that in \mbox{$p$+$p$}\xspace collisions scaled by the average number of binary $NN$ collisions~\cite{STAR:2013kwk,CMS:2017uuv,ATLAS:2018hqe,CMS:2018zza,ALICE:JpsiRAA2019,ALICE:2021UpsForward,ATLAS:2022xso,CMS-PAS-HIN-21-007}.
One of the most remarkable signatures is the ordered suppression of \Ups 1, \Ups 2, and \Ups 3 mesons by their binding energies reported in LHC~\cite{ALICE:2021UpsForward,ATLAS:2022xso,CMS:2018zza,CMS-PAS-HIN-21-007}.

To better understand the in-medium effects of quarkonia in $A$+$A$ collisions in a sophisticated way, it is important to study the ``cold nuclear matter" (CNM) effects which are typically probed using proton-nucleus (\mbox{$p$+$A$}\xspace) collisions. Modification of parton distribution functions in the nucleus~\cite{Vogt:2015uba}, energy loss~\cite{Arleo:2015lja} or nucleus absorption~\cite{McGlinchey:2012bp,Arleo:1999af}, and interactions with comoving particles~\cite{Capella:1996va,Ferreiro:2014bia,Ferreiro:2018wbd} are examples of CNM effects.
On the other hand, various experiments have reported capital results, suggesting a QGP-like behavior of the created medium also in smaller collision systems, such as the observation of long-range collective azimuthal correlations in high multiplicity regions~\cite{CMS:2010ifv,ATLAS:2015hzw,CMS:2016fnw,ALICE:2021nir,CMS:2012qk,ALICE:2012eyl,LHCb:2015coe,CMS:2015yux,CMS:2019wiy,STAR:2015kak,PHENIX:2018lia}.
Therefore, sophisticated phenomenological studies in such interactions become the subject that is sensitive to understanding the quarkonium production in small collision systems.

In this paper, we report a detailed study of the in-medium effects for \Ups 1, \Ups 2, and \Ups 3 mesons in proton-lead (\mbox{$p$+Pb}\xspace), proton-oxygen (\mbox{$p$+O}\xspace), and oxygen-oxygen (\mbox{O+O}\xspace) collisions.
Theoretical calculations for dissociation of \Ups n~\cite{Hong:2019ade} are incorporated with the SONIC framework~\cite{Romatschke:2015gxa} to describe the time evolution of the medium.
The dissociation component is constraint in potential non-relativistic QCD (pNRQCD) limits, and coupled into the Boltzmann equation. The thermal width is calculated based on hard thermal loop (HTL) perturbation theory using the Bethe-Salpeter amplitude.
We report the nuclear modification factors and the second-order Fourier coefficient (\mbox{$v_2$}\xspace) of the azimuthal distribution of \Ups n mesons in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions, and the contribution of feed-down from higher excited states are considered to compare with the experimental data properly.
For the demonstration of the framework, we also present the results in \mbox{Pb+Pb}\xspace collisions and compare them with the experimental results.

\section{Simulation Framework}
\label{sec:Simul}

The simulation framework is composed of two parts, hydrodynamics simulation for background medium and medium response of quarkonia.
For the first part, we follow the procedure described in Ref.~\cite{Lim:2018huo}.
Initial energy density in the transverse plane of in the collision of nucleus-nucleus or proton-nucleus is obtained with the MC-Glauber framework~\cite{Miller:2007ri}.
In the MC-Glauber, a nucleon-nucleon inelastic cross section of 72 mb~\cite{ATLAS:2016tc} is employed for small system collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace= 8 TeV at the LHC.
A Gaussian of width $\sigma=0.4$ fm is used to describe the energy deposition of each nucleon participating in at least one inelastic collision.
The deposited energy distribution from all wounded nucleons in each event is converted into energy density for the hydrodynamic simulation (SONIC).
A single scale factor is used for all events in a certain collision system, and the scale factor is determined to match the charged particle multiplicity ($dN_{ch}/d\eta$) at mid-rapidity in \mbox{$p$+Pb}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8.16 TeV~\cite{CMS:2017shj,ALICE:2018wma}.
The same scale factor is used for \mbox{$p$+O}\xspace and \mbox{O+O}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV by assuming the scale factor does not change much in the collision systems with a similar number of participants.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.9\linewidth]{Figures/SONIC_mult_small.pdf}
    \caption{Charged particle multiplicity distribution in \mbox{$p$+O}\xspace, \mbox{O+O}\xspace, and \mbox{$p$+Pb}\xspace collisions at $\mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace=8$ TeV from SONIC.}
    \label{fig:sonic_mult_small}
\end{figure}

Figure~\ref{fig:sonic_mult_small} shows the charged particle multiplicity distributions of unbiased (0--100\%) \mbox{$p$+O}\xspace, \mbox{O+O}\xspace, and \mbox{$p$+Pb}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV from SONIC.
The mean charged particle multiplicity for \mbox{$p$+O}\xspace and \mbox{O+O}\xspace with the scale factor obtained from \mbox{$p$+Pb}\xspace collisions is 10.0 and 31.2, respectively.
We also check the elliptic flow of charged particle in 0--5\% high multiplicity events as shown in Fig.~\ref{fig:sonic_v2_small}, and the \mbox{$p$+Pb}\xspace result is slightly lower than \mbox{$p$+O}\xspace and \mbox{O+O}\xspace results at higher \mbox{$p_T$}\xspace.
Note that we use a shear viscosity to entropy density ratio of 0.08 and a bulk viscosity ratio of zero in the SONIC calculations.
There are other frameworks for initial conditions~\cite{Moreland:2014oya,Schenke:2012wb}, but we use only the MC-Glauber initial condition as a simple case in this study.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.9\linewidth]{Figures/SONIC_v2_small.pdf}
    \caption{Elliptic flow as a function of \mbox{$p_T$}\xspace for charged particles in 0--5\% central \mbox{$p$+O}\xspace, \mbox{O+O}\xspace, and \mbox{$p$+Pb}\xspace collisions at $\mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace=8$ TeV from SONIC.}
    \label{fig:sonic_v2_small}
\end{figure}

During the SONIC simulation, temperature distributions at multiple time steps are stored for medium response of Upsilons.
At the beginning of the second part, we generate \Ups n mesons, and the $x$ and $y$ positions are determined based on the initial energy density distribution.
We used a Tsallis fit to the \mbox{$p_T$}\xspace distribution of \Ups 1 in $\sqrt{s}=5.02$ TeV~\cite{CMS:2018zza} to sample \mbox{$p_T$}\xspace of Upsilons, and $p_{x}$ and $p_{y}$ are determined based on a randomly assigned azimuthal angle. The initial \mbox{$p_T$}\xspace distribution is assumed to be the same for 5.02 TeV and 8 TeV. Also, any difference in the \mbox{$p_T$}\xspace shape, if any, is expected to be negligible because the results are obtained in fine \mbox{$p_T$}\xspace bins.
The medium response is simulated based on the procedure described in Ref.~\cite{Hong:2019ade}.
The fraction of survived Upsilons for a certain time step ($\Delta t$) is calculated as:

\begin{align}
\label{eq:rate}
    \frac{N(t+\Delta t,p_{T})}{N(t,p_{T})} = e^{- \int^{t+\Delta t}_{t} d t^{\prime} \Gamma_\mathrm{diss}(t^{\prime},p_{T})},
\end{align}
where $\Gamma_\mathrm{diss}$ is the thermal width depending on the medium temperature and Upsilon \mbox{$p_T$}\xspace, and $\Delta t$ is set to 0.02 fm/$c$ in the hydrodynamic simulation.
The medium temperature is from the SONIC simulation according to $x$, $y$, and $t$ values. The temperature in which the quarkonium's in-medium binding energy reaches zero, i.e. dissociation temperature, is set as 600 MeV, 240 MeV, and 190 MeV for \Ups 1, \Ups 2, and \Ups 3, respectively~\cite{MPS:2013qqgp,HS:2006cdqb}.
The thermal width from numerical calculations considering gluo-dissociation and inelastic parton scattering is obtained from Ref.~\cite{Hong:2019ade}.
Figure~\ref{fig:thermal_width} shows thermal width as function of temperature for \Ups 1 (left), \Ups 2 (middle), and \Ups 3 (right), and each line represents a distribution for different \mbox{$p_T$}\xspace.
It generally increases with \mbox{$p_T$}\xspace, and the thermal width for excite states is larger.
Note that thermal width for \Ups 3 is obtained at the limited \mbox{$p_T$}\xspace range, so we estimated values for $\mbox{$p_T$}\xspace>2~\mathrm{GeV}/c$ using the \mbox{$p_T$}\xspace dependence of \Ups 2 as:

\begin{align}
    \Gamma_\mathrm{diss}^{\Upsilon(3S)}&(p_{T}) = \Gamma_\mathrm{diss}^{\Upsilon(3S)}(2~\mathrm{GeV}/c) \frac{\Gamma_\mathrm{diss}^{\Upsilon(2S)}(p_{T})}{\Gamma_\mathrm{diss}^{\Upsilon(2S)}(2~\mathrm{GeV}/c)}.
\end{align}

A different formation time ($\tau_\mathrm{form}$) of each states is used as 0.5, 1.0, and 1.5 for \Ups 1, \Ups 2, and \Ups 3, respectively, based on values in Ref.~\cite{Du:2017qkv}.
We do not consider dissociation at the pre-resonance stage, so the medium response on the $\Upsilon$ states is turned off for $\tau< \gamma  \tau_\mathrm{form}$, where $\gamma$ is the Lorentz factor.
We found that the overall medium response is sensitive to the choice of the formation time and pre-resonance response.
The setting that can well reproduce heavy-ion results is used for small systems consistently.
After incorporating the regeneration effect in heavy-ion collisions, a detailed systematic study will be performed later.

In each time step after the formation time, a survival rate is calculated for each $\Upsilon$ based on Eq.~\ref{eq:rate}, and the $\Upsilon$ is removed when a random number from the uniform distribution within 0--1 is greater than the survival rate.
If the $\Upsilon$ is survived, the position for the next time step is calculated based on the momentum and time step ($\Delta t$).
We repeat the survival rate calculation with a different temperature at the new position until the temperature is lower than a critical temperature of 170 MeV.
For the simulation results presented in later sections, we use 1000 events for MC-Glauber and SONIC.
The number of generated \Ups n for a certain event scales with the number of inelastic collisions from the MC-Glauber.
In this study, we do not consider a contribution from the regeneration effect which is not expected to be significant in small collision systems.
In addition, there could be other nuclear effects such as modification of parton density, initial-state energy loss, and other final-state effects, but we focus on evaluating the dissociation effect.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/ThermalWidth_Y1S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/ThermalWidth_Y2S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/ThermalWidth_Y3S.pdf}
    \caption{Thermal width as a function of temperature for different \mbox{$p_T$}\xspace of \Ups 1 (top), \Ups 2 (middle), and \Ups 3 (bottom).}
    \label{fig:thermal_width}
\end{figure}

\subsection{Nuclear modification factor}
The nuclear modification factors (\mbox{$R_{pA}$}\xspace and \mbox{$R_{AA}$}\xspace) are defined as the ratios of the production cross sections of $\Upsilon$ mesons in \mbox{$p$+$A$}\xspace and $A$+$A$ collisions to the expected cross sections extrapolated from \mbox{$p$+$p$}\xspace collisions. To make an equivalent comparison, we took the ratio of generated to survived number of $\Upsilon$s in the full medium response simulation as the nuclear modification factor.

\subsection{Elliptic flow}
The elliptic flow ($v_2$) is calculated from the azimuthal angle distribution of survived \Ups n with respect to the event plane angle

\begin{align*}
    \frac{dN}{d(\phi - \Psi)} \propto 1 + 2v_{2}(p_T) \cos (2 (\phi - \Psi) ),
\end{align*}
where $\phi$ is the azimuthal angle of the \Ups n, and $\Psi$ is the event plane angle.
The event plane angle for each event is calculated with the initial energy density profile from the MC-Glauber.

\subsection{Feed-Down correction}
\label{sec:FDCorr}
To simulate the medium response of inclusive bottomonium production, the  contributions from feed-down decays, i.e. decays from higher excited states, need to be carefully considered.
Here we consider only strong or electromagnetic decay modes of excited states in the feed-down decays, as the decays from $H$, $Z$, and $W$ bosons to $\Upsilon$ states have negligible effects to the inclusive yields.
In general, the feed-down component from a $Q_{m}$ state to a $Q_{n}$ state ($Q_{m}$ being a higher excited state) can be calculated as

\begin{equation}
\label{eq:FD}
    \mathcal{F}^{Q_{m}}_{Q_{n}} = \mathcal{B}(Q_{m}\rightarrow Q_{n})\frac{\sigma_{Q_{m}}}{\sigma_{Q_{n}}},
\end{equation}
where $\mathcal{B}$ denotes the branching ratio of $Q_{m}$ into $Q_{n}$ and $\sigma_{Q}$ referring the cross section of the corresponding state.
We used publish experimental data in \mbox{$p$+$p$}\xspace collisions at LHC by CMS~\cite{CMS:2015FD} and LHCb~\cite{LHCb:2015FD,LHCb:2015FDchiB} to estimate the feed-down fraction.
Figure~\ref{fig:FeedDown} shows the feed-down fraction from higher excited states for each \Ups 1, \Ups 2, and \Ups 3 states.
The feed-down fractions from CMS are calculated using the presented cross section ratios multiplied by the PDG world-average branching ratios~\cite{Zyla:2020zbs} as in Eq.~\ref{eq:FD}.
The data points of each feed-down component are fitted with an empirical function and their sum drawn as the black solid line which is the total amount of feed-down fraction for the given $\Upsilon$ state.
For \Ups 2 and \Ups 3, the fit function is parameterized to be the same as for the P-wave state feed-down function of \Ups 1 because of the absence of experimental measurements at low-\mbox{$p_T$}\xspace. Afterwards, the function is normalized to match the data points at high-\mbox{$p_T$}\xspace for each \Ups 2 and \Ups 3 meson, respectively.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/FeedDown/FeedDownTo1S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/FeedDown/FeedDownTo2S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/FeedDown/FeedDownTo3S.pdf}
    \caption{Feed-down fraction for \Ups 1 (top), \Ups 2 (middle), and \Ups 3 (bottom) states.}
    \label{fig:FeedDown}
\end{figure}

The feed-down contribution for the nuclear modification factor and elliptic flow of inclusive \Ups n is considered by calculating an weighted average of the quantities for possible states as

\begin{align*}
    R_{n}(p_T) = \sum_{i} R_{i}(p_T)  \mathcal{F}^{Q_{i}}_{Q_{n}} (p_T),
\end{align*}
where $R_{n}$ is the weighted averaged value for a certain \Ups n state, $R_{i}$ is the value for a certain state contributing the the \Ups n state, and $\mathcal{F}^{Q_{i}}_{Q_{n}}$ is the feed-down fraction. The effect of possible \mbox{$p_T$}\xspace shifts by the decays are ignored as the daughter particle carries mostly the mother $\Upsilon$'s momentum because of their similar masses. Since we do not have the nuclear modification factor and elliptic flow for $\chi_b$ states, it is assumed that $R_{\Upsilon(2S)} \approx R_{\chi_{b}(1P)}$ and $R_{\Upsilon(3S)} \approx R_{\chi_{b}(2P)} \approx R_{\chi_{b}(3P)}$ like the study in Ref.~\cite{Hong:2019ade}.

\clearpage

\section{Results and Discussions}
\label{sec:Results_Disscussion}

This section presents and discuss the obtained results of nuclear modification factors (\mbox{$R_{pA}$}\xspace, \mbox{$R_{AA}$}\xspace) and elliptic flow (\mbox{$v_2$}\xspace) for \Ups 1, \Ups 2, and \Ups 3 mesons in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV.
The performance of our framework is tested in nucleus-nucleus collisions and compared with the experimental data at LHC as described in Sec.~\ref{sec:res_aa}.
The results of the nuclear modification factors and \mbox{$v_2$}\xspace in the three small collision systems are shown and discussed in Sec.~\ref{sec:res_rpa} and Sec.~\ref{sec:res_vtwo}, respectively.
For all results, the uncertainty represents the statistical uncertainty that arises from the number of generated events. Also, the feed-down corrections are applied in all calculations as described in Sec.~\ref{sec:FDCorr}.

\subsection{Framework demonstration in \mbox{Pb+Pb}\xspace }
\label{sec:res_aa}
The medium response of the full simulation is demonstrated in \mbox{Pb+Pb}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 5.02 TeV. Figure~\ref{fig:raa_npart_pbpb} shows the \mbox{$R_{AA}$}\xspace curves for \Ups 1, \Ups 2, and \Ups 3 as a function of $\langle \mbox{$N_{\rm part}$}\xspace \rangle$ in \mbox{Pb+Pb}\xspace collision at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 5.02 TeV together with the measurements from CMS~\cite{CMS:2018zza,CMS-PAS-HIN-21-007}.
The calculations agree well with the experimental data for \Ups 1 after applying feed-down corrections, while this correction seems to have restricted effects for \Ups 2 and \Ups 3. The results of the excited $\Upsilon$ states show consistency with data in peripheral collisions, although deviations are found towards central \mbox{Pb+Pb}\xspace collisions.
The discrepancy in central collisions might be present due to the exclusion of recombination processes in this paper, which is expected to be more prominent at large $\langle \mbox{$N_{\rm part}$}\xspace \rangle$ values.
\begin{figure}[htb]
    \centering
        \includegraphics[width=0.9\linewidth]{Figures/RAA_Npart_PbPb_5TeV.pdf}
    \caption{Calculated \mbox{$R_{AA}$}\xspace for $\Upsilon$ states as a function of $\langle \mbox{$N_{\rm part}$}\xspace \rangle$ in \mbox{Pb+Pb}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 5.02 TeV. The uncertainties of the curves represent the statistical uncertainty to the corresponding simulated events. The data points are taken from results by CMS~\cite{CMS:2018zza,CMS-PAS-HIN-21-007}.}
    \label{fig:raa_npart_pbpb}
\end{figure}

Figure~\ref{fig:v2_aa_pt} shows the computed results of \mbox{$v_2$}\xspace as a function of \mbox{$p_T$}\xspace for \Ups 1 mesons under the same conditions for the \mbox{Pb+Pb}\xspace \mbox{$R_{AA}$}\xspace calculations. The results indicate very small \mbox{$v_2$}\xspace values in \mbox{Pb+Pb}\xspace collisions for \Ups 1 mesons which is consistent with the measurements from ALICE~\cite{ALICE:2019Yv2} and CMS~\cite{CMS:2021Yv2}.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.9\linewidth]{Figures/v2/v2_Pt_PbPb.pdf}
    \caption{Calculated \mbox{$v_2$}\xspace for \Ups 1 as a function of \mbox{$p_T$}\xspace in \mbox{Pb+Pb}\xspace collsions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 5.02 TeV. The uncertainties of the curves represent the statistical uncertainty to the corresponding simulated events. The data points are taken from results by ALICE~\cite{ALICE:2019Yv2} and CMS~\cite{CMS:2021Yv2}.}
    \label{fig:v2_aa_pt}
\end{figure}

\subsection{Nuclear modification factor}
\label{sec:res_rpa}

The calculated nuclear modification factor as a function of $dN_{ch}/d\eta$ for \Ups 1, \Ups 2, and \Ups 3 mesons in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV are shown in Fig.~\ref{fig:rpa_mult}.
The nuclear modification factor shows a gradual decrease with increasing event multiplicity for all three $\Upsilon$ states in all three collision systems. Also, the amount of suppression towards higher multiplicity events is found to be sequentially ordered following the magnitude of the corresponding $\Upsilon$ binding energy.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_mult_pPb.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_mult_pO.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_mult_OO.pdf}
    \caption{Nuclear modification factors for \Ups 1, \Ups 2, and \Ups 3 as a function of $dN_{ch}/d\eta$ in \mbox{$p$+Pb}\xspace (top), \mbox{$p$+O}\xspace (middle), and \mbox{O+O}\xspace (bottom) collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV. The uncertainties represent the statistical uncertainty to the corresponding simulated events.}
    \label{fig:rpa_mult}
\end{figure}

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_mult_1S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_mult_2S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_mult_3S.pdf}
    \caption{Calculated nuclear modification factor as a function of $dN_{ch}/d\eta$ in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV each for \Ups 1 (top), \Ups 2 (middle), \Ups 3 (bottom). The uncertainties represent the statistical uncertainty to the corresponding simulated events.}
    \label{fig:rpa_mult_coll}
\end{figure}

To diagnose the strength of modification among different collision systems, the nuclear modification factors are separately displayed for each $\Upsilon$ state as shown in Fig.~\ref{fig:rpa_mult_coll}.
In the low multiplicity region ($dN_{ch}/d\eta < 25)$, the suppression level of $\Upsilon$ states is found to be similar in \mbox{$p$+Pb}\xspace and \mbox{$p$+O}\xspace collisions, while a slightly less suppression is seen in \mbox{O+O}\xspace collisions in particular for \Ups 2 and \Ups 3.
Despite the medium size in \mbox{O+O}\xspace collisions being the largest and that in \mbox{$p$+O}\xspace collisions the smallest, the amount of suppression is not proportional to the size of the created system.
The fact that the system size is larger in a given multiplicity event rather implies a smaller energy density of the medium. Therefore, we conclude that the weaker suppression in \mbox{O+O}\xspace collisions is accounted for the smaller energy density compared to \mbox{$p$+O}\xspace and \mbox{$p$+Pb}\xspace collisions.
At higher multiplicity region, where the results in \mbox{$p$+Pb}\xspace and \mbox{O+O}\xspace collisions are compared, $\Upsilon$ states are stronger suppressed in \mbox{$p$+Pb}\xspace collisions than in \mbox{O+O}\xspace collisions which is consistent with the finding at low multiplicity because of the same reason.

Figure~\ref{fig:rpa_pt} shows the same quantities presented as a function of \mbox{$p_T$}\xspace. The ordering of the nuclear modification factors is similar at low-\mbox{$p_T$}\xspace as in Fig.~\ref{fig:rpa_mult_coll}, but found to be reversed between \mbox{$p$+Pb}\xspace and \mbox{O+O}\xspace collisions at high-\mbox{$p_T$}\xspace.
Since the formation time of the $\Upsilon$ states is delayed towards the higher \mbox{$p_T$}\xspace region, the effective interaction time with the medium is also reduced. Such effect is more prominent in \mbox{$p$+Pb}\xspace collisions because of the smaller initial medium size.
The relation between the medium size and energy density affects the $\Upsilon$ suppression in the opposite direction for low- and high-\mbox{$p_T$}\xspace.
These results show that the shape of the \mbox{$p_T$}\xspace-dependent nuclear modification factor is sensitive to the formation time and initial collision geometry.
\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_Pt_1S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_Pt_2S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_Pt_3S.pdf}
    \caption{Nuclear modification factors for \Ups 1, \Ups 2, and \Ups 3 as a function of \mbox{$p_T$}\xspace in \mbox{$p$+O}\xspace, \mbox{$p$+Pb}\xspace, and \mbox{O+O}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV. The uncertainties represent the statistical uncertainty to the corresponding simulated events.}
    \label{fig:rpa_pt}
\end{figure}

We also compare our calculations with measured data by CMS in \mbox{$p$+Pb}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 5.02 TeV~\cite{CMS:2022wfi} as shown in Fig.~\ref{fig:rpa_pt_cms}.
The results are in good agreement with data for \Ups 1, whereas deviations are seen for \Ups 2 and \Ups 3.
Note that our \mbox{$R_{pA}$}\xspace values are for 8 TeV, in which the average multiplicity is about 15\% higher than that at 5.02 TeV~\cite{CMS:2017shj}. This difference affects the modification of $\Upsilon$ yields more significantly at low-\mbox{$p_T$}\xspace, where they mostly experience the full medium evolution.
Although the mean \mbox{$p_T$}\xspace of the highest \mbox{$p_T$}\xspace bin of the CMS measurements for the \Ups 2 and \Ups 3 meson is expected to be slightly lower than the center of the bin, our calculations still overshoot the experimental data.
As discussed in Fig.~\ref{fig:rpa_pt}, this disagreement is possibly related to the formation time and initial collision geometry.
Note that initial-state model calculations considering modification of nuclear parton distribution functions and energy loss show about 10\% suppression of \Ups 1 at mid-rapidity~\cite{CMS:2017shj,Vogt:2015uba,Arleo:2014oha}.

\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/RpPb/RpPb_Pt_pPb_CMS.pdf}
    \caption{Calculated nuclear modification factor for \Ups 1, \Ups 2, and \Ups 3 as a function of \mbox{$p_T$}\xspace in \mbox{$p$+Pb}\xspace  collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV compared with the CMS results at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 5.02 TeV.}
    \label{fig:rpa_pt_cms}
\end{figure}

\subsection{Elliptic flow}
\label{sec:res_vtwo}

The \mbox{$v_2$}\xspace is calculated in the same setting as in the calculation for the nuclear modification factors in Sec.~\ref{sec:res_rpa} and shown in Fig.~\ref{fig:v2_pt}. For all three $\Upsilon$ states, the \mbox{$v_2$}\xspace values are consistent with zero in the overall \mbox{$p_T$}\xspace region in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions.
Also, their \mbox{$v_2$}\xspace values are very similar among the three collision systems. It may indicate that the elongated formation time of $\Upsilon$s towards high-\mbox{$p_T$}\xspace is not as significant as implemented in the employed model.
\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/v2/V2_Pt_Y1S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/v2/V2_Pt_Y2S.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/v2/V2_Pt_Y3S.pdf}
    \caption{Calculated $v_2$ in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions as a function of \mbox{$p_T$}\xspace for \Ups 1 (top), \Ups 2 (middle), \Ups 3 (bottom) at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV. The uncertainties represent the statistical uncertainty to the corresponding simulated events.}
    \label{fig:v2_pt}
\end{figure}
To test the relation between \mbox{$v_2$}\xspace and the amount of suppression, we computed the \mbox{$v_2$}\xspace in the 0--5\% selected high-multiplicity events.
For all $\Upsilon$ states, the \mbox{$v_2$}\xspace values are found to be larger at high-multiplicity only with a small magnitude.
\begin{figure}[htb]
    \centering
        \includegraphics[width=0.8\linewidth]{Figures/v2/V2_Pt_Y1S_HM.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/v2/V2_Pt_Y2S_HM.pdf}
        \includegraphics[width=0.8\linewidth]{Figures/v2/V2_Pt_Y3S_HM.pdf}
    \caption{Calculated $v_2$ in \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions as a function of \mbox{$p_T$}\xspace for \Ups 1 (top), \Ups 2 (middle), \Ups 3 (bottom) at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV in 0--5\% high multiplicity events. The uncertainties represent the statistical uncertainty to the corresponding simulated events}
    \label{fig:v2_pt_hm}
\end{figure}

Although the yields of $\Upsilon$ mesons are more strongly suppressed for excited states as well as for higher multiplicity events, their \mbox{$v_2$}\xspace values show an overall consistency with each other in all studied cases.
In particular, since the low-\mbox{$p_T$}\xspace $\Upsilon$ mesons are traversing very slowly, it is unlikely expected that they escape the medium before it reaches the chemical freeze-out temperature. Therefore, such $\Upsilon$ mesons are not able to capture the anisotropy of the initial collision geometry despite their significant suppression. It is understood that the path-length dependent suppression could lead to a non-zero \mbox{$v_2$}\xspace which can be probed at high-\mbox{$p_T$}\xspace for $\Upsilon$ mesons. However, with our current simulation settings, we did not observe any firm non-zero \mbox{$v_2$}\xspace that deviates from that at low-\mbox{$p_T$}\xspace in small collision systems for all $\Upsilon$ states.

\pagebreak

\section{Summary}
\label{sec:Sum}

We have performed a Monte Carlo simulation study for the medium response of bottomonia in \mbox{$p$+$A$}\xspace and $A$+$A$ collisions. The simulation framework has been developed based on the theoretical calculation of the thermal width of \Ups n~\cite{Hong:2019ade} and the publicly available codes to describe the initial condition and evolution of heavy-ion collisions.
In this initial work, only dissociation effect is considered, and we also considered the contribution of feed-down from higher excited states. To demonstrate the framework, we calculate the nuclear modification factor and elliptic flow of \Ups n and compare to experimental results.
In the nuclear modification factor as a function of the number of participants, the model and data are comparable for \Ups 1 at the entire region, whereas the suppression for \Ups 2 and \Ups 3 are stronger in the model.
In addition to the formation time, which could directly affect the magnitude of nuclear modification, the comparison for excited states would be improved when incorporating the regeneration effect.
In case of the elliptic flow, the model expect $v_2<0.01$ which agrees with the experimental results~\cite{ALICE:2019Yv2,CMS:2021Yv2}.

We extend the framework to small systems, \mbox{$p$+Pb}\xspace, \mbox{$p$+O}\xspace, and \mbox{O+O}\xspace collisions at \mbox{$\sqrt{s_{_\mathrm{NN}}}$}\xspace = 8 TeV. Generally, a stronger modification is observed for higher states ($R_{pA, AA}^{\Upsilon(1S)} > R_{pA, AA}^{\Upsilon(2S)} > R_{pA, AA}^{\Upsilon(3S)}$) even in small system. In low multiplicity events where the system size is very small, less suppression of \Ups 3 is observed due to the late formation time. In comparing different systems at the same multiplicity, a similar suppression is seen in $dN_{ch}/d\eta<25$. At higher multiplicity, the modification in \mbox{O+O}\xspace is weaker than that in \mbox{$p$+Pb}\xspace because the energy density (temperature) in \mbox{O+O}\xspace is lower due to the larger system size. Regarding the elliptic flow results, we obtain very small ($<0.01$) elliptic flow for all three systems, even in high multiplicity events. It will be very interesting to compare with experimental results from the upcoming LHC run with the oxygen ion. It can provide valuable information on sources of nuclear effects on bottomonia production in small systems.

\section{Acknowledgments}

The work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) under Project No. 2018R1A5A1025563 and 2020R1C1C1004985. We also acknowledge technical support from KIAF administrators at KISTI.

\clearpage

\newpage


\section{Introduction}

The usefulness of feature engineering for applications of deep learning has long been considered a settled question in the negative, as neural networks are on their own universal function approximators \cite{GoodBengCour16}. However, even in the context of deep learning, tabular features are often treated with some form of encoding for preprocessing. Automunge \cite{anonymous_github} is a platform for encoding dataframes developed by the authors. This python library was originally built for a simple use case of basic encoding conventions for numeric and categoric features, like z-score normalization and one-hot encodings. Along the iterative development journey we began to flesh out a full library of encoding options, including a series of options for numeric and categoric features that now include scenarios for normalization, binarization, hashing, and missing data infill under automation. Although it was expected that these range of encoding options would be superfluous for deep learning, that does not rule out their utility in other paradigms which could range from simple regression, support vector machines, decisions trees, or as will be the focus of this paper, gradient boosting.

The purpose of this work is to present the results of a benchmarking study between alternate encoding strategies for numeric and categoric features for gradient boosted tabular learning. We were particularly interested in validating the library’s default encoding strategies, and found that in both primary performance metrics of tuning duration time and model performance the current defaults under automation of categoric binarization and numeric z-score normalization demonstrated merit to serve as default encodings for the Automunge library. We also found that in addition to our default binarization, even a frequency sorted variant of ordinal encoding on average outperformed one hot encoding.

\section{Gradient Boosting}

Gradient boosting \cite{friedman2000greedy} refers to a paradigm of decision tree learning \cite{quinlan:induction} similar to random forests \cite{10.1023/A:1010933404324} but in which the optimization is boosted by recursively training an iteration’s model objective to correct the performance of the preceding iteration’s model. It is commonly implemented in practice by the XGBoost library \cite{Chen_2016} for GPU acceleration, although there are architecture variations available for different fortes, like LightGPM \cite{NIPS2017_6449f44a} which may train faster on CPU’s than XGBoost (with a possible performance tradeoff).

Gradient boosting has traditionally been found as a winning solution for tabular modality competitions on the Kaggle platform, and its competitive efficacy has even been demonstrated for more sophisticated applications like time series sequential learning when used for window based regression \cite{https://doi.org/10.48550/arxiv.2101.02118}. Recent tabular benchmarking papers have found that gradient boosting may still mostly outperform sophisticated neural architectures like transformers \cite{gorishniy2021revisiting}, although even a vanilla multi layer perceptron neural network could have capacity to outperform gradient boosting with comprehensively tuned regularizers \cite{kadra2021welltuned}. Gradient boosting can also be expected to have higher latency inference than neural networks \cite{https://doi.org/10.48550/arxiv.2110.01889}.

Conventional wisdom is that one can expect gradient boosting models to have capacity for better performance than random forests for tabular applications but with a tradeoff of increased probability of overfitting without hyperparameter tuning \cite{Howard:20}. With both more sensitivity to tuning parameters and a much higher number of parameters in play than random forest, gradient boosting usually requires more sophistication than a simple grid or random search for tuning. One compromise method available is for a sequential grid search through different subsets of parameters \cite{procedural_tuning}, although more automated and even parallelized methods are available by way of black box optimization libraries like Optuna \cite{https://doi.org/10.48550/arxiv.1907.10902}. There will likely be more improvements to come both in libraries and tuning conventions, this is an active channel of industry research.

\section{Feature Encodings}

Feature encoding refers to feature set transformations that serve to prepare the data for machine learning. Common forms of feature encoding preparations include normalizations for numeric sets and one hot encodings for categoric, although some learning libraries may accept categoric features in string representations for internal encodings. Before the advent of deep learning, it was common to supplement features with alternate representations of extracted information or to combine features in some fashion. Such practices of feature engineering are sometimes still applied in gradient boosted learning, and it was one of the purposes of these benchmarks to evaluate benefits of the practice in comparison to directly training on the data.

An important distinction of feature encodings can be considered as those that can be applied independent of an esoteric domain profile verses those that rely on external structure. An example could be the difference between supplementing a feature with bins derived based on the distribution of populated numeric values verses extracting bins based on an external database lookup. In the case of Automunge, the internal library of encodings follows almost exclusively the former, that is most encodings are based on inherent numeric or string properties and do not consider adjacent properties that could be inferred based on relevant application domains. (An exception is made for date-time formatted features which under automation automatically extract bins for weekdays, business hours, holidays, and redundantly encodes entries based on cyclic periods of different time scales \cite{time_encoding}.) The library includes a simple template for integrating custom univariate transformations \cite{anonymous_customtransforms} if a user would like to integrate into a pipeline alternate conventions.

\subsection{Numeric}

Numeric normalizations in practice are most commonly applied similar to our default of z-score `nmbr' (subtract mean and divide by standard deviation) or min-max scaling `mnmx' (converting to range between 0–1). Other variations that may be found in practice include mean scaling `mean'(subtract mean and divide by min max delta), and max scaling `mxab' (divide by feature set absolute max). More sophisticated conventions may convert a distribution shape in addition to the scale, such as the box-cox power law transformation `bxcx' \cite{https://doi.org/10.1111/j.2517-6161.1964.tb00553.x} or Scikit-Learn’s \cite{JMLR:v12:pedregosa11a} quantile transformer `qttf', which both may serve the purpose of converting a feature set to closer resemble a Gaussian distribution. In general, numeric normalizations are more commonly applied for learning paradigms other than those based on decision trees, where for example in neural networks they serve the purpose of normalizing gradient updates across features. We did find that the type of normalizations applied to numeric features appeared to impact performance, and we will present these findings below.

\subsection{Categoric}

Categoric encodings are most commonly derived in practice as a one hot encoding, where each unique entry in a received feature is translated to boolean integer activations in a dedicated column among a returned set thereof. The practice of one hot encoding has shortcomings in the high cardinality case (where a categoric feature has an excessive number of unique entries), which in the context of gradient boosting may be particularly impactful as an inflated column count impairs latency performance of a training operation — or when the feature is targeted as a classification label may even cause training to exceed memory overhead constraints. The Automunge library attempts to circumvent this high cardinality edge case in two fashions, first by defaulting to a binarization encoding instead of one hot, and second by distinguishing highest cardinality sets for a hashed encoding \cite{NIPS1988_82161242} \cite{10.1145/1553374.1553516} \cite{anonymous_customtransforms} which may stochastically consolidate multiple unique entries into a shared ordinal representation for a reduced number of unique entries.

The library default of categoric binarization `1010' refers to translating each unique entry in a received feature to a unique set of zero, one, or more boolean integer activations in a returned set of boolean integer columns. Where one hot encoding may return a set of n columns for n unique entries, binarization will instead return a smaller count of log2(n) rounded up to nearest integer. We have previously seen the practice discussed in the blogging literature, such as \cite{onehot_blog}, although without validation as offered herein.

A third common variation on categoric representations includes ordinal encodings, which simply refers to returning a single column encoding of a feature with a distinct integer representation for each unique entry. Variations on ordinal encodings in the library may sort the integer representations by frequency of the unique entry `ord3' or based on alphabetic sorting `ordl'.

Another convention for categoric sets unique to the Automunge library we refer to as parsed categoric encodings `or19' \cite{anonymous_Strings}. Parsed encodings search through tiers of string character subsets of unique entries to identify shared grammatical structure for supplementing encodings with structure derived from a training set basis. Parsed encodings are supplemented with extracted numeric portions of unique entries for additional information retention in the form received by training.

\section{Benchmarking}

The benchmarking sought to evaluate a range of numeric and categoric encoding scenarios by way of two key performance metrics, training time and model performance. Training was performed over the course of ~1.5 weeks on a Lambda workstation with AMD 3970X processor, 128Gb RAM, and two Nvidia 3080 GPUs. Training was performed by way of XGBoost tuned by Optuna with 5-fold fast cross-validation \cite{NIPS2013_f33ba15e} and early stopping criteria of 50 tuning iterations without improvement. Performance was evaluated against a partitioned 25\% validation set based on a f1 score performance metric, which we understand is a good default for balanced evaluation of bias and variance performance of classification tasks \cite{stevens2020learning}. This loop was repeated and averaged across 5 iterations and then repeated and averaged across 31 tabular classification data sets sourced from the OpenML benchmarking repository \cite{Vanschoren_2014}. Rephrasing for clarity, the reported metrics are averages of 5 repetitions of 31 data sets for each encoding type as applied to all numeric or categoric features for training. The distribution bands shown in the figures are across the five repetitions. The data sets were selected for diverse tabular classification applications with in-memory scale training data and tractable label cardinality.

We found that these benchmarks gave us comfort in the Automunge library's defaults of numeric z-score normalization and categoric binarization. An interesting result was the outperformance of categoric binarization in comparison to one-hot encoding, as the latter is commonly used in mainstream practice as a default. Further dialogue for the interpretation of the results presented in Figures \ref{fig1:test} and \ref{fig2:test} are provided in the Appendix.

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.8\linewidth]{numeric_time_031422.png}
  \caption{Numeric tuning time comparison}
  \label{fig1:sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.8\linewidth]{numeric_performance_031422.png}
  \caption{Numeric model performance comparison}
  \label{fig1:sub2}
\end{subfigure}
\caption{Numeric Results}
\label{fig1:test}
\end{figure}

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.8\linewidth]{categoric_time_031422.png}
  \caption{Categoric tuning time comparison}
  \label{fig2:sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.8\linewidth]{categoric_performance_031422.png}
  \caption{Categoric model performance comparison}
  \label{fig2:sub2}
\end{subfigure}
\caption{Categoric Results}
\label{fig2:test}
\end{figure}

\section{Conclusion}

We hope that these benchmarks may have provided some level of user comfort by validating the default encodings applied under automation by the Automunge library of z-score normalization and categoric binarization, both from a training time and model performance standpoint. If you would like to try out the library we recommend the tutorials folder found on GitHub \cite{anonymous_github} as a starting point.

\newpage