§ INTRODUCTION

In the study of high-energy particle physics,  general-purpose event generator stands as a powerful tool, providing a platform to simulate the intricate processes of hadron formation and decay [1, 2, 3, 4]. One of the natural steps in the simulation of high-energy collisions in  is the emergence of observable hadronic states through fragmentation and decay of quarks and gluons at the end of a sequential parton shower via an overall non-perturbative cluster hadronization model [5]. In this context, it would be productive to transmit the polarization information of the heavy quarks appearing at the end of the parton shower to the hadronization handler to accurately account for the net polarization observed in excited mesons and heavy baryons. This information ensures that the interplay between quark spins and flavours is appropriately considered during the hadronization process, leading to a more realistic depiction of the polarization properties of the resulting particles. This is possible by implementing the effects of the Heavy Quark Effective Theory (HQET) [6, 7, 8, 9].

HQET is a framework that provides a systematic way to describe the behaviour and interactions of heavy quarks within hadrons. In this context, a quark is considered “heavy" if its mass is much larger than the QCD energy scale, $\Lambda_{\rm QCD}$. This heavy quark is then treated as moving at a sub-relativistic velocity. HQET takes advantage of this mass hierarchy and velocity separation to simplify the mathematical description of heavy quark interactions within hadrons, assuming that the velocity and the spin of the heavy quarks can be treated as independent variables, resulting in a decoupling between the dynamics of the heavy quark and the lighter degrees of freedom. This allows for a separation of scales, where the interactions involving the heavy quark can be analyzed using a perturbative expansion while the interactions involving the lighter degrees of freedom are treated non-perturbatively. The emergence of this heavy quark spin-flavour symmetry and the associated effects allow for model-independent predictions of certain properties of heavy hadrons, such as their spectra and decay rates.

The above endeavour involves modelling the decays of excited heavy hadrons, specifically those containing a single charm or bottom quark. In the framework of heavy quark symmetry, these decays are not only determined by the possibilities of available decay modes but also exhibit inherent relationships between the associated couplings. The limited experimental data concerning many of these decays necessitates a reliance on these symmetries to determine the decay modes, widths and branching ratios. Therefore, it would be pertinent to refine the existing prescriptions for the electromagnetic and strong-isospin-violating decay modes of excited heavy mesons while employing the effects of HQET. The latter modes become more significant in cases where strong isospin-conserving decays are either disallowed or significantly kinematically suppressed.

Undoubtedly, incorporating the heavy quark spin-flavour symmetry and consequent enhancements to the heavy excited hadronic states, particularly the isospin-violating decay modes, will trigger a cascade of effects that can reverberate through 's hadronization and decay handlers, thereby influencing the production rates of all types of hadronic states. Furthermore, the introduction of new model-independent parameters into 's cluster hadronization model, with the aim of refining the kinematic mass threshold that governs cluster splittings holds the potential to ripple across both baryon and meson production rates. This prompts the necessity of undertaking a general tuning for hadronization parameters, alongside the pertinent parton shower parameters.

In this study, we outline the advancements made in predictive modelling of heavy hadron production, focusing on enhancing  event generator. Section <ref> provides a thorough examination of HQET and Spin-Flavour Symmetry, breaking it down further into the kinematics of fragmentation, and the polarization mechanisms for excited heavy mesons and baryons. Section <ref> discusses the decay processes of these excited heavy mesons, laying the groundwork for Section <ref>, which scrutinizes relevant observables. In Section <ref>, a comprehensive tuning methodology for  is introduced, driven by a brute-force, multilayered approach to parameter optimization. Section <ref> then evaluates the efficacy of these modifications by analyzing heavy hadron production rates using an extensive dataset. A brief summary and conclusions are presented in Section <ref>. The paper also includes an appendix detailing the signature of heavy hadrons in the $\Upsilon(4S)$ continuum region, providing further context for the validation of our study.

The calculations pertaining to induced polarizations of excited heavy mesons and heavy baryons closely follow the methodology of [10], while the enhancements to decay modes build upon the notation of [11]. The results showcased herein are generated using , incorporating internally generated matrix elements and 's  parton shower [12, 13]. The subsequent analysis is conducted employing  (v3.1.8) [14]. Our tuning endeavours draw inspiration from the methodology outlined in  (v2.3.7) [15], albeit substantially adapted to accommodate the tuning of $12$ independent parameters. The outcomes of this comprehensive exploration, embedded within the  general-purpose event generator, will be made publicly accessible alongside the anticipated release of .

§ HEAVY QUARK EFFECTIVE THEORY AND SPIN-FLAVOUR SYMMETRY

A hadronic state can be considered as a non-recoiling source of colour, if its constituents can be divided into separate heavy ($Q$) and light ($q$) degrees of freedom, while $m_Q \gg \Lambda_{\rm QCD}$ [16, 17, 18, 19, 20]. In such a scenario, the effects associated with the heavy quark colour magnetic moment decouple from hadronic properties and a heavy quark spin-flavour symmetry emerges [10]. This may be used in producing model-independent predictions of heavy hadron spectra, weak matrix elements and strong decay rates. A particularly interesting side-effect of having this spin-flavour symmetry is the ability to distinguish between short- and long-distance interactions in a given hadronization sequence, say at the end of a procedurally generated parton shower.

Generally speaking, the strongly interacting particles produced during the hard processes and the subsequent parton showers in a simulated event can be considered to have relatively large momenta. Henceforth, the corresponding perturbative interactions during these stages are expected to take place in short time scales. Meanwhile, the non-perturbative fragmentation processes that wrap these individual branches are inherently long-distance and therefore will take place over longer time scales. On the other hand, the ensuing hadronization processes, regardless of their model dependencies, are expected to take place at length scales of the order of $\mathcal{O}(\Lambda_{\rm QCD}^{-1})$. This means that the resulting redistribution of energies would be much smaller than $m_Q$ and it would be completely plausible to assume that the velocity, mass and spin of the heavy quark state (which are determined by the calculable short-distance physics) remain unchanged, suggesting a decoupling between the dynamics of the heavy and light quark stats [21, 22, 23].

Interestingly, the above argument can be extended to the production of excited heavy mesons and heavy baryons in a given hadronization event, wherever it is possible to tag the produced hadrons with the spin polarization information of their heavy quark ingredients [10]. Such a provision usually comes with a couple of strings attached; (i) If the phase-space of the light-quark state allows its angular momentum to become of the order of $\mathcal{O}(m_Q^{-1})$, the process can redistribute this angular momentum, meaning that the out-going heavy quark polarization will depend on the polarization of the light degree of freedom created in the fragmentation process. (ii) While conserving parity, fragmentation can in principle produce anisotropic light degrees of freedom along its axis. The alignment of a light degree of freedom with spin $j$ may be characterized in the form of a model-dependent dimensionless parameter, $\omega_j$ [10].

§.§ Kinematics of Fragmentation

To properly translate the implications of HQET and the above-mentioned spin-flavour symmetry, one has to start by projecting the kinematics and the time scales of heavy quark fragmentation. Here we closely follow the footsteps of [10], by assuming to work in the rest frame of the heavy quark while identifying the preferred direction as the momentum of the parton shower progenitor. In this setup, the only interaction that can manipulate the spin of the heavy quark, i.e. the colour magnetic moment of the light degrees of freedom, is sufficiently suppressed by a factor of $m_Q^{-1}$. Collectively, the hadron can possess a total spin of $j_{\pm} = j_q \pm 1/2$, with $j_q$ being the spin of the light state. We can, therefore, identify the $j_{\pm}$ configurations with a spin multiplet ($H$, $H^\star$), and assume a small mass splitting $\Delta m = m_{H^\star} - m_{H}$.

By virtue of having a heavy quark constituent, one can consider the $H$ and $H^\star$ states to be unstable. Hence, they can decay weakly to lighter hadrons, where it is safe to assume an identical decay rate

\begin{equation}
\Gamma(H \to X) = \Gamma(H^\star \to X).
\end{equation}

Alternatively, $H^\star$ can strongly or radiately decay to $H$, with rate

\begin{equation}
\gamma(H^\star \to H X) \propto \Phi_{\rm phase-space}^{H^\star \to H X} \times |\mathcal{M}(H^\star \to H X)|^2 \sim \mathcal{O}(m_q^{-(2+n)}),
\label{gamma}
\end{equation}

where $n \geq 1$, suggesting $\Delta m \gg \gamma$. Meanwhile, the relation between $\Gamma$, $\Delta m$ and $\gamma$ would determine the nature of the fragmentation. Particularly, the case of $\Gamma \gg \Delta m \gg \gamma$ would allow for rapid decay of heavy hadrons while decoupling the colour magnetic moments of the heavy and the light degrees of freedom. This would be valid for both strong or weak decays and even for $\Gamma \sim \Lambda_{\rm QCD}$, where the heavy quark can go through partial hadronization before it decays. Both other possible cases, namely $\Delta m \gg \Gamma \gg \gamma$ and $\Delta m \gg \gamma \gg \Gamma$ result in depolarization of the heavy quark from its initial orientation. The above argument suggests that under most conditions, the angular distribution of decay products gives no information on the polarization of the heavy quark, unless the condition $\Gamma \gg \Delta m \gg \gamma$ is realized. This is the so-called Falk-Peskin “no-win” theorem, suggesting that the polarization of heavy quarks can remain intact only though strong, weak or radiative decays of heavy excited mesons and heavy baryons [10].

§.§ Polarization of Excited Heavy Mesons

To explore the ramifications of HQET and spin-flavour symmetry regarding the polarizations of heavy excited mesons, we initially focus on the charm sector, with particular attention directed at the observed excited charmed mesons $D_0^{\star}$, $D_1$, $D_1'$ and $D_2^{\star}$. Note that this discussion can be readily extended to the case of bottom mesons. Within this framework, the spins of the light and heavy degrees of freedom in the first and second excited $D$ mesons can be comprehended as follows: The first excited $D$ mesons are a result of coupling a spin-${1 \over 2}$ light quark with a spin-${1 \over 2}$ heavy quark, thereby generating a composite angular momentum of $j = 1$. This leads to the potential quantum numbers $j^P = 0^+$ or $j^P = 1^+$, respectively corresponding to $D_0^{\star}$ and $D_1$[At order $1/m_c$, there is mixing between the $D_1$ and the $D_1'$ states as they possess identical quantum numbers.] states. In contrast, the second excited $D$ meson engages higher angular momentum states for both light and heavy components, resulting in quantum numbers $j^P = (3 / 2)^+$ for the second excited $D$ meson. This corresponds to the $D_2^{\star}$ charmed state.

Assuming that the initial $c$ quark entering the fragmentation processes has left-handed spin polarization, $j_Q = -1/2$, a heavy charmed meson can be formed by combining with light degrees of freedom with $j_q = {1/2}$. One would expect that the colour magnetic interaction becomes decoupled form in such a process, due to its time scale, leaving the spin-orientation of the light degrees independent of the charm quark and distributed uniformly, i.e. $j_q^{(3)} = \pm 1 / 2$ with equal probabilities:

\begin{equation}
{\left|\downarrow\,\rangle \right.}_c \; {\left|\downarrow\,\rangle \right.}_q
\quad {\rm and } \quad
{\left|\downarrow\,\rangle \right.}_c \; {\left| \uparrow\,\rangle \right.}_q
\nonumber
\end{equation}

This scenario becomes more intriguing when considering the $D_1$ and $D_2^{\star}$ meson states. In conjunction with the assumed left-handed polarization of the charm quark, we now have light degrees of freedom with $j_q = 3/2$, which can manifest in any of the four possible helicity states. Parity invariance dictates that the probability of forming a specific helicity state cannot rely on the sign of its helicity, $j_3$, although probabilities may differ for states with distinct helicity magnitudes, $|j_q^{(3)}|$. Defining the parameter $\omega_j$ ($0 \leq \omega_j \leq 1$), that is the likelihood of fragmentation leading to a state with the maximum value of $|j_q^{(3)}|$ in a system with light degrees of freedom with spin $j_q$, one may derive the probabilities for the emergence of different helicity states of the light degrees as

$j_q^{(3)}$     -3/2     -1/2     1/2     3/2

Probability     ${1 \over 2} \omega_{3 \over 2}$
${1 \over 2} (1-\omega_{3 \over 2})$
${1 \over 2} (1-\omega_{3 \over 2})$
$\begin{matrix}  \\ \\ \end{matrix} {1 \over 2} \omega_{3 \over 2} \begin{matrix}  \\ \\ \end{matrix}$

Probability distribution of the helicity of the light degrees in excited heavy mesons.

The combination of the left-handed $c$ spin state with a specific light degrees of freedom helicity $j_q^{(3)}$ results in a coherent linear superposition of the charmed states with helicity $j = j_q + j_Q$. This rationale leads to the following populated helicity states and their corresponding probabilities for $D$, $D_0^{\star}$, $D_1$ and $D_2^{\star}$ [10]:

$j^{(3)}$     $-2$     $-1$     $0$     $+1$      $\begin{matrix}  \\ \\ \end{matrix} +2 \begin{matrix}  \\ \\ \end{matrix}$

$D$     $-$     $-$     ${1 \over 4}$     $-$     $\begin{matrix}  \\ \\ \end{matrix} - \begin{matrix}  \\ \\ \end{matrix}$

$D^{\star}$     $-$
${1\over 2}$
${1\over 4}$
$\begin{matrix}  \\ \\ \end{matrix} - \begin{matrix}  \\ \\ \end{matrix}$

$D_1$     $-$
${1 \over 8} (1-\omega_{3 \over 2})$
${1 \over 4} (1-\omega_{3 \over 2})$
${3 \over 8}(1-\omega_{3 \over 2})$
$\begin{matrix} \\ \\ \end{matrix}-\begin{matrix}  \\ \\ \end{matrix}$

$D_2^{\star}$     ${1 \over 2} \omega_{3 \over 2}$
${3 \over 8}(1-\omega_{3 \over 2})$
${1 \over 4} (1-\omega_{3 \over 2})$
${1 \over 8} \omega_{3 \over 2}$
$\begin{matrix} \\ \\ \end{matrix} 0 \begin{matrix} \\ \\ \end{matrix} $

Probabilities for the population of the possible helicity states of $D$, $D^{\star}$, $D_1$ and $D_2^{\star}$.

Note that since $j^{(3)} = +1$ ($j^{(3)} = +2$) is only conceivable for the $D^{\star}$ ($D_2^{\star}$) state, whereas based on the initial assumption $j_Q^{(3)} = -1/2$, its correcponding probability would be zero.

To evaluate the parameter $\omega_j$, we consider the amplitude for the production of a pion at $\theta,\phi$ from a $H^{\star} \to H\pi$ type meson decay, which is proportional to the spherical harmonics $Y_{j}^{\ell}(\theta, \phi)$ ($\ell$ being the angular momentum quantum number of $H^{\star}$)

\begin{equation}
{d\Gamma(H^{\star} \to H\pi) \over d\cos\theta} \propto \int d\phi \sum_{j}  P_{H^{\star}}(j)  \bigl| Y_{j}^{\ell}(\theta,\phi) \bigr|^2,
\label{diff_width}
\end{equation}

where $P_{H^{\star}}(j)$ are the probabilities given in Table <ref>. Additionally, $\theta$ and $\phi$ are the angles of emission for the produced pion. Considering the case of $D_2^{\star} \to D\pi$, Eq. (<ref>) can be rewritten as

\begin{equation}
{1 \over \Gamma}{d\Gamma(D_2^{\star} \to D\pi) \over d\cos\theta} = {1 \over 4}
\left[ 1 + 3 \cos^2 \theta - 6 \; \omega_{3 \over 2} \left( \cos^2 \theta - {1 \over 3} \right)\right].
\label{diff_width_10}
\end{equation}

In Ref. [10], a comprehensive analysis over differential decay width Eq. (<ref>) and experimental data from [24, 25, 26, 27, 28] resulted in predicting $\omega_{3 \over 2} < 0.24$ with in 90% CL.

In order to implement the above arguments in the context of HQET and spin-flavour symmetry principles in , firstly we generalise Table <ref> for an arbitrary heavy quark helicity, $\rho_Q$:

$\hat{\rho}$     $\rho_{0,0}$     $\rho_{1,1}$     $\rho_{2,2}$     $\rho_{3,3}$     $\rho_{4,4}$

$D$     $1$     $-$     $-$     $-$     $\begin{matrix}  \\ \\ \end{matrix} -\begin{matrix}  \\ \\ \end{matrix}$

$D^{\star}$     ${1\over 2}(1-\rho_Q)$
${1\over 2}$
${1\over 2}(1+\rho_Q)$
$\begin{matrix}  \\ \\ \end{matrix} -\begin{matrix}  \\ \\ \end{matrix}$

$D_1$     ${1\over 16}[1-\rho_Q + \omega_{3 \over 2}(3-5\rho_Q)]$
${1 \over 4}(1-\omega_{3 \over 2})$
${1\over 16}[1-\rho_Q + \omega_{3 \over 2}(3+5\rho_Q)]$
$\begin{matrix} \\ \\ \end{matrix}-\begin{matrix}  \\ \\ \end{matrix}$

$D_2^{\star}$     ${1 \over 4} \omega_{3 \over 2} (1-\rho_Q)$
${3 \over 16}(1-\rho_Q) - {1 \over 8}\omega_{3 \over 2}(1-\rho_Q) $
${1 \over 4} (1-\omega_{3 \over 2})$
${3 \over 16}(1+\rho_Q) - {1 \over 8}\omega_{3 \over 2}(1+\rho_Q) $
$\begin{matrix} \\ \\ \end{matrix}
{1 \over 4} \omega_{3 \over 2} (1+\rho_Q)
\begin{matrix} \\ \\ \end{matrix} $

Possible polarization states of charmed mesons, $D$, $D^{\star}$, $D_1$ and $D_2^{\star}$.

Here, $\hat{\rho}$ is the diagonal polarization matrix of the produced mesons, with its components $\rho_{i,i}$, $i=0,1,2,3,4$ running through the most negative to most positive valid helicity states. The value $\omega_{3 \over 2} = 0.20$ is assumed as defualt.

In , the above implementation is remitted through the introduction of a so-called  transitional class. In particular, its  sub-module systematically handles the assignment of spin information and polarization to generated mesons based on their spin characteristics and heavy constituent quark flavours. Starting with checks on the meson's parentage and its constituents,  employs spin information from the heavy quark to construct the spin properties of the meson, pertinent to the excited heavy mesons (in addition to their complex conjugates):

(i) $D^{\star +}$, $D^{\star 0}$, $D_s^{\star +}$, $B^{\star 0}$, $B^{\star +}$, $B_s^{\star 0}$

(ii) $K_1^0$, $K_1^+$, $h'_1$, $D_1^+$, $D_1^0$, $D_{s1}^+$, $D_1^{' +}$, $D_1^{' 0}$, $D_{s1}^{' +}$, $B_1^0$, $B_1^+$, $B_{s1}^0$

(iii) $K^{\star 0}_2$, $K^{\star +}_2$, $f'_2$, $D^{\star +}_2$, $D^{\star 0}_2$, $D_{s2}^+$, $B_2^0$, $B_2^+$, $B_{s2}^0$

The above items correspond to the 2nd, 3rd and 4th row of Table <ref>, respectively. the  sub-module also calculates the spin polarization of the quark, updates the average polarization for the specific flavour and then assigns spin density matrix elements to the meson according to its spin value. It accommodates different spin types and meson categories, tailoring the spin behaviour and polarization calculations to match the distinct characteristics of various excited heavy mesons.

§.§ Polarization of Heavy Baryons

The ground state of a heavy baryon is characterized by a heavy quark combined with a $j_{qq} = 0$ helicity arrangement of a diquarks system. In this configuration, there is no angular momentum available to be transferred to the heavy quark, resulting in the preservation of the initial polarization without dilution. This implies that the polarization of the initial heavy quarks can indeed influence the ground state of heavy baryons. Here, the relative probabilities of finding these states in heavy sector fragmentation are still governed by two parameters, $\omega_a$ and $\omega_j$. The parameter $\omega_a$ represents the relative probability of producing a $j_{qq} = 1$ diquark as opposed to the ground state $j_{qq} = 0$ configuration. In the context of 's cluster hadronization model, we take $\omega_a = 1$ as in general no preference exists between spin-0 and spin-1 diquarks.

Now, similar to the case of excited heavy mesons, we consider the fragmentation of a $c$ quark with an assumed left-handed polarization. Given values of $\omega_a$ and $w_1$, the various helicity states of the charmed baryons can be populated by fragmentation
according to the following table [10]:

$j^{(3)}$     $-3/2$
$\begin{matrix}  \\ \\ \end{matrix} +3/2 \begin{matrix}  \\ \\ \end{matrix}$

$\Lambda_c$     $-$
${1\over 1+\omega_a}$
$\begin{matrix}  \\ \\ \end{matrix}-\begin{matrix}  \\ \\ \end{matrix}$

$\Sigma_c$     $-$
${(1-\omega_1) \omega_a \over 3 (1+\omega_a)}$
${ \omega_1 \omega_a \over 3 (1+\omega_a)}$
$\begin{matrix} \\ \\ \end{matrix}-\begin{matrix}  \\ \\ \end{matrix}$

$\Sigma_c^\star$     ${\omega_1 \omega_a \over 2 (1+\omega_a)}$
${2 (1-\omega_1) \omega_a \over 3 (1+\omega_a)}$
${\omega_1 \omega_a \over 6 (1+\omega_a)}$
$\begin{matrix} \\ \\ \end{matrix} 0 \begin{matrix} \\ \\ \end{matrix} $

Probabilities for the population of the possible helicity states of $\Lambda_c$, $\Sigma_c$ and $\Sigma_c^\star$.

Again, zero probabilities in Table <ref> are the helicity states that are viable but vanishing with the initial assumption of $j_Q^{(3)}=-1/2$. Additionally, the parameter $\omega_1$ can be estimated in a similar fashion as for $\omega_{3 \over 2}$ in Section <ref>, expanding and normalising Eq. (<ref>) for the observed decay modes $\Sigma_c \to \Lambda_c \pi$ and $\Sigma_c^{\star} \to \Lambda_c \pi$:

\begin{equation}
{1 \over \Gamma}{d\Gamma(\Sigma_c \to \Lambda_c \pi) \over d\cos\theta} = {1 \over 2}\ ,
\end{equation}

\begin{equation}
{1 \over \Gamma}{d\Gamma(\Sigma_c^{\star} \to \Lambda_c \pi) \over d\cos\theta} = {1 \over 4}
\left[ 1 + 3 \cos^2 \theta - {9 \over 2} \omega_{1} \left( \cos^2 \theta - {1 \over 3} \right)\right],
\end{equation}

that results in $\omega_{1} = {2/3}$.

The above arguments can be trivially generalised for arbitrary heavy quark helicity:

$\hat{\rho}$     $\rho_{0,0}$     $\rho_{1,1}$     $\rho_{2,2}$     $\rho_{3,3}$

$\Lambda_c$     ${1\over 2}(1-\rho_Q)$
${1\over 2}(1+\rho_Q)$
$\begin{matrix}  \\ \\ \end{matrix} -\begin{matrix}  \\ \\ \end{matrix}$

$\Sigma_c$     ${1\over 2}(1-\rho_Q) + \omega_{1} \rho_Q$
${1\over 2}(1+\rho_Q) - \omega_{1} \rho_Q$
$\begin{matrix} \\ \\ \end{matrix}-\begin{matrix}  \\ \\ \end{matrix}$

$\Sigma_c^\star$     ${3\over 8}\omega_{1}(1-\rho_Q)$
${1\over 2}(1-\rho_Q) - {1\over 8}\omega_{1}(3-5\rho_Q)$
${1\over 2}(1-\rho_Q) - {1\over 8}\omega_{1}(3+5\rho_Q)$
$\begin{matrix} \\ \\ \end{matrix}
{3\over 8}\omega_{1}(1+\rho_Q)
\begin{matrix} \\ \\ \end{matrix} $

Possible polarization states of charmed baryons, $\Lambda_c$, $\Sigma_c$ and $\Sigma_c^\star$.

Afterward, we are able to pass on the polarization information of the heavy baryons in the transitional , through the  sub-module. Commencing with checks on the baryon's parentage and its constituents, the module utilizes spin information from the heavy quark to establish the spin characteristics of the baryon, especially relevant to excited heavy baryons:

(i) $\Lambda_c^{+} ,\; \Lambda_c^{0} ,\; \Sigma_c^{++} ,\; \Xi_c^{0} ,\; \Xi_c^{+} ,\; \Xi_c^{'0} \; {\rm and \; C.C.}$

(ii) $\Sigma_c^{+} ,\; \Sigma_c^{0} ,\; \Sigma_c^{+} ,\; \Sigma_c^{0} ,\; \Xi_c^{+} ,\; \Xi_c^{0} ,\; \Omega_c^{0} ,\; \Omega_c^{0} \; {\rm and \; C.C.}$

(iii) $\Sigma_c^{++} ,\; \Xi_c^{++} ,\; \Xi_c^{+} ,\; \Omega_c^{+} ,\; \Xi_c^{+} ,\; \Xi_c^{++} ,\; \Omega_c^{++} ,\; \Omega_c^{+} ,\; \Xi_c^{+} \; {\rm and \; C.C.}$

The above items respectively correspond to rows 1 to 3 of Table <ref>.

§ DECAY OF THE EXCITED HEAVY MESONS

In the preceding section, we delineated a systematic methodology for transmitting the spin information of heavy hadron constituents, aiming to simulate the polarization states of excited heavy mesons and heavy baryons within , using the principles of HQEF and spin-flavour symmetry. It is imperative to highlight that, given the limited available experimental data on these hadrons, our reliance on HQET to determine their decay modes, widths, and branching ratios is considerable. Notably, HQEF does more than just earmark potential decay modes; it also provides a foundational basis for relationships among the associated couplings. This amplifies the importance of accurate modelling of these hadron decay modes. Historically, HQEF has been utilized for the strong and radiative decays of heavy baryons, as evidenced by works in references[34, 35] and its incorporation into  [1, 2]. Yet, there remains a pressing need to model the decays of excited heavy mesons to optimistically align predictions with the handful of extant experimental observations.

In this section, we detail the interactions specifically for charm mesons, though the driven results possess analogous applicability to the bottom meson system. Focusing initially on the $s-$ and $p-$wave mesons, we identify three meson multiplets under the auspices of heavy quark symmetry, as presented in [11]:

(i) The ground state ($J^P=0^-, 1^-$ doublet), encompassing $D$ and $D^{\star}$ mesons, characterized by $l=0$ for the light degrees of freedom.

(ii) The $J^P=1^+, 2^+$ doublet, consisting of $D_1$ and $D_2^{\star}$, distinguished by $l=1$ and $j=\frac32$ for the light degrees of freedom.

(iii) The $J^P=0^+, 1^+$ doublet, incorporating $D_0^{\star}$ and $D_1^\prime$, defined by $l=1$ and $j=\frac12$ for the light degrees of freedom.

Furthermore, one must consider the potential for mixing between $D_1$ and $D_1^\prime$, which arises due to sub-leading corrections within the heavy quark limit.

To ascertain the pertinent matrix elements associated with the decay of these charmed mesons, we invoke the methodology and notation detailed in Refs. [11, 36]. We aim to strike a sound balance between a cogent representation of the experimental data and the preservation of theoretical accuracy. Accordingly, we prioritize the leading terms within the heavy quark limit for the interactions, while also accommodating the mixing between the $D_1$ and $D_1^\prime$ mesons, $\theta_q$. This yields the following expression for the decay matrix elements:

\begin{align}
\mathcal{M}(D^{\star}\to D\pi) &= -{2g \over f_\pi} \left( m_Dm_{D^{\star}} \right)^{1\over 2} p_0\cdot \epsilon_0 ,
\label{m1}
\\
\mathcal{M}(D_2^{\star}\to D\pi) &= -{2h \over f_\pi \Lambda} \left( m_{D_2}m_D^{\star} \right)^{1\over 2} \epsilon_0^{\mu\nu} p_{0,\mu} p_{0,\nu} ,
\label{m2}
\\
\mathcal{M}(D_2^{\star}\to D^{\star}\pi) &= -i{2h \over f_\pi \Lambda} \left( {m_{D^{\star}} \over m_{D_2}} \right)^{1\over 2}  \epsilon^{\alpha\beta\mu\nu} \epsilon^0_{\alpha\gamma} p_0^\gamma p_{0,\mu} p_{1\nu}\epsilon_{1\beta} ,
\label{m3}
\\
\mathcal{M}(D_1\to D^{\star}\pi) &=  {h \over f_\pi\Lambda}   \left( {2 \over 3} {m_{D_1}m_D} \right)^{1\over 2}\left[ \epsilon_0\cdot\epsilon_1\left(p_0^2-\left[{p_0\cdot p_1 \over m_0}\right]^2\right)-3\epsilon_0\cdot p_0\epsilon_1\cdot p_0\right] ,
\label{m4}
\\
\mathcal{M}(D_0^{\star}\to D\pi) &= {f^{\prime\prime} \over f_\pi}\left( {m_{D_0^{\star}}m_D} \right)^{1\over 2} \,p_0\cdot\left({p_1 \over m_{D^{\star}_0}}+{p_2 \over m_D}\right) ,
\label{m5}
\\
\mathcal{M}(D_1^\prime \to D^{\star}\pi) &=-{f^{\prime\prime} \over f_\pi}\left( {m_{D_1^\prime}m_D} \right)^{1\over 2}\left[-p_0\cdot\left({p_1 \over m_{D^{\star}_0}}+{p_2 \over m_D}\right)\epsilon_0\cdot\epsilon_1\right.
\nonumber\\
&  \left. +{1 \over m_{D'_1}}\epsilon_1\cdot p_1\epsilon_0\cdot p_0
+{1 \over m_{D}}\epsilon_0\cdot p_2\epsilon_1\cdot p_0 \right].
\label{m6}
\end{align}

Here, $p_i$ and $\epsilon_i$ are the momenta and the polarization vectors of the hadrons, with $i=0$ pointing to the parent hadron and $i=1,2$ to the heavy and light child hadrons respectively. $m_H$ is the mass of hadron $H$ and $g$, $h$, $\Lambda$, $f_{\pi}$ and $f^{''}$ are decay parameters. We also introduce $\theta_q$ as the mixing angles between ($D_1$, $D'_1$) and ($D_{s1}$, $D'_{s1}$) mesons.

Having the matrix elements of the decays, we can calculate the partial widths using the Fermi golden rule for two-body decays,

\begin{equation}
\Gamma(H^{\star} \to H \pi) = {\Theta(H^{\star} \to H \pi) \over 2m_{H^{\star}}}
\left| \mathcal{M}(H^{\star} \to H \pi) \right|^2,
\end{equation}

with $\Theta$ being the two-body phase-space factor. Furthermore, to ensure consistency between the theoretical calculations and our implementation in , it becomes necessary to retain some of the sub-leading terms when considering the heavy quark limit. Explicitly, we use the following prescription for our calculations:

\begin{equation}
\Gamma(H^{\star} \to H \pi) = {1 \over 8 \pi \; m_{H^{\star}}^2}
\left| \mathcal{M}(H^{\star} \to H \pi) \right|^2 p_{\rm CM},
\end{equation}

where $p_{\rm CM}$ denotes the momentum of the parent hadron in the kinematical centre-of-mass frame of the two-body decay. The resulting partial widths are outlined below:

\begin{align}
\Gamma (D^{\star}\to D\pi) & =  \frac{g^2}{6\pi f^2_\pi}  \frac{m_D}{m_{D^{\star}}} p_{\text{cm}}^3 \\
\Gamma (D_2^{\star}\to D^{\star}\pi) & = \frac{h^2}{15\pi f^2_\pi\Lambda^2} \frac{m_{D^{\star}}}{m_{D^{\star}_2}} p_{\text{cm}}^5\\
\Gamma (D_2^{\star}\to D\pi) &=  \frac{h^2}{10 \pi f^2_\pi\Lambda^2}  \frac{m_{D}}{m_{D^{\star}_2}}  p_{\text{cm}}^5\\
\Gamma (D_1\to D^{\star}\pi) &=  \frac{h^2}{144 \pi f^2_\pi\Lambda^2}  \frac{\left[-2m^2_{D^1} (m^2_{D^{\star}}-5 m_\pi^2)+(m_\pi^2-m^2_{D^{\star}})^2+25 m^4_{D^1}\right]}{m_{D^{\star}}m^3_{D^1}}   p_{\text{cm}}^5  \\
\Gamma (D_0^{\star}\to D\pi) &= \frac{ f^{\prime\prime2}}{32 f_\pi^2\pi}         \frac{\left(m_{D^{\star}_0}-m_{D}\right)^2     }{m_{D}m^3_{D^{\star}_0}} \left[ \left(m_{D^{\star}_0}+m_{D}\right)^2 -m_\pi^2\right]^2  p_{\text{cm}} \\
\Gamma (D_1^\prime \to D^{\star}\pi) &= \frac{ f^{\prime\prime2}}{32 f_\pi^2\pi }  \frac{\left(m_{D^\prime_1}-m_{D^{\star}}\right)^2}{m_{D^{\star}}m^3_{D^\prime_1}}  \left[ \left(m_{D^\prime_0}+m_{D^{\star}}\right)^2 -m_\pi^2\right]^2p_{\text{cm}}
\end{align}

Note that in  uses both the amplitude and the decay width expressions simulation of the decay modes to form a precise phase space.

To extract the numerical values of the decay parameters $g$, $h$, $\Lambda$, $f_{\pi}$ and $f^{''}$ in Eqs. (<ref>) through (<ref>), we make use of the latest observed masses and decay widths of the charmed mesons from BaBar [37, 29] and LHCb [38] collaborations. The complete list of utilized masses and widths for the $(0^-,1^-)$, $(1^+,2^+)$ and $(0^+,1^+)$  multiplets are given in Table <ref>.

Multiplet     State     Meson     Mass [GeV]     Width     $\Delta m$

$(0^-,1^-)$     $0^-$     $D^+$          $ 1.8697 \pm  0.0001$     N/A     $   4.82 \pm    0.07$

$D^0$          $ 1.8648 \pm  0.0001$     N/A     $   0.00 \pm    0.00$

$D^+_s$        $ 1.9683 \pm  0.0001$     N/A     $ 103.51 \pm    0.09$

$1^-$     $D^{*+}$       $ 2.0103 \pm  0.0001$     $( 0.0834 \pm  0.0018)\times10^{-3}$     $   3.41 \pm    0.07$

$D^{*0}$       $ 2.0069 \pm  0.0001$     $< 0.0021$     $   0.00 \pm    0.00$

$D^{*+}_s$     $ 2.1122 \pm  0.0004$     $< 0.0019$      $ 105.35 \pm    0.40$

$(1^+,2^+)$     $1^+$     $D_1^+$          $ 2.4232 \pm  0.0024$     $ 0.0250 \pm  0.0060$     $   2.40 \pm    2.45$

$D_1^0$          $ 2.4208 \pm  0.0005$     $ 0.0317 \pm  0.0025$     $   0.00 \pm    0.00$

$D_{s1}^+$       $ 2.5351 \pm  0.0001$     $ 0.0009 \pm  0.0001$     $ 114.31 \pm    0.50$

$2^+$     $D_2^{*+}$          $ 2.4654 \pm  0.0013$     $ 0.0467 \pm  0.0012$     $   4.70 \pm    1.36$

$D_2^{*0}$          $ 2.4607 \pm  0.0004$     $ 0.0475 \pm  0.0011$     $   0.00 \pm    0.00$

$D_{s2}^{*+}$       $ 2.5691 \pm  0.0008$     $ 0.0169 \pm  0.0007$     $ 108.40 \pm    0.89$

$(0^+,1^+)$     $0^+$     $D_0^+$          $ 2.3490 \pm  0.0070$     $ 0.2210 \pm  0.0180$     $  49.00 \pm   20.25$

$D_0^0$          $ 2.3000 \pm  0.0190$     $ 0.2740 \pm  0.0400$     $   0.00 \pm    0.00$

$D_{s0}^+$       $ 2.3178 \pm  0.0005$     $< 0.0038$     $  17.80 \pm   19.01$

$1^+$     $D_1^{\prime+}$     2c|Not observed

$D_1^{\prime0}$          $ 2.4270 \pm  0.0400$     $   0.38^{+   0.13}_{-   0.11}$     $   0.00 \pm    0.00$

$D_{s1}^{\prime+}$       $ 2.4595 \pm  0.0006$     $< 0.0035$     $  32.50 \pm   40.00$

Observed masses and widths of the charm mesons, from BaBar [37, 29] and LHCb [38] collaborations.

Figure <ref> shows the helicity angle distributions for the $D_{s1}(2536)\to D^{*+}K^0$ and $D_{1}^0(2420)\to D^{*+}\pi^-$ decays modes. $s$- and $d$-wave contributions are shown independently and the mixed sums are fitted against the corresponding experimental data. The fitting is performed using  python module [39] to find the best values for the $g$, $h$, $\Lambda$, $f_{\pi}$, $f^{''}$ and $\theta_q \;(q=u,d,s)$ parameters. These best-fit values are recorded in Table <ref>.

Helicity Angle distributions for the decays: (left panel) $D_{s1}(2536)\to D^{*+}K^0$ and (right panel) $D_{1}^0(2420)\to D^{*+}\pi^-$.

Parameter          Fitted Value
$f''$              $ -0.465 \pm 0.017$
$f_\pi$            $  0.130 \pm 0.001$ [GeV]
$h$                $  0.824 \pm 0.007$
$\Lambda$          $  1.000 \pm 0.000$ [GeV]
$g$                $  0.565 \pm 0.006$
$\theta_{u,d}$     $  0.000 \pm 0.100$
$\theta_s$         $ -0.047 \pm 0.002$

Fitted values of the decay parameters $g$, $h$, $\Lambda$, $f_{\pi}$, $f^{''}$ and $\theta_q \;(q=u,d,s)$, using the charm mesons masses and decay widths from Refs. [37, 29, 38].

Besides the strong isospin-conserving decays, certain excited mesons also undergo electromagnetic decays, as well as strong-isospin-violating ones. These latter decays are particularly prominent in cases where the strong isospin-conserving decays are either kinematically inhibited or strongly suppressed, either because of the lack of sufficient rest energy, threshold effects, angular momentum conservation, particular selection rules in the given event handler, hierarchy of coupling constants or presence of other dominant channels. These are predominantly observed in the following cases:

(i) $D^{\star}$ mesons: Here, the strong isospin-conserving decays are kinematically limited, making radiative modes crucial.

(ii) $B^{\star}$ mesons: In this case, the strong isospin-conserving decays cannot occur due to kinematic constraints. Thus, only the radiative mode emerges as a possibility.

(iii) $D^{+}_s$, $D^{+}_{s0}$, and $D^{+}_{s1}(2460)$ mesons: For these, both radiative and isospin violating decay modes hold significance since the strong isospin conserving $DK$ modes are kinematically proscribed.

(iv) $B_s^{*0}$ meson: Only the radiative mode is feasible from a kinematic perspective.

Again, our focus will be on the $D^{\star}$ system, as it represents the most intricate scenario, exhibiting a plethora of excited mesons below the strong decay threshold. Nevertheless, our arguments can be effortlessly extended to other cases.

The amplitude for the radiative decay of the $D^{\star}$ mesons, as detailed in [40], is given by:

\begin{equation}
\mathcal{M}(D^{\star}\to D\gamma) = A \left[ 64\pi \frac{m_D}{m_{D^{\star}}} \right]^{1 \over 2} \epsilon^{\alpha\beta\mu\nu}\epsilon_{\gamma\alpha}p_{\gamma\beta}\epsilon_{D^{\star}\mu}p_{D^{\star}\nu}.
\end{equation}

Here, the coupling, $A$, is expressed as:

\begin{equation}
A= \frac{e_Q}{4m_Q}\alpha(m_Q)^{\frac{1}{2}} +\frac{c_H}{\Lambda}e_q\alpha(\Lambda)^{\frac12},
\end{equation}

with $e_Q$ and $e_q$ being the electric charges of the heavy and light quarks respectively $m_Q$ being the mass of the heavy quark. $\alpha$ is the electromagnetic coupling constant while $c_H = -1.058$ is electromagnetic coefficient for heavy meson decays. The scale $\Lambda$ is the same as the case of strong decays. Consequently, the partial width becomes [41, 42]:

\begin{equation}
\Gamma(D^{\star}\to D\gamma) = \frac{2 m_D|A|^2 }{3 m_{D^{\star}}} \left(\frac{m_{D^{\star}}^2-m_D^2}{m_{D^{\star}}}\right)^3.
\end{equation}

To implement the strong and radiative decays of excited heavy mesons in , we introduced two specialized classes:  and . The decay parameters highlighted in this section are defined as user-adjustable variables. This design provides flexibility, facilitating potential tuning and refinement based on future insights or requirements. In the subsequent section, we will evaluate the robustness of our implementation by juxtaposing it against the available experimental data.

§ EXAMINATION OF THE OBSERVABLES

The integration of the  into  enabled more accurate polarization details for excited heavy mesons and baryons, based on their constituent heavy quarks' spin. Additionally, the  and  classes were incorporated, further enhancing 's ability to model both strong and radiative decays, thereby providing a comprehensive framework for heavy hadron decay dynamics. To evaluate these improvements, we scrutinized existing experimental data, particularly those that are sensitive to hadron polarization states. Limited but valuable insights were gained from measurements like the angular decomposition of the $D_{s1}^+ \to D^{\star +} K_{s0}$ decay mode, as observed by BaBar [29] and BELLE [30] and from $\Lambda_b$ baryon polarization in hadronic $Z^0$ decays by ALEPH [43] and OPAL [32].

Efficiency-corrected decay rates of the $D_{s1}^+$ meson: (left panel) displays the $D_{s1} \to D^{\star +} K^0$ decay mode rates as a function of lab-frame angle $\theta$. (right panel) shows the same but in the $D_{s1}^+$ center-of-mass frame, denoted by $\theta'$. Data sourced from Ref. [29]. Red histograms represent  predictions with only the  active, while blue histograms indicate the combined effect of  and  classes. Green histograms correspond to  predictions, which lack HQET enhancements.

Further analysis of efficiency-corrected decay rates for $D_{s1}^+$: (left panel) rates of $D_{s1} \to D^{\star +} K^0$ as a function of lab-frame angle $\alpha$. (right panel) the same rates with respect to angle $\beta$. (bottom panel) focuses on the angle $\gamma$ between $\pi^+$ and $K^0$ in the $D^{\star +}$ rest frame. The data is sourced from [30] and the histograms are annotated as in Figure <ref>.

In Figures <ref> and <ref>, we investigate the efficacy of using HQET and spin-flavor symmetry to predict $e^-e^+$ data concerning polarization-sensitive measurements of the $D_{s1}$ meson decays. The figures display normalized, efficiency-corrected rates for the $D_{s1} \to D^{\star +} K^0$ decay mode as functions of different decay angles. We evaluate the contributions of the  and  classes in , activating them individually to study their distinct effects on the predictions. For comparison, we also include predictions from the previous version, , which lacks both spin-flavor symmetry and HQET enhancements. As anticipated,  fails to accurately predict either the existence or the behaviour of angle-dependent decay rates. On the other hand,  correctly predicts the mean values of the $s$-wave contributions when only the  class is activated, but does not capture the full angle-dependent behaviour. Importantly, with both the  and  classes enabled,  achieves an excellent agreement with the observed data.

Average polarization of $\Lambda_b$ baryons in hadronic $Z^0$ decays at LEP by ALEPH [43] (left panel) and OPAL [32] (right panel) for $\sqrt{s}=91.2$ GeV. Histograms are annotated as in Figure <ref>.

For heavy baryons, we validate our implementation of HQET and spin-flavor symmetry by comparing with measurements of the average polarization of $\Lambda_b$ baryons in hadronic $Z^0$ decays at LEP, as reported by the ALEPH [43] and OPAL [32] collaborations, shown in Figure <ref>. Notably, although the predictions from  are small yet distinct from zero, the combined use of  and  in  leads to predictions nearly an order of magnitude larger, while consistent with the experimental uncertainty.

After incorporating modifications to the hadronization and decays of heavy excited mesons and heavy baryons in , it is advisable to revisit the production rates of these hadrons in the simulation. This ensures that the framework remains consistent, particularly if the new modifications could significantly influence kinematic distributions or other dependent observables. With this prospect in mind, we initially undertake a general tuning effort, trying to optimise the relevant existing and newly defined parameters to sin synergy with a large pool of experimental data. In the next section, we will describe our methodology for this general tune and its outcome, before re-evaluating the production rates of the heavy and excited hadrons for a rigorous validation process.

§ GENERAL TUNE FOR HERWIG-7.3.0

For the first general tune of  since the 7.2.0 release, we utilized $e^-e^+$ data from LEP, PETRA, SLAC, SLC and TRISTAN, spanning just over 9,200 separate data bins. The data were weighted around both light and heavy hadron production rates and multiplicities, with a focus on more dominant processes. Initial attempts to use  yielded inconsistencies due to the complexity of the tuning and the high number of parameters. We then resorted to a multi-layered, brute-force approach using the  module from , minimizing the $\chi^2$ as an indicator of the best tune.

In this context, we analyzed a total of 12 parameters, comprising 10 for cluster hadronization and 2 for the parton shower—specifically  and . Below is a detailed breakdown of these parameters:

* : Maximum allowable cluster mass for light quarks.

* : Power exponent for the mass of clusters with light quarks.

* : Parameter affecting the mass splitting for clusters with light quarks.

* : Probability for a $s\bar{s}$ quark pair to be spawned during cluster splittings.

* : Probability for quarks forming a di-quark.

* : Weighting factor for singlet baryons in hadronization.

* : Weighting factor for decouplet baryons in hadronization.

* : Exponential factor in the  probability function.

* : Offset in the  probability function.

* : Adjustment to the kinetic threshold in .

* : Initial value for the strong coupling constant at $M_{Z^0}=91.1876$ GeV[The reader should note that after the  release, the  has been split into two separate parameters,  and  that control the strong coupling constants for initial- and final-state radiations respectively. Henceforth, the tuned value we suggest in Table <ref> points to , as it is more relevant to our $e^+e^-$ tuning strategy. Meanwhile, the tune proposed in [44, 45] has been adopted for the intial-state radiation coupling constant =$0.1185$.].

* : Minimum transverse momentum in parton shower.

Notably, the parameters , , and  are introduced only in . The first two parameters belong to the newly introduced  Probability Function, $P_{\rm cluster}$, which serves as a decision-making utility that evaluates whether a given cluster meets a dynamically calculated kinematic threshold:

\begin{equation}
P_{\rm cluster} = \frac{1}{1 + \left| \frac{M - \delta}{M_{\rm th}} \right|^r}, \qquad P_{\rm cluster} > P_{\rm random}(0,1),
\end{equation}

Here, $M$ denotes the effective mass of the cluster while $M_{\rm th}$ is the mass threshold for the cluster, defined as the sum of the masses of the constituent quarks plus the mass of the spawned di-quark. $\delta$ and $r$ are the  and  parameters, respectively. $P_{\rm cluster}$ is then compared against a randomly generated number between 0 and 1, denoted as $P_{\rm random}(0,1)$, to either accept or reject the cluster splitting. This new formulation provides a smooth Gaussian distribution for dynamic threshold cuts based on the cluster mass and its kinematic properties. Furthermore, the  kinematic threshold for 2-body cluster decay in  was:

\begin{equation}
M > M_1 + M_2, \quad M_1 >  m + m_1, \quad M_2 > m + m_2,
\end{equation}

This has been updated to a new,  threshold in :

\begin{equation}
M^2 > M_1^2 + M_2^2 , \quad M_1^2 > m^2 + m_1^2 + \delta_{\rm th}, \quad M_2^2 > m^2 + m_2^2 + \delta_{\rm th}.
\end{equation}

In this new formulation, $M_1$ and $M_2$ are the masses of the child clusters, $m_1$ and $m_2$ are the masses of the constituent quarks of the parent cluster, and $m$ represents the mass of the spawned quark-antiquark pair. $\delta_{\rm th}$ is the  parameter. Note that both  and dynamic threshold options are available in , with the default being the  choice:

Following the above description and methodology, we tuned the hadronization and parton shower parameters in . The results are given in Table <ref>.

Tuned Parameter

[GeV]    3.529     3.649

1.849     2.780

0.914     0.899

0.374     0.292

0.331     0.298

0.891     0.740

0.416     0.620

6.486     $-$

-0.879     $-$

[GeV]    0.088     $-$

0.102     0.126

[GeV]    0.655     0.958

The values of tuned parameters in  compared to .

Despite its computational intensity, our brute-force approach can be deemed a reasonable success, achieving a substantial improvement in the overall $\chi^2$ values. Specifically, the tuned version of  shows a $\sim$50.75% improvement in $\chi^2$ compared to its untuned counterpart and a $\sim$12.76% improvement when compared to . In the following section, we will look into the effects of this new tune on the production rates of heavy mesons and baryons.

§ HEAVY HADRON PRODUCTION RATES

In the realm of heavy hadron experimental physics, most of the available data on the production of excited states primarily stems from experiments focusing on the production of $B$ mesons from the $\Upsilon(4S)$ resonance, captured through data in the continuum region, i.e. at an energy range just above the production threshold for $B$ mesons but below the $\Upsilon(4S)$ peak, where $B$ mesons can be produced directly through $e^+e^-$ annihilation without necessarily forming a $\Upsilon$ state[This is in contrast to the “resonance region", which is the energy range near the $\Upsilon(4S)$ peak, where the $e^+e^-$ collision is likely to form an $\Upsilon(4S)$  state before subsequently decaying into $B$ mesons.]. To thoroughly evaluate these phenomena, we undertook a substantial effort to collect and compile all accessible experimental data [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80] into a single  analysis, see Appendix <ref>. This comprehensive dataset served as a benchmark to scrutinize the production rates of all heavy mesons and baryons in . Particular attention was devoted to iso-spin multiplets; we meticulously calibrated their production rates to be nearly identical, thereby ensuring a robust and accurate hadronization and decay framework within the  simulation.

Comparison of heavy meson production rates, collated from various experimental data sets [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59]. The observed rates are compared against the predictions from both the tuned (represented by red lines) and untuned (depicted by blue lines) versions of . For reference, the results from  are indicated by green lines.

Comparison of heavy baryon production rates, collated from various experimental data sets [60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78]. Histograms are annotated as in Figure <ref>.

Figures <ref> and <ref> illustrate the impact of incorporating our spin-flavour symmetry approach into  and HQET hadronic decay prescriptions. These figures compare the production rates of heavy mesons and baryons, as predicted by both the tuned and untuned versions of . For a comprehensive view, we also include corresponding results from . A notable enhancement in predictive accuracy is evident in the tuned  simulations.

The improvements introduced in the previous sections, combined with the new general tune of  have made tangible improvements in the predictive modelling of heavy hadron production. This updated framework is poised to serve as an invaluable asset for both theoretical and experimental research in the high-energy physics community, especially for those focused on the complex dynamics of heavy hadron formation.

§ CONCLUSION

In this study, we have made substantial advancements in the predictive modelling of heavy hadron production through targeted modifications to the  event generator. Specifically, we have overhauled its native cluster hadronization model by integrating HQET and spin-flavour symmetry. This was realized by adding three key modules: the  class, which delivers polarization information for excited heavy mesons and heavy baryons based on their constituent heavy quarks' helicity states; and the  and  classes, which improve the handling of heavy meson decays by enhancing their two-body decay matrix elements. This effectively complements the pre-existing heavy baryon decayers. Additionally, we formulated a new dynamic kinematic regime and included a decision-making probability function for smoother cluster splitting distributions, which improves the overall handling power of  in the domain of heavy hadron production.

In the next phase of our efforts, we tuned  using a brute-force, multilayered approach aimed at minimizing the $\chi^2$ statistical parameter. To facilitate this, we assembled an extensive collection of experimental data tailored to our specific objectives. This data-driven methodology allowed us to fine-tune 12 parameters concurrently, including two parton shower-related parameters essential for the integration of the newly introduced electroweak parton shower. Consequently, we have achieved a notable elevation in 's capacity to make accurate predictions, particularly for the complex case of heavy hadron production.

The inclusion of HQET and spin-flavour symmetry, underpinned by robust theoretical foundations, serves as a significant milestone in improving the predictive capabilities of the cluster hadronization model. These methodological enhancements are poised to become invaluable assets for both theoretical and experimental research in high-energy physics, especially in areas requiring a nuanced understanding of heavy hadron dynamics and precision measurements. Given the notable enhancements in predictive accuracy, we strongly recommend the utilization of  in future high-precision experimental analyses. We are pleased to announce that these improvements will be made publicly available in the forthcoming  release, underscoring its potential as a reliable computational tool for the broader high-energy physics community.

§ ACKNOWLEDGEMENTS

We thank our fellow Herwig authors for useful discussions. This work has received funding from the European Union’s Horizon 2020 research and innovation program as part of the Marie Skłodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). MRM is also supported by the UK Science and Technology Facilities Council (grant numbers ST/P001246/1).

§ SIGNITURE OF HEAVY HADRONS IN THE $\UPSILON(4S)$ CONTINUUM REGION

We have compiled all existing experimental data on heavy mesons and baryons in the $\Upsilon(4S)$ “continuum region" [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80], where $B$ mesons are produced directly from $e^+e^-$ annihilation, bypassing the formation of an $\Upsilon(4S)$ state. This is distinct from the “resonance region", where $e^+e^-$ collisions are more likely to form an $\Upsilon(4S)$ state that then decays into $B$ mesons. This data was utilized to study the effects of employing the HQET and the spin-flavour symmetry on the overall balance of heavy hadron production rate in  and have been particularly beneficial in our general tuning efforts.

The tables that follow offer a detailed look at various measurements. Specifically, Tables <ref> and <ref> display results for non-strange and strange heavy charmed mesons, respectively. Tables <ref> and <ref> focus on the experimental data for non-strange and strange charm baryon production in the $\Upsilon$ energy region. Lastly, Table <ref> provides insights into the observations of charmonium states in the continuum region.

Hadron     Energy/GeV     Experiment     Observable     Notes

10.47     ARGUS [49]     $\sigma$

10.52     BELLE  [56]     $\sigma$

10.5     CLEO [55]     $\sigma$

10.47     ARGUS [49]     $\sigma$

10.52     BELLE  [56]     $\sigma$

10.5     CLEO [55]     $\sigma$

10.52     BELLE [56]     $\sigma$

10.5     CLEO [55]     $\sigma$

10.47     ARGUS [49]     $\sigma$

10.52     BELLE [56]     $\sigma$

10.5     CLEO [55]     $\sigma$

10.58     ARGUS [47]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.55$

10.58     CLEO [54]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

10.58     ARGUS [46]     $\sigma\times\text{Br}$

10.58     ARGUS [47]     $\sigma\times\text{Br}$

10.58     CLEO [53]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.55$

10.58     CLEO [54]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

Measurements of non-strange charm meson production in the $\Upsilon$ energy region.

Hadron     Energy/GeV     Experiment     Observable     Notes

5*$D_s^+$     10.58     BaBar [57]     $\sigma$

10.52     BELLE [56]     $\sigma$

10.5     CLEO [58]     $\text{Br}\times\text{d}\sigma/\text{d}x_p$

3*$D_s^{*+}$     10.58     BaBar [57]     $\sigma$

10.5     CLEO [58]     $\text{Br}\times\text{d}\sigma/\text{d}x_p$

2*$D_{s2}^+$     10.58     ARGUS [52]     $\sigma\times\text{Br}$     for $x_p>0.7$ + extrapolation

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.6$

2*$D_{s0}^{\star}(2317)^+$     10.6     BaBar [59]     $\sigma\times\text{Br}$     for $p*>3.2$ GeV

$1/N\text{d}N/\text{d}p^{\star}$     for $p*>3.2$ GeV

2*$D_{s1}^{\star}(2460)^+$     10.6     BaBar [59]     $\sigma\times\text{Br}$     for $p*>3.2$ GeV

$1/N\text{d}N/\text{d}p^{\star}$     for $p*>3.2$ GeV

3*$D_{s1}^{\star}(2536)^+$     10.6     ARGUS [48]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

10.4     ARGUS [51]     $\sigma\times\text{Br}$

2*$D_{s1}^{\star}(2536)^+$     10.6     BaBar [59]     $\sigma\times\text{Br}$     for $p*>3.2$ GeV

10.6     CLEO [50]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

Measurements of strange charm meson production in the $\Upsilon$ energy region.

Hadron     Energy/GeV     Experiment     Observable     Notes

8*$\Lambda_c^+$     10.54     BaBar [66]

10.52     BELLE [56]

10.52     BELLE [78]     $\sigma$

10.5     CLEO [70]     $\text{Br}\times\text{d}\sigma/\text{d}x_p$

5*$\Lambda_c(2595)^+$     10.4     ARGUS [61]
$\sigma\times\text{Br}$     for $x_p>0.7$ + extrapolation

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

10.6    CLEO [74]     $\sigma/\sigma(\Lambda_c^+)$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

10.52     BELLE [78]     $\sigma$

4*$\Lambda_c(2625)^+$     10.4     ARGUS [62]
$\sigma\times\text{Br}$     extrapolation to whole region

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

10.6    CLEO [74]

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

4*$\Sigma_c^0$     10.58      ARGUS [63]     $\sigma/\sigma(\Lambda_c^+)$

$1/N\text{d}N/\text{d}x_p$     averaged over $\Sigma_c^{0,++}$

10.52     BELLE [78]     $\sigma$

3*$\Sigma_c^{++}$     10.58      ARGUS [63]     $\sigma/\sigma(\Lambda_c^+)$     averaged over $\Sigma_c^{0,++}$

$1/N\text{d}N/\text{d}x_p$     averaged over $\Sigma_c^{0,++}$

10.58     CLEO [72]      $\sigma/\sigma(\Lambda_c^+)$     averaged over $\Sigma_c^{0,++}$

2*$\Sigma_c^{*0}$     10.52     BELLE [78]     $\sigma$

2*$\Sigma_c^{*++}$     10.5    CLEO [73]    $\sigma/\sigma(\Lambda_c^+)$     averaged over $\Sigma^{*0,++}_c$

$\text{d}\sigma/\text{d}x_p$     $x_p>0.5$ averaged over $\Sigma^{*0,++}_c$

2*$\Sigma_c(2800)$     10.52     BELLE [77]
$\sigma\times\text{Br}$     extrapolation to whole region

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

Measurements of non-strange charm baryon production in the $\Upsilon$ energy region.

Hadron     Energy/GeV     Experiment     Observable     Notes

2*$\Xi_c^0$     10.5     ARGUS [60]
$\sigma\times\text{Br}$     for $x_p>0.5$ + extrapolation

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.55$ averaged over $\Xi_c^{0,+}$

4*$\Xi_c^0$     10.58     BaBar [65]

10.52     BELLE [78]     $\sigma\times\text{Br}$

2*$\Xi_c^+$     10.5     ARGUS [60]
$\sigma\times\text{Br}$     for $x_p>0.5$ + extrapolation

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.55$ averaged over $\Xi_c^{0,+}$

4*$\Xi_c^{\prime0}$     10.5     BaBar [67]
$\sigma\times\text{Br}$      Br is for the $\Xi_c^0$ decay

$1/N\text{d}N/\text{d}x_p$     also includes $\Upsilon(4S)$ decays

CLEO [76]     $\sigma/\sigma(\Xi_c^0)$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$ averaged over $\Xi_c^{\prime0,+}$

4*$\Xi_c^{\prime+}$     10.5     BaBar [67]
$\sigma\times\text{Br}$      Br is for the $\Xi_c^+$ decay

$1/N\text{d}N/\text{d}x_p$     also includes $\Upsilon(4S)$ decays

10.5     CLEO [76]     $\sigma/\sigma(\Xi_c^+)$

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$ averaged over $\Xi_c^{\prime0,+}$

2*$\Xi_c^{*+}$     10.6     CLEO [75]

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

2*$\Xi_c^{*0}$     10.6     CLEO [71]

$1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

1*$\Xi_{c1}^+$     10.6     CLEO [64]     $1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

1*$\Xi_{c1}^0$     10.6     CLEO [64]     $1/N\text{d}N/\text{d}x_p$     for $x_p>0.5$

4*$\Omega_c^0$      10.52     BELLE [78]     $\sigma\times\text{Br}$

10.5     BaBar [69]     $\sigma\times\text{Br}$

$1/N\text{d}N/\text{d}x_p$     also includes $\Upsilon(4S)$ decays

1*$\Omega_c^{*0}$     10.5     BaBar [68]     $\sigma/\sigma(\Omega_c)$     $x_p>0.5$

Measurements of strange charm baryon production in the $\Upsilon$ energy region.

Hadron     Energy/GeV     Experiment     Observable

3*$J/\psi$     10.57     BaBar [79]     $\sigma$

10.6      BELLE [80]     $\sigma$

2*$\psi(2S)$     10.6      BELLE [80]     $\sigma$

Measurements of charmonium production in the $\Upsilon$ energy region.

[1]
M. Bahr, S. Gieseke, M. A. Gigg, D. Grellscheid, K. Hamilton, O. Latunde-Dada, S. Platzer, P. Richardson, M. H. Seymour and A. Sherstnev, et al.

Eur. Phys. J. C 58 (2008), 639-707.

[2]
J. Bellm, S. Gieseke, D. Grellscheid, S. Plätzer, M. Rauch, C. Reuschle, P. Richardson, P. Schichtel, M. H. Seymour and A. Siódmok, et al.

Eur. Phys. J. C 76 (2016) no.4, 196.

[3]
J. Bellm, S. Gieseke, D. Grellscheid, P. Kirchgaeßer, F. Loshaj, G. Nail, A. Papaefstathiou, S. Plätzer, R. Podskubka and M. Rauch, et al.

[arXiv:1705.06919 [hep-ph]].

[4]
J. Bellm, G. Bewick, S. Ferrario Ravasio, S. Gieseke, D. Grellscheid, P. Kirchgaeßer, M. R. Masouminia, G. Nail, A. Papaefstathiou and S. Platzer, et al.

Eur. Phys. J. C 80 (2020) no.5, 452.

[5]
A. Kupco,

[arXiv:hep-ph/9906412 [hep-ph]].

[6]
A. F. Falk, H. Georgi, B. Grinstein and M. B. Wise,

Nucl. Phys. B 343 (1990), 1-13.

[7]
M. Neubert,

Phys. Rept. 245 (1994), 259-396.

[8]
A. F. Falk, M. E. Luke and M. J. Savage,

Phys. Rev. D 49 (1994), 3367-3378.

[9]
I. I. Y. Bigi, M. A. Shifman, N. G. Uraltsev and A. I. Vainshtein,

Phys. Rev. D 50 (1994), 2234-2246.

[10]
A. F. Falk and M. E. Peskin,

Phys. Rev. D 49 (1994), 3320-3332.

[11]
A. F. Falk and M. E. Luke,

Phys. Lett. B 292 (1992), 119-127.

[12]
M. R. Masouminia and P. Richardson,

JHEP 04 (2022), 112.

[13]
N. Darvishi and M. R. Masouminia,

Nucl. Phys. B 985 (2022), 116025.

[14]
A. Buckley, J. Butterworth, D. Grellscheid, H. Hoeth, L. Lonnblad, J. Monk, H. Schulz and F. Siegert,

Comput. Phys. Commun. 184 (2013), 2803-2819.

[15]
A. Buckley, H. Hoeth, H. Lacker, H. Schulz and J. E. von Seggern,

Eur. Phys. J. C 65 (2010), 331-357.

[16]
E. V. Shuryak,

Phys. Lett. B 93 (1980), 134-136.

[17]
W. E. Caswell and G. P. Lepage,

Phys. Lett. B 167 (1986), 437-442.

[18]
E. Eichten,

Nucl. Phys. B Proc. Suppl. 4 (1988), 170.

[19]
G. P. Lepage and B. A. Thacker,

Nucl. Phys. B Proc. Suppl. 4 (1988), 199.

[20]
N. Isgur and M. B. Wise,

Phys. Lett. B 232 (1989), 113-117.

[21]
A. F. Falk and B. Grinstein,

Phys. Lett. B 249 (1990), 314-320.

[22]
T. Mannel and Z. Ryzak,

Phys. Lett. B 247 (1990), 412-417.

[23]
T. D. Cohen and J. Milana,

Phys. Lett. B 306 (1993), 134-138.

[24]
N. Isgur and M. B. Wise,

Phys. Rev. Lett. 66 (1991), 1130-1133.

[25]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 221 (1989), 422-426.

[26]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 232 (1989), 398-404.

[27]
P. Avery et al. [CLEO],

Phys. Rev. D 41 (1990), 774.

[28]
M. Lu, M. B. Wise and N. Isgur,

Phys. Rev. D 45 (1992), 1553-1556.

[29]
J. P. Lees et al. [BaBar],

Phys. Rev. D 83 (2011), 072003.

[30]
V. Balagura et al. [Belle],

Phys. Rev. D 77 (2008), 032001.

[31]
B. Aubert et al. [BaBar],

Phys. Rev. D 78 (2008), 034008.

[32]
G. Abbiendi et al. [OPAL],

Phys. Lett. B 444 (1998), 539-554.

[33]
L. K. Li et al. [Belle],

Sci. Bull. 68 (2023), 583-592.

[34]
M. A. Ivanov, J. G. Korner, V. E. Lyubovitskij and A. G. Rusetsky,

Phys. Rev. D 60 (1999), 094002.

[35]
M. A. Ivanov, J. G. Korner and V. E. Lyubovitskij,

Phys. Lett. B 448 (1999), 143-151.

[36]
A. F. Falk and T. Mehen,

Phys. Rev. D 53 (1996), 231-240.

[37]
P. del Amo Sanchez et al. [BaBar],

Phys. Rev. D 82 (2010), 111101.

[29]
J. P. Lees et al. [BaBar],

Phys. Rev. D 83 (2011), 072003.

[38]
R. Aaij et al. [LHCb],

JHEP 09 (2013), 145.

[39]
F. James and M. Roos,

Comput. Phys. Commun. 10 (1975), 343-367

[40]
P. L. Cho and H. Georgi,

Phys. Lett. B 296 (1992), 408-414
[erratum: Phys. Lett. B 300 (1993), 410].

[41]
W. A. Bardeen, E. J. Eichten and C. T. Hill,

Phys. Rev. D 68 (2003), 054024.

[42]
P. L. Cho and M. B. Wise,

Phys. Rev. D 49 (1994), 6228-6231.

[43]
D. Buskulic et al. [ALEPH],

Phys. Lett. B 365 (1996), 437-447.

[44]
G. Bewick, S. Ferrario Ravasio, P. Richardson and M. H. Seymour,

JHEP 04 (2020), 019.

[45]
G. Bewick, S. Ferrario Ravasio, P. Richardson and M. H. Seymour,

JHEP 01 (2022), 026.

[46]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 221 (1989), 422-426.

[47]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 232 (1989), 398-404.

[48]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 230 (1989), 162.

[49]
H. Albrecht et al. [ARGUS],

Z. Phys. C 52 (1991), 353-360.

[50]
J. P. Alexander et al. [CLEO],

Phys. Lett. B 303 (1993), 377-384.

[51]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 297 (1992), 425-431.

[52]
H. Albrecht et al. [ARGUS],

Z. Phys. C 69 (1996), 405-408.

[53]
P. Avery et al. [CLEO],

Phys. Lett. B 331 (1994), 236-244
[erratum: Phys. Lett. B 342 (1995), 453-453].

[54]
T. Bergfeld et al. [CLEO],

Phys. Lett. B 340 (1994), 194-204.

[55]
M. Artuso et al. [CLEO],

Phys. Rev. D 70 (2004), 112001.

[56]
R. Seuster et al. [Belle],

Phys. Rev. D 73 (2006), 032002.

[57]
B. Aubert et al. [BaBar],

Phys. Rev. D 65 (2002), 051101.

[58]
M. Bishai et al. [CLEO],

Phys. Lett. B 482 (2000), 203-209.

[59]
B. Aubert et al. [BaBar],

Phys. Rev. D 73 (2006), 091101.

[60]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 247 (1990), 121-126.

[61]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 402 (1997), 207-212.

[62]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 317 (1993), 227-232.

[63]
H. Albrecht et al. [ARGUS],

Phys. Lett. B 211 (1988), 489-492.

[64]
J. P. Alexander et al. [CLEO],

Phys. Rev. Lett. 83 (1999), 3390-3393.

[65]
B. Aubert et al. [BaBar],

Phys. Rev. Lett. 95 (2005), 142003.

[66]
B. Aubert et al. [BaBar],

Phys. Rev. D 75 (2007), 012003.

[67]
B. Aubert et al. [BaBar],

[arXiv:hep-ex/0607086 [hep-ex]].

[68]
B. Aubert et al. [BaBar],

Phys. Rev. Lett. 97 (2006), 232001.

[69]
B. Aubert et al. [BaBar],

Phys. Rev. Lett. 99 (2007), 062001.

[70]
P. Avery et al. [CLEO],

Phys. Rev. D 43 (1991), 3599-3610.

[71]
P. Avery et al. [CLEO],

Phys. Rev. Lett. 75 (1995), 4364-4368.

[72]
T. J. V. Bowcock et al. [CLEO],

Phys. Rev. Lett. 63 (1989), 2189-2192.

[73]
G. Brandenburg et al. [CLEO],

Phys. Rev. Lett. 78 (1997), 2304-2308.

[74]
K. W. Edwards et al. [CLEO],

Phys. Rev. Lett. 74 (1995), 2650-2653.

[75]
L. Gibbons et al. [CLEO],

Phys. Rev. Lett. 77 (1996), 810-813.

[76]
C. P. Jessop et al. [CLEO],

Phys. Rev. D 58 (1998), 032003.

[77]
R. Mizuk et al. [Belle],

Phys. Rev. D 70 (2004), 031104.

[78]
M. Niiyama et al. [Belle],

Phys. Rev. D 95 (2017) no.3, 032003.

[79]
B. Aubert et al. [BaBar],

Phys. Rev. Lett. 87 (2001), 162002.

[80]
K. Abe et al. [Belle],

Phys. Rev. Lett. 88 (2002), 052001.