\section{Introduction}

The baryon asymmetry of the Universe (BAU) is one of the most mysterious unsolved problems in the Standard Model (SM) of particle physics. Leptogenesis is a very attractive explanation~\cite{Fukugita:1986hr,Luty:1992un,Davidson:2008bu}, as it links the BAU to the origin of neutrino masses. In that mechanism, the $CP$ violating decay of the heavy right-handed neutrinos (RHN) to the SM leptons generates a lepton asymmetry, which is then converted to the BAU via the electroweak (EW) sphaleron process. The same interaction also accounts for the tiny neutrino masses~\cite{Davis:1994jw, Fukuda:1998mi,Eguchi:2002dm} via the Type-I seesaw mechanism~\cite{Minkowski:1977sc}. In the conventional thermal leptogenesis formalism, the $CP$ violating effects are related to the RHN mass $M_N$, and the explanation of the BAU requires $M_N\gtrsim10^9$ GeV~\cite{Davidson:2002qv}, making it very challenging to test the mechanism experimentally. However, the constraints on $M_N$ can be relaxed if at least two of the RHNs are highly degenerate, so that the $CP$ asymmetry is resonantly enhanced, and hence even $\mathcal{O}({\rm TeV})$ RHNs can generate the observed BAU~\cite{Flanz:1996fb,Pilaftsis:1997jf,Pilaftsis:2003gt,Iso:2010mv,Dev:2017wwc}. In particular, if such a resonant leptogenesis mechanism is embedded into a gauge theory in which the RHNs exist naturally for anomaly cancellation, then the $\mathcal{O}({\rm TeV})$ leptogenesis is testable at the colliders via searches for RHNs or new gauge/scalar bosons as well as the leptonic charge asymmetries~\cite{Blanchet:2009bu,Okada:2012fs,Heeck:2016oda,Dev:2017xry, BhupalDev:2014hro, Chun:2017spz, BhupalDev:2015khe, Dev:2019rxh}.

In this article, we perform a comprehensive study of the collider phenomenology of the resonant leptogenesis mechanism in the gauged $U(1)_{B-L}$ extended SM (i.e. the so-called $B-L$ model~\cite{Davidson:1978pm,Marshak:1979fm,Mohapatra:1980qe,Davidson:1987mh}), focusing on the interplay between leptogenesis and the $Z'$ gauge boson of the $U(1)_{B-L}$ group, since the $Z'$-mediated scattering process greatly impacts the generated BAU. In particular, we use the $Z'\to NN\to\ell^\pm\ell^\pm+{\rm jets}$ channel to probe the $Z'$, $N$ particles as well as the $CP$ violation. Compared with previous studies on similar topics~\cite{Blanchet:2009bu,Okada:2012fs,Heeck:2016oda,Dev:2017xry, BhupalDev:2014hro, Chun:2017spz, BhupalDev:2015khe, Dev:2019rxh}, our work includes not only the newest constraints from the LHC and the corresponding projected reach at the future HL-LHC, but also the first projections at the future multi-TeV muon colliders. The study on the physics potential of muon colliders started around three decades ago~\cite{Barger:1995hr,Barger:1996jm}, and receives a renewed interest recently~\cite{Han:2012rb,Chakrabarty:2014pja,Ruhdorfer:2019utl,DiLuzio:2018jwd,Delahaye:2019omf,Long:2020wfp,Buttazzo:2018qqp,Costantini:2020stv,Han:2020uid,Capdevilla:2020qel,Han:2020uak,Han:2020pif,Bartosik:2020xwr,Chiesa:2020awd,Yin:2020afe,Buttazzo:2020eyl,Lu:2020dkx,Huang:2021nkl,Liu:2021jyc,Cheung:2021iev,Han:2021udl,Capdevilla:2021fmj,Huang:2021biu,Han:2021kes,AlAli:2021let,Asadi:2021gah,Franceschini:2021aqd,Buarque:2021dji,Han:2021lnp,Chiesa:2021qpr}. Due to the high energy and precision measurement environment, the future multi-TeV muon colliders offer us the opportunity to probe both SM and beyond the SM physics very accurately.

This paper is organized as follows. The $B-L$ model and the corresponding (resonant) leptogenesis is described in Section~\ref{sec:blz}, where the relation between BAU and the size of $CP$ asymmetry is derived. Then we study the collider phenomenology in Section~\ref{sec:mc}, including the reach for $Z'$ boson and the $CP$ asymmetry (via RHNs) at the HL-LHC and muon colliders. Finally, the conclusion is given in Section~\ref{sec:con}.

\section{Resonant leptogenesis in the $B-L$ model}\label{sec:blz}

\subsection{The model}

In the $B-L$ model, three generations of RHNs (with $B-L=-1$) are needed naturally for gauge anomaly cancellation. Besides, the model also contains an extra gauge boson $Z'$ and a complex scalar $\Phi=(\phi+i\eta)/\sqrt2$ with $B-L=2$ that breaks the $U(1)_{B-L}$ spontaneously. The relevant Lagrangian reads
\begin{equation}\label{LB-L}\begin{split}
\mathcal{L}_{B-L}=&~\sum_{i}\bar\nu_R^ii\slashed{D}\nu_R^i-\frac12\sum_{i,j}\left(\lambda_N^{ij}\bar\nu_R^{i,c}\Phi\nu_R^j+{\rm h.c.}\right)
-\sum_{i,j}\left(\lambda_D^{ij}\bar\ell_L^i\tilde H\nu_R^j+{\rm h.c.}\right)\\
&~+D_\mu\Phi^\dagger D^\mu\Phi-\lambda_\phi\left(|\Phi|^2-\frac{v_\phi^2}{2}\right)^2-\frac14Z'_{\mu\nu}Z'^{\mu\nu},
\end{split}\end{equation}
with $i$, $j$ being the family indices, $D_\mu=\partial_\mu-ig_{B-L}XZ_\mu'$ the covariant derivative ($X$ is the $B-L$ quantum number), $\ell_L$ the SM left-handed lepton doublet and $H$ the SM Higgs doublet. Without loss of generality, we assume $\lambda_N^{ij}={\rm diag}\{\lambda_{N_1},\lambda_{N_2},\lambda_{N_3}\}$, and define the four-component Majorana RHNs as $N_i=\nu_R^i+(\nu_R^i)^c$. The minimum of the scalar potential is at $\ave{|\Phi|}=v_\phi/\sqrt{2}$, breaking the $U(1)_{B-L}$ symmetry and providing masses for the particles
\begin{equation}
M_{Z'}=2g_{B-L}v_\phi,\quad M_{N_i}=\lambda_{N_i}\frac{v_\phi}{\sqrt2},\quad M_\phi=\sqrt{2\lambda_\phi}v_\phi,
\end{equation}
and the imaginary part $\eta$ of $\Phi$ is absorbed to be the longitudinal mode of $Z'$. We are interested in the parameter region $v_\phi\sim\mathcal{O}({\rm TeV})$.

The usage of Yukawa interactions between $\ell$ and $N$ is two-fold. On one hand, they account for the tiny left-handed neutrino mass via the Type-I seesaw
\begin{equation}
m_{\nu_L}\sim\frac{\lambda_D^2v^2}{M_N}\sim0.06~{\rm eV}\times\left(\frac{\lambda_D}{10^{-6}}\right)^2\left(\frac{1~{\rm TeV}}{M_N}\right),
\end{equation}
while on the other hand they trigger the RHN decay $N\to\ell H/\bar\ell H^*$, which is the crucial process in thermal leptogenesis. Due to the $CP$ violation phase in the Yukawa couplings, the widths of $N\to\ell H$ and $N\to\bar\ell H^*$ are different, and hence a $CP$ asymmetry can be defined as~\cite{Flanz:1996fb,Pilaftsis:1997jf,Pilaftsis:2003gt,Iso:2010mv}
\begin{equation}\label{epsilon}
\epsilon_i=\frac{\sum_j\Gamma_{N_i\to\ell_jH}-\Gamma_{N_i\to\bar\ell_jH^*}}{\sum_j\Gamma_{N_i\to\ell_jH}+\Gamma_{N_i\to\bar\ell_jH^*}}=-\sum_{j\neq i}\frac{M_{N_i}\Gamma_{N_j}}{M_{N_j}^2}\left(\frac{V_{ij}}{2}+S_{ij}\right)\frac{{\rm Im}(\lambda_D\lambda_D^\dagger)^2_{ij}}{(\lambda_D\lambda_D^\dagger)_{ii}(\lambda_D\lambda_D^\dagger)_{jj}},
\end{equation}
where
\begin{equation}
V_{ij}=2\frac{M_{N_j}^2}{M_{N_i}^2}\left[\left(1+\frac{M_{N_j}^2}{M_{N_i}^2}\right)\ln\left(1+\frac{M_{N_j}^2}{M_{N_i}^2}\right)-1\right],\quad S_{ij}=\frac{M_{N_j}^2(M_{N_j}^2-M_{N_i}^2)}{(M_{N_j}^2-M_{N_i}^2)^2+M_{N_i}^2\Gamma_{N_j}^2},
\end{equation}
are respectively the vertex correction and RHN self-energy correction to the decay process, and the tree level width is
\begin{equation}
\Gamma_{N_j}=\frac{M_{N_j}}{8\pi}(\lambda_D\lambda_D^\dagger)_{jj}.
\end{equation}
If there is a mass hierarchy among the three RHNs, i.e. $M_{N_1}\ll M_{N_2}\ll M_{N_3}$, $V_{ij}$ and $S_{ij}$ are comparable, and $\epsilon_i$ is proportional to the lightest RHN mass but typically $\gtrsim10^{-6}$. Therefore, a sizable BAU requires a RHN with mass $\gtrsim10^9$ GeV~\cite{Davidson:2002qv}. However, if at least two RHNs are highly degenerate that $|M_{N_j}^2-M_{N_i}^2|\sim M_{N_i}\Gamma_{N_j}$, then $S_{ij}\sim M_{N_j}/\Gamma_{N_j}\gg1$, and $\epsilon_i$ can reach $\mathcal{O}(1)$~\cite{Flanz:1996fb,Pilaftsis:1997jf,Pilaftsis:2003gt,Iso:2010mv,Dev:2017wwc}. In this case, even $\mathcal{O}(\rm TeV)$ RHNs can generate a successful BAU. This is the scenario under consideration in this article.

\subsection{Thermal leptogenesis}

We assume first two flavors of RHNs are near-degenerate, and the third flavor RHN is much heavier. In that case, if there is a hierarchy between the decay widths of the two near-degenerate RHNs, then the generated BAU will be dominated by the RHN that with larger decay width~\cite{Iso:2010mv}. For simplicity we will consider such a scenario and hence the one-flavor analysis is sufficient. We then denote $\epsilon_1$ as $\epsilon$, and $M_{N_1}\approx M_{N_2}$ as $M_N$ from now on. In the radiation dominated era, the energy and entropy densities of the Universe are respectively
\begin{equation}
\rho=\frac{\pi^2}{30}g_*T^4,\quad s=\frac{2\pi^2}{45}g_*T^3,
\end{equation}
where $g_*=106.75$ is the number of relativistic degrees of freedom. The Hubble constant is then derived by the first Friedmann equation $H^2=(8\pi/3M_{\rm Pl}^2)\rho$, with $M_{\rm Pl}=1.22\times10^{19}$ GeV the Planck scale. Defining the dimensionless parameter $z\equiv M_N/T$, the Boltzmann equations for the RHN and net $B-L$ number yield in the thermal bath read~\cite{Iso:2010mv,Plumacher:1996kc}
\begin{eqnarray}\label{leptogenesis}
\frac{s_NH_N}{z^4}\frac{dY_N}{dz}&=&-\left(\frac{Y_N}{Y_N^{\rm eq}}-1\right)(\gamma_D+2\gamma_{h,s}+4\gamma_{h,t})-\left(\frac{Y_N^2}{(Y_N^{\rm eq})^2}-1\right)2\gamma_{Z'},\\
\frac{s_NH_N}{z^4}\frac{dY_{B-L}}{dz}&=&-\left[\frac{Y_{B-L}}{2Y_\ell^{\rm eq}}-\epsilon\left(1-\frac{Y_N}{Y_N^{\rm eq}}\right)\right]\gamma_D-\frac{Y_{B-L}}{Y_\ell^{\rm eq}}\left[2(\gamma_{N,s}+\gamma_{N,t}+\gamma_{h,t})+\frac{Y_N}{Y_N^{\rm eq}}\gamma_{h,s}\right],\nonumber
\end{eqnarray}
where $s_N$ and $H_N$ are the entropy density and Hubble constant at $z=1$, respectively, and the abundances are defined as number density to entropy density ratios (e.g. $Y_N=n_N/s$), with
\begin{equation}
Y_N^{\rm eq}=\frac{45z^2}{2\pi^4g_*}K_2(z),\quad Y_\ell^{\rm eq}=\frac{3\zeta(3)}{2\pi^2}\frac{M_N^3}{s_N},
\end{equation}
the equilibrium abundances, and $K_i(z)$ is the modified Bessel function of the $i$th kind.

The reaction rates in \Eq{leptogenesis} are defined as~\cite{Plumacher:1996kc},
\begin{equation}
\gamma_D=\frac{M_N^3}{\pi^2}\Gamma_N\frac{K_1(z)}{z},
\end{equation}
for the $N\to\ell H/\bar\ell H^*$ decay, and
\begin{equation}
\gamma_{ab\to cd}\equiv\ave{\sigma_{ab\to cd}v}n_a^{\rm eq}n_b^{\rm eq}=\frac{M_N}{64\pi^4z}\int_{s_{\rm min}}^\infty ds\,\hat\sigma_{ab\to cd}(s)\sqrt{s}K_1\left(\frac{\sqrt{s}}{M_N}z\right),
\end{equation}
for the $2\to2$ scattering, where $s_{\rm min}=\max\{(M_a+M_b)^2,~(M_c+M_d)^2\}$, and the dimensionless reduced scattering cross section is
\begin{equation}
\hat\sigma_{ab\to cd}(s)\equiv 2\,\sigma_{ab\to cd}(s)\cdot s\cdot\left[1-2\left(\frac{M_a^2}{s}+\frac{M_b^2}{s}\right)+\left(\frac{M_a^2}{s}-\frac{M_b^2}{s}\right)^2\right].
\end{equation}
The correspondence between the reaction rates in \Eq{leptogenesis} and the $2\to2$ processes are listed in Table~\ref{tab:reaction_rates}, except the case of $\gamma_{Z'}$, which corresponds to the scattering $NN\to Z'\to f\bar f$ ($f$ denotes the SM fermions) and is highlighted below
\begin{equation}
\hat\sigma_{Z'}(s)=\frac{13g_{B-L}^4}{6\pi}\frac{\sqrt{x(x-4)^3}}{(x-M_{Z'}^2/M_N^2)^2+M_{Z'}^2\Gamma_{Z'}^2/M_N^4},
\end{equation}
where
\begin{equation}
\label{eq:gammaz}
\Gamma_{Z'}=\frac{g_{B-L}^2}{24\pi}M_{Z'}\left[13+2\left(1-\frac{4M_N^2}{M_{Z'}^2}\right)^{3/2}\theta(M_{Z'}-2M_N)\right],
\end{equation}
with $x=s/M_N^2$ and $\theta$ the Heaviside step function. We also assume $M_{Z'}$, $M_\phi\gtrsim M_N$ so that scatterings $NN\to Z'Z'$ or $\phi\phi$ are suppressed.\footnote{The impact of $NN\to Z'Z'$ and $NN\to\phi\phi$ can be found in Ref.~\cite{Heeck:2016oda}.}

\begin{table}\small\renewcommand\arraystretch{1.8}\centering
\begin{tabular}{|c|c|c|c|c|c|}\hline
Reaction rates & $2\to2$ scattering & $\hat\sigma(s)$ \\ \hline
$\gamma_{h,s}$ & $N_1\ell\to\bar tq$ & $\frac{3y_t^2\lambda_D^2}{4\pi}\left(\frac{x-1}{x}\right)^2$ \\ \hline
$\gamma_{h,t}$ & $N_1q\to t\ell$ or $N_1\bar t\to\bar q\ell$ & $\frac{3y_t^2\lambda_D^2}{4\pi}\left(\frac{x-1}{x}+\frac1x\ln\frac{x-1+M_h^2/M_N^2}{M_h^2/M_N^2}\right)$ \\ \hline
$\gamma_{N,s}$ & $\ell H\to\bar\ell H^*$ & $\frac{\lambda_D^4}{2\pi}\left[1+\frac{2}{D_1(x)}+\frac{x}{2D_1^2(x)}-\left(1+\frac{2(x+1)}{D_1(x)}\right)\frac{\ln(1+x)}{x}\right]$ \\ \hline
$\gamma_{N,t}$ & $\ell\ell\to H^*H^*$ & $\frac{\lambda_D^4}{2\pi}\left(\frac{x}{2(x+1)}+\frac{\ln(1+x)}{x+2}\right)$ \\ \hline
\end{tabular}
\caption{The reaction rates and the reduced cross sections taken from Refs.~\cite{Iso:2010mv,Plumacher:1996kc}, where $x=s/M_N^2$ and $D_1(x)=x-1+(\Gamma_N^2/M_N^2)/(x-1)$.}\label{tab:reaction_rates}
\end{table}

For a $\mathcal{O}(\rm TeV)$ leptogenesis, $\lambda_D\sim10^{-6}$, making the $\gamma$'s in Table~\ref{tab:reaction_rates} rather small, as they are proportional to $\lambda_D^2$ or $\lambda_D^4$. Only the reaction rate $\gamma_{Z'}$ can be sizable since it is $\propto g_{B-L}^4$. Therefore we can omit the reaction rates in \Eq{leptogenesis} except $\gamma_D$ and $\gamma_{Z'}$, and simplify \Eq{leptogenesis} into a single equation,
\begin{equation}\label{leptogenesis_sim}
\frac{dY_{B-L}}{dz}+\left(\frac{z^4}{s_NH_N}\frac{\gamma_D}{2Y_\ell^{\rm eq}}\right)Y_{B-L}\approx\epsilon \frac{dY_N^{\rm eq}}{dz}\frac{\gamma_D}{\gamma_D+4\gamma_{Z'}},
\end{equation}
where we have approximated $Y_N/Y_N^{\rm eq}+1\approx2$ and $dY_N/dz\approx dY_N^{\rm eq}/dz$, since $N$ is not far away from equilibrium. \Eq{leptogenesis_sim} can be solved analytically as
\begin{equation}\label{eq:kappa}
Y_{B-L}(z)\approx\epsilon\int_{z_{\rm in}}^zdz'\frac{dY_N^{\rm eq}}{dz'}\frac{\gamma_D(z')}{\gamma_D(z')+4\gamma_{Z'}(z')}\exp\left\{-\int_{z'}^zdz''\frac{z''^4}{s_NH_N}\frac{\gamma_D(z'')}{2Y_\ell^{\rm eq}}\right\},
\end{equation}
where we adopt $z_{\rm in}=1$ as the lower limit of the integral. It is very clear in above equation that $\epsilon$ and $\gamma_D$ generate the lepton asymmetry, while $\gamma_{Z'}$ tends to washout this asymmetry because it tends to push $N$ back to the equilibrium. The final generated baryon asymmetry is
\begin{equation}
Y_B=\frac{28}{79}\,Y_{B-L}(z_{\rm sph}),
\end{equation}
where $z_{\rm sph}=M_N/T_{\rm sph}$, with $T_{\rm sph}\approx130$ GeV the decoupled temperature of the EW sphaleron~\cite{Burnier:2005hp}. We have checked that \Eq{eq:kappa} matches the numerical solution of the complete equation set \Eq{leptogenesis} very well.

\begin{figure}
\includegraphics[scale=0.35]{Figures/asy_blz_717.png}\quad
\includegraphics[scale=0.35]{Figures/mzp_gBL_limit.png}
\caption{Left: The required $CP$ asymmetry $\epsilon$ to explain the observed BAU in the mass range $T_{\rm sph.}<M_N<M_{Z'}$, for a fixed $g_{B-L}=0.8$. Only the $\epsilon$ contours for $M_{Z'}>M_N$ are shown, because for $M_{Z'}<M_N$ new scattering channels open and \Eq{eq:kappa} needs to be modified. Right: The upper limits on $g_{B-L}$ as a function of $M_{Z^\prime}$ from the LEP~\cite{Abdallah:2005ph}, ATLAS dilepton searches~\cite{ATLAS:2019erb, Chiang:2019ajm}, and projected reach from the 3 TeV muon colliders~\cite{Huang:2021nkl}, and the 10 (30) TeV muon colliders via rescaling.}
\label{fig:blz}
\end{figure}

Given the value of $g_{B-L}$, one is able to derive the $CP$ asymmetry $\epsilon$ as a function of $(M_{Z'},M_N)$ via \Eq{eq:kappa} by the observed BAU $Y_B^{\rm obs.}\approx10^{-10}$~\cite{Ade:2015xua, Planck:2018vyg}. This is shown in the left panel of Fig.~\ref{fig:blz}, where $g_{B-L}=0.8$ is fixed. Near the line $M_{Z'}=2M_N$, the washout process $NN\to Z'\to f\bar f$ is resonantly enhanced and hence a large $\epsilon$ is needed to realize $Y_B^{\rm obs.}$. The region with $\epsilon>1$ is forbidden in the leptogenesis mechanism. We can see that there is plenty of parameter space allowed by leptogenesis for $\mathcal{O}({\rm TeV})$ $Z'$ and $N$. Since this mass region is accessible at current or near future colliders~\cite{Ilten:2018crw, Cepeda:2019klc, Bagnaschi:2019djj, Deppisch:2019ldi, Deppisch:2019kvs}, some of the parameter space is already excluded by the LEP~\cite{Abdallah:2005ph} and LHC~\cite{ATLAS:2019erb, Chiang:2019ajm} searches for $Z'\to\ell^+\ell^-/jj$, as plotted in the shaded region in the right panel of Fig.~\ref{fig:blz}. Also plotted in the figure is the projected reach for $\mu^+\mu^-\to Z'^*\to\ell^+\ell^-$ and $\mu^+\mu^-\to Z'(\to\ell^+\ell^-/\nu_\ell\bar\nu_\ell)\gamma$ at the multi-TeV muon colliders taken/rescaled from Ref.~\cite{Huang:2021nkl}.
We are mostly interested in the parameter space that is allowed by current data but can be probed at the HL-LHC and future muon colliders; for this sake we fix $g_{B-L}=0.8$ and focus on $Z^\prime$ with $6~{\rm TeV} < M_{Z^\prime} <  30$ TeV to perform the collider phenomenology study.

\section{Collider phenomenology}\label{sec:mc}

\subsection{Same-sign dilepton final state}

Since the SM fermions are charged under the $U(1)_{B-L}$ group, $Z'$ can be produced at the LHC via quark fusion $q\bar q\to Z'$ or at the muon colliders via $\mu^+\mu^-\to Z'\gamma/Z'Z$ (the so-called radiative return) and $\mu^+\mu^-\to Z'$. Since the collision energy $\sqrt{s}$ of a muon collider is fixed, the $Z'\gamma/Z'Z$ channel is suited for probing $M_{Z'}<\sqrt{s}$, while the $\mu^+\mu^-\to Z'$ is more appropriate to probe the off-shell region $M_{Z'}>\sqrt{s}$. As for the decay of $Z'$, we focus on the channel
\begin{equation}\label{ssdl}
Z'\to NN\to\ell^\pm\ell^\pm+{\rm jets},
\end{equation}
channel as it is directly related to the essential ingredients of the leptogenesis mechanism: the new particles $Z'$, $N$ and the $CP$ asymmetry. By reconstructing the invariant masses of the decay products, we can find  clues for the $Z'$ and $N$ resonances; while by counting the event number difference between the $\ell^+\ell^+$ and $\ell^-\ell^-$ final states, we can probe the $CP$ violation.

\begin{figure}\centering
	\includegraphics[scale=0.35]{Figures/csnn.png}
	\caption{The production rates of $Z'\to NN$ at the 13 TeV LHC and the 10 (30) TeV muon colliders, with $M_N = M_{Z^\prime}/3$ fixed.}
	\label{fig:cszp}
\end{figure}

The cross sections of $Z'\to NN$ for various production channels are plotted in Fig.~\ref{fig:cszp} at the 13 TeV LHC, 10 and 30 TeV muon colliders where $M_{N} = M_{Z'}/3$ is fixed.\footnote{We adopt the {\tt FeynRules}~\cite{Alloul:2013bka} model file from Refs.~\cite{Amrith:2018yfb, Deppisch:2018eth}, which is also publicly available  in the {\tt FeynRules} model database, \url{https://feynrules.irmp.ucl.ac.be/wiki/B-L-SM}. The model file is then interfaced with the {\tt MadGraph5aMC@NLO-v2.8.2}~\cite{Alwall:2014hca} package for a parton-level simulation.} We can see the LHC cross section drops rapidly as $M_{Z'}$ increases, and it is only $\sigma\sim10^{-4}$ fb for $M_{Z'}=12$ TeV, yielding less than 1 event even at the HL-LHC. In contrast, the $\mu^+\mu^-\to Z'\gamma/Z'Z$ cross section increases when $M_{Z'}$ is close to the collision energy of the muon collider; for $M_{Z'}>6$ TeV, $\sigma\gtrsim10~(1)$ fb at a 10 (30) TeV muon collider, providing $\gtrsim10^5~(9\times10^4)$ events at an integrated luminosity of $10~(90)~{\rm ab}^{-1}$. The 10 TeV muon collider can probe $M_{Z'}>10$ TeV via the off-shell production of $\mu^+\mu^-\to Z'^*\to NN$. In this case, the cross section drops when $M_{Z'}$ increases, showing a tendency similar to the LHC but with much higher rates (for example $\sigma\sim10$ fb for $M_{Z'}=12$ TeV).

The RHN decays to a lepton and a vector/Higgs boson, and due to the Goldstone equivalence theorem the branching ratios are approximately
\begin{equation}
\frac12{\rm Br}(\ell^+W^-)\approx{\rm Br}(\bar\nu_\ell Z)\approx{\rm Br}(\bar\nu_\ell h),\quad \frac12{\rm Br}(\ell^-W^+)\approx{\rm Br}(\nu_\ell Z)\approx{\rm Br}(\nu_\ell h).
\end{equation}
The difference between ${\rm Br}(\ell^+W^-)$ and ${\rm Br}(\ell^-W^+)$ connects to the $CP$ violation $\epsilon$, as defined in \Eq{epsilon}. For simplicity, in the following we assume the $N$ decays exclusively to the first generation of leptons, and focus on the $NN\to e^\pm e^\pm W^\mp W^\mp$ final state, which is clean due to its same-sign dilepton feature. As we are considering $M_{Z'}\gtrsim6$ TeV, the $W^\pm$ bosons from cascade decay are typically boosted, and hence can be treated as boosted fat jets.

To trigger the signal events, we require the final state to have exactly two electrons within
\begin{equation}\begin{split}
&p_T^e>100~{\rm GeV},\quad |\eta_e|<2.5,\quad {\rm (for~LHC)}\\
&p_T^e>30~{\rm GeV},\quad |\eta_e|<2.43;\quad {\rm (for~muon~colliders)}\\
\end{split}\end{equation}
and two jets within
\begin{equation}\begin{split}
&p_T^W>500~{\rm GeV},\quad |\eta_W|<2,\quad {\rm (for~LHC)}\\
&p_T^W>500~{\rm GeV},\quad |\eta_W|<2.43.\quad {\rm (for~muon~colliders)}\\
\end{split}\end{equation}
The electrons are further required to be same-sign, and the jets are required to be $W$-tagged. At the LHC, the main backgrounds come from the dilepton decay of $t\bar t$ (with one lepton's charge misidentified and two QCD jets mistagged as $W$) and $W^\pm W^\pm jj$ (with jets mistagged)~\cite{Blanchet:2009bu}. While at the muon colliders, the $t\bar t$ background is subdominant; instead, the main backgrounds are $\mu^+ \mu^- \rightarrow e^+ e^- W^+ W^-$ and $\mu^+ \mu^- \rightarrow e^+ e^- W^+ W^- \gamma/Z$ (with lepton charge misidentified), $\mu^+ \mu^- \rightarrow W^+ W^- j j $ (with jets mistagged). The backgrounds from charge misidentification are significantly suppressed by the misidentification rate, which we adopt as 0.1\% based on the LHC detector performance~\cite{ATLAS:2019jvq}. Note that this is a conservative estimate as the charge identification efficiency is expected to be improved at the muon colliders. In our parton-level study, the $W$-jet is simulated by a parton-level $W$ boson multiplying by the hadronic decay branching ratio 67.4\%~\cite{ParticleDataGroup:2020ssz} and the tagging efficiency 60\%. The backgrounds from mistagged QCD jets are then suppressed by the mistag rate, which we adopt as 5\%. The $t\bar t$ background is further suppressed by requiring the final state objects' invariant mass to be above 6 TeV, while the same cut almost does not affect the signal rate. The background rates are listed in Table~\ref{tab:back}, where the signal cross sections for $M_{Z'}=8$ TeV, $M_N=500$ GeV, and $g_{B-L}=0.8$ are also given. We can see that at the LHC the backgrounds are under control, while at the muon colliders the backgrounds are negligible after the selection cuts.

\begin{table}
	\small\centering\renewcommand\arraystretch{1.2}
	\begin{tabular}{c|ccc}
		\hline\hline
		LHC  & Trigger cut [fb] &   Same-sign lepton [fb] & $W$-jet [fb]
		\\ \hline
		Signal & $\sim10^{-3}$  & $\sim$10$^{-3}$ & $\sim$10$^{-4}$   \\
		$t\bar t$ & $\sim 10^{-4}$ ($*$)  & $\lesssim10^{-7}$ &  $\lesssim10^{-10}$  \\
		$W^\pm W^\pm jj$  & $\lesssim10^{-2}$  & $\lesssim10^{-4}$ & $\lesssim10^{-7}$  \\ \hline\hline
		10 TeV muon collider  & Trigger cut [fb]   & Same-sign lepton [fb] & $W$-jet [fb]
		\\ \hline
		Signal & $\sim 1$  & $\sim 1$ & $\sim10^{-1}$  \\
		$\mu^+ \mu^- \rightarrow e^+ e^- W^+ W^-$ & $\sim10^{-2}$  & $\sim10^{-5}$ &  $\sim10^{-6}$ \\
		$\mu^+ \mu^- \rightarrow e^+ e^- W^+ W^- \gamma/Z$ & $\sim10^{-2}$  & $\sim10^{-5}$ & $\sim10^{-6}$  \\
		$\mu^+ \mu^- \rightarrow W^+ W^- j j $ & $\sim 10^{-1}$  & $\sim10^{-6}$ & $\sim10^{-9}$ \\

		\hline\hline
	\end{tabular}
	\caption{The cross sections of the signal and main backgrounds at the 13 TeV LHC and 10 TeV muon collider after the trigger cuts, the same-sign dilepton and $W$-jet requirements. ($*$): The decay products of $t \bar{t}$ are further required to have invariant masses $>6$ TeV. The signal process is generated at $M_{Z'}=8$ TeV, $M_N=500$ GeV, and $g_{B-L}=0.8$.}
	\label{tab:back}
\end{table}

\begin{figure}
	\includegraphics[scale=0.23]{Figures/wpt.png}
	\includegraphics[scale=0.23]{Figures/hnlmass.png}
	\includegraphics[scale=0.233]{Figures/zpmass.png}
	\caption{$p^W_T$ and the invariant mass distribution of the signal and $\mu^+ \mu^- \rightarrow e^+ e^- W^+ W^- \gamma$ at the 10 TeV muon collider for $(M_{Z^\prime} =$ 8 TeV, $M_N = $ 500 GeV).  }
	\label{fig:inv}
\end{figure}

At the muon colliders, we can further remove the backgrounds also by putting cuts on the invariant mass of the final state particles. This is shown in Fig.~\ref{fig:inv}, where $p^W_T$ and the invariant masses of the signal and $\mu^+ \mu^- \rightarrow e^+ e^- W^+ W^- \gamma$ at the 10 TeV muon collider are illustrated for one benchmark point with $M_{Z^\prime} =$ 8 TeV and $M_N = $ 500 GeV. The jet four-momentum is smeared according to a jet energy resolution of $\Delta E/E = 10 \%$. For the reconstruction of $N$, we pair the two leptons and $W$-jets by minimizing $\chi_j^2=(M_{l_1^\pm W_1} - M_N)^2+(M_{l_2^\pm W_2} - M_N)^2$. From the figure, it is clear to see that the signal and background have very different kinematical distributions especially when reconstructing the $N$; a cut on the $\ell^\pm W$ invariant mass can significantly remove the backgrounds, even in case of $M_{Z'}>\sqrt{s}$.

\subsection{Testing Leptogenesis}

\begin{figure}\centering
	\includegraphics[scale=0.35]{Figures/senmc.png}
	\caption{The regions where $\epsilon>\epsilon_{\rm min}$ that the resonant leptogenesis can be tested at the HL-LHC and the 10~(30) TeV muon collider with an integrated luminosity of 10~(90) ab$^{-1}$. $g_{B-L} =0.8$ is fixed. The dark blue region labeled as ``no leptogenesis'' is for $\epsilon>1$ or $M_N<T_{\rm sph}$.}
	\label{fig:asy_zp}
\end{figure}

We have shown in the previous subsection that the backgrounds for the same-sign dilepton channel are negligible after imposing a set of selection cuts, which however keep the signal rates almost unchanged. There are also other channels available for the $Z'\to NN$ process (e.g. $N\to \nu h/\nu Z$), which can further increase the signal significance of the model. Based on this consideration, we make the zero background assumption in this subsection to derive the projected sensitivities for the thermal leptogenesis parameter space of the $B-L$ model. The $CP$ asymmetry $\epsilon$ defined in \Eq{epsilon} is related to the asymmetry between the positive and negative same-sign dileptons as~\cite{Blanchet:2009bu}
\begin{align}
\epsilon=\frac{1}{2}\abs{\frac{N_+-N_-}{N_++N_-}},
\end{align}
where $N_\pm$ denotes the event number of the $e^\pm e^\pm W^\mp W^\mp$ final state. If we observe no asymmetry of the same-sign dileptons at colliders, upper limits can be put on the $CP$ asymmetry. This is done by assuming the number of the signal events to follow a Poisson distribution, and hence at 1$\sigma$ confidence level~\cite{ParticleDataGroup:2020ssz}
\begin{align}
N_+ =  \ave{N_+} \pm \Delta N_+=\ave{N_+} \pm \sqrt{ \ave{N_+}},
\end{align}
yielding a sensitivities of
\begin{equation}\label{limit}
\epsilon_{\rm min}=\frac{1}{2 \sqrt{\ave{N_+}}}.
\end{equation}
The $CP$ asymmetry is detectable if $\epsilon>\epsilon_{\rm min}$.

With \Eq{limit} in hand, one is able to test resonant leptogenesis at a specific collider environment. Given a set of $(M_{Z'},M_N,g_{B-L})$, the $\epsilon$ required by leptogenesis is derived by \Eq{eq:kappa}. On the other hand, the corresponding $N_+$ and hence $\epsilon_{\rm min}$ is also available by collider simulation of the $Z'\to NN$ process. If $\epsilon>\epsilon_{\rm min}$, then this parameter setup is detectable at the collider. The projected detectable parameter space with a fixed $g_{B-L}=0.8$ is plotted in Fig.~\ref{fig:asy_zp} for the HL-LHC (13 TeV, 3 ab$^{-1}$) and two setups of muon colliders (10 TeV, 10 ab$^{-1}$ and 30 TeV, 90 ab$^{-1}$) with different colors. The region shaded in dark blue is for $\epsilon>1$ or $M_N<T_{\rm sph}$ that the leptogenesis is not feasible. As we can see, the HL-LHC is able to probe leptogenesis for $M_{Z'}\lesssim12$ TeV with a RHN mass $M_N\lesssim600$ GeV, corresponding to the $CP$ asymmetry $\epsilon\sim0.2$. While for the muon colliders, due to the higher production rates and integrated luminosities, smaller $\epsilon\sim10^{-3}$ can be probed and a much higher mass region can be reached. At the 10 TeV muon collider, if the RHN is as light as a few hundreds of GeV, leptogenesis with $M_{Z'}$ up to 28 TeV can be probed via the off-shell process $\mu^+\mu^-\to Z'^*\to NN$. The 30 TeV muon collider can cover a wider parameter space via the on-shell $Z'Z$ or $Z'\gamma$ production for $M_{Z'}$ up to 30 TeV, and via the off-shell production for $M_{Z'}$ up to 100 TeV. We also notice that the parameter space for $M_{Z'}\sim\mathcal{O}(10~{\rm TeV})$ and $M_N\sim2$ TeV is not reachable for all collider setups, because the $\epsilon$ needed for BAU is too small (see the contours in Fig.~\ref{fig:blz}). However, a band in this region can be filled up by the 10 TeV muon collider for $M_{Z'}\sim10$ TeV, due to the resonant enhancement of the $Z'Z/Z'\gamma$ process. In summary, a considerable fraction of parameter space can be probed at the future muon colliders.

\section{Conclusion}\label{sec:con}

Both the BAU and tiny neutrino mass can be explained simultaneously within the framework of the thermal leptogenesis via the RHNs. If two of the RHNs are nearly degenerate, the scale of leptogenesis can be as low as $\sim$ 100 GeV and hence accessible at current and future colliders. In this article, we take the $B-L$ model as an example to perform a study on the collider probe of the leptogenesis mechanism. In this model, the new gauge boson $Z'$ mediates a scattering process that tends to washout the BAU. Therefore, if both $Z'$ and $N$ are at $\mathcal{O}({\rm TeV})$ scale, a sizable $CP$ asymmetry $\epsilon$ is needed to explain the observed BAU. In other words, successful (resonant) leptogenesis can be realized via TeV scale $Z'$, $N$ and a sizable $\epsilon$, which are testable at the colliders.

We choose the channel $Z'\to NN\to \ell^\pm\ell^\pm+{\rm jets}$ to probe leptogenesis, as this process involves the three key ingredients of the scenario: a new gauge boson $Z'$, the RHN $N$, and the $CP$ violation $\epsilon$ from the asymmetry between $N(\ell^+\ell^+)$ and $N(\ell^-\ell^-)$. Three collider setups are considered: 13 TeV LHC with 3 ab$^{-1}$ and 10 (30) TeV muon colliders with 10 (90) ab$^{-1}$. At the LHC, the $Z'$ is produced via the Drell-Yan process; while at the muon colliders, $Z'$ can be produced in association with a $Z'$ or $\gamma$, or via the off-shell $s$-channel fusion. We have shown that the backgrounds are safely negligible after a set of selection cuts, and hence the sensitivities for $\epsilon$ can be derived under the zero background assumption. Our quantitative study shows that leptogenesis in the $B-L$ model can be probed at the LHC for $6~{\rm TeV}\lesssim M_{Z'}\lesssim12$ TeV and $M_N\sim{\rm TeV}$. At the muon colliders, due to the higher signal rates, the projected probe limits can cover $M_{Z'}$ up to 30 TeV or even higher. For the 30 TeV muon collider, the leptogenesis can be probed up to $M_{Z^\prime} \sim 100$ TeV via the off-shell production of $Z'$. Our work demonstrates that the muon colliders can serve as a machine to efficiently probe the early Universe dynamics.

\section*{Acknowledgements}
We thank Frank Deppisch and Suchita Kulkarni for useful early discussions. We are grateful to Steve Blanchet, Zackaria Chacko and Rabindra Mohapatra for very helpful discussions. KPX is supported by the University of Nebraska-Lincoln.

\bibliographystyle{JHEP}


\section{Introduction}

\subsection{Knowledge Preservation}
A basic goal of most approaches to artificial general intelligence is the ability to learn multiple tasks without overwriting what has previously been learned (i.e. catastrophic forgetting, or CF). The potential utility of learning without forgetting, or 'knowledge preservation', is twofold: first, the system can return to previously learned tasks after learning something new. In addition, bootstrapping previous knowledge may allow for faster learning of a novel task through transfer learning. Both of these potential benefits of preserving learned information may serve to conserve resources. Allowing the same computational hardware to be ‘multiplexed’ for another task, rather than requiring allocation of distinct resources which reduces the required energy resources involved in the learning process. In the former, multiplexed information is stored in such a way that the system can switch between already learned tasks while retaining performance. In the latter case, information is accessed for learning new tasks in a contextually relevant manner. Multiplexed information is thus stored and retrieved in a way that facilitates the conservation and use of prior learning, broadening the ability of a learning system to multiple task domains.

\subsection{Previous Approaches}
A significant body of recent work has explored approaches to CF and bootstrapping in relation to multi-layer neural networks or deep learning (DL) \cite{xie2020artificial}. Learning in these systems consists of updating the weights between nodes according to their respective partial derivatives calculated from a loss function (i.e. backpropagation). With the recent implementation of auto-differentiation and parallel computing devices, progress in backpropagation-based DL has rapidly emerged as the dominant approach in AI \cite{schmidhuber2015deep}. Previous approaches to overcoming CF in DL models have involved strategies with regard to the learning process itself. Most commonly in the case of bootstrapping, ensemble networks are re-trained on resampled datasets with little improvement \cite{nixon2020bootstrapped}. Transfer learning is also used as a bootstrapping method in which previously-learned features are imported as the starting weights for later update \cite{weiss2016survey}. CF approaches are more diverse, mostly involving modification to the weight update methods. Notably, Elastic Weight Consolidation has shown promising results in overcoming CF in sequences of Atari tasks by slowing down updates on previously encoded weights \cite{kirkpatrick2017overcoming,aich2021elastic}. However, these approaches have not been widely adopted, likely because the techniques do not demonstrate knowledge preservation on a level necessary for general intelligence. We suspect this is in part due to the fundamental nature of backpropagation in DL models to rewrite learning; where weights are updated most heavily towards current learning at the expense of weights learned in more distant tasks.

\subsection{Network Topology}
Another potential approach that has thus far been largely neglected  is to consider the topology of the learning system — that is, its ‘architecture’. While there is a significant body of work on architecture optimization for other metrics of learning \cite{lukovsevivcius2012reservoir,elsken2019neural,goudarzi2012emergent,vanschoren2018meta,lee2014evolution,williams2017evolution,scholkopf2021toward}, to date there has been little work considering architectural choices in relation to CF and/or bootstrapping. Fortunately, hundreds of million years of evolution has generated highly conserved, vertebrate and mammalian brain architectures that have been optimized with regard to intelligent functions including, it is reasonable to assume,  the preservation of learning without CF.

One potentially promising property of biological  nervous system architectures, which may serve to overcome the pitfalls of CF, is the presence of modular architectures such as the cortex and subcortical structures. Modularity is a highly conserved property that is present across many organisms with general learning abilities. Complex neural architectures are thought to have emerged in the Cambrian explosion, 540 million years ago, with records of both camera eyes and the emergence of animals with complex behaviors \cite{trestman2013cambrian}. In contemporary species, the four basic brain structures (Telencephalon, Cerebellum, Diencephalon, Mesencephalon) have been preserved across nearly all mammalian species and most avian, reptile, shark brains, and some fish brains \cite{avianB,Yopak12946,Kotrschal2020,larsell1926cerebellum,b9}. Insects have their own highly preserved modular nervous-system  structures as well \cite{ITO2014755}.

In humans, modularity is highly pervasive. For example, the cortex processes incoming sensory information while the basal ganglia controls motor movements in behavioral learning \cite{ring2002neuropsychiatry,sutherland1988contributions}. These structures are specialized not only with regard to the types of information they  processed (e.g. photon patterns for V1 stream vs. layers 5/6 cortical output from cortex), but also with regard to their neuron types as well as their morphological structures. Communication between the two structures occurs in specialized. limited routes, cortico-basal ganglia networks, which serve to integrate information between the two regions \cite{lalo2008patterns}. The result of this particular organization is a cortex that is optimized for sensory information, while the basal ganglia is optimized to interpret cortical information for motor movement implementation \cite{coward2005system}.

Here, we propose that such modular architectures may also serve to preserve learning during the update process by separating different types of computation (specialization) and storing the learning across un-updatable boundaries (segregation) which communicate through specific upstream pathways (integration). In artificial systems, this modularity  may be modeled in the form of segregation of weight updates based on local update rules within a module. Previous work on artificial systems has looked at the effects of modularity for generalization, but using different learning metrics \cite{goyal2019recurrent,goyal2020inductive,werle2020investigating}. Promising results on the effects of modularity for generalizing to out-of-distribution data has been shown in the Badger architecture, which relies on nested self-sufficient computational units which communicate for a shared objective \cite{rosa2019badger}. Badger is an agent based model in the sense that each agent, or small self-sufficient network, updates in an internal loop, while the outer loop optimizes the collection of homogeneous agents. However, badger does not incorporate within-agent modularity which we believe is essential to solving the knowledge preservation problem.

While we have attempted to make the specific case for modularity as an example of a potentially important design principle in natural systems, we propose that knowledge preservation should serve as a general guiding framework in architecture design that may give rise to dynamical learning systems capable of some level of generality. The key idea is modeled from the human brain where it is likely that network topology is perhaps at least, if not more important than later learning algorithms in achieving this goal. This allows us to avoid ending up with a system that is doing the bulk of the learning process after it has been configured. In a system that prioritizes initial connectivity, this manifests as selecting appropriate modularity over a variety of different domain tasks. Rewriting existing learning prevents sharing of computational resources for novel tasks. This is often the problem in designing AGI from theoretical systems that don’t take resource constraints into account, such as universal Turing machines and AIXI \cite{turing2009computing,hutter2003gentle}.

\section{Possible Implementation Approaches}
Many methods could be conceived for discovering the inductive biases that modularity has to offer learning system architectures. Crucially, there is a tradeoff between initial connectivity and later weight modification. A strong dynamical learning system would utilize initial topology search as a prior learning experience to later be exploited during the weight modification learning phase.

Consider regions separated by un-updatable boundaries that delineate the extent of backpropagation. These modules could optimize without overwriting each other, while still being able to share resources when relevant (i.e. bootstrapping). This approach also limits rewriting prior learning for the sake of the current learning, where modules are activated at different times depending on their relevance to the current tasks. Furthermore, modules dedicated to computational tasks, instead of behavioral tasks (e.g. memory, sensory, self-activation, etc.) may serve as a more biologically inspired approach, since this is largely what we observe between cortical and subcortical structures.

Evolutionary searches are particularly useful in discovering unique solutions to optimization problems, which can be extended here. Although computationally extravagant, one promising candidate would be an evolutionary approach that couples network topology and embodiment search to different learning environments that also evolve over time. Derivatives of this approach would be useful, especially those focusing on modularity and learning. The advantage to this method is that appropriate submodules will likely be paired with the particular sensory receptors given to the agent and the types of learning that needs to be delegated throughout the system in ways that may not be currently known.

Other approaches include intelligently designed architectures where architecture is random, learned, selected by an operator, and most commonly in combination. Random networks, such as \cite{gershenson2004introduction}, have seen interesting results and may offer some potential insights into the core properties of learning, however they will likely  have to be harnessed in a constructive way to get desired behavior, as in \cite{lukovsevivcius2012reservoir}. Similar to random learning, methods exist for crystallizing order from disorder in unorganized networks as in the case of \cite{gutowitz1991cellular,hebb1949organisation}, leading to module-like structures. Still, there are other conceivable approaches involving learning functions that select sequences of computational architecture pathways and more. More extensive studies could iteratively search through all possible configurations of N-node networks to discover what modularity benefits particular learning environments. Even hand-designed simple models that start to test the questions of module design and learning without forgetting would be useful in answering the questions presented here.

\section{Benchmarking Generality}
In order to assess a model's knowledge preservation, it must be able to demonstrate that it can learn a sequence of tasks in a variety of domains and return to those tasks with improvements in computational costs and/or learning time. Below we describe initial tests to measure generality in terms of learning without forgetting and bootstrapping. However, as generality scales, more advanced testing is likely necessary.

\subsubsection{Knowledge Preservation Metrics}
Learning without catastrophic forgetting is typically  measured by performance on a sequence of tasks. Here, we propose it is imperative to test on sequences in different domains (e.g. vision classification and time-series forecasting) to demonstrate generality over narrow-learning.

\subsubsection{Boostrapping A-Priori Knowledge Metrics}
Measuring bootstrapping a-priori knowledge is slightly more challenging due to the effects seen in transfer learning. The goal here is to see an increased speed in learning task B after learning task A. Again, in trying to capture generality, it is useful to compare tasks from different domains so the distribution of features in the data is different enough that it's possible to measure using contextual knowledge instead of the re-use of features. This later point is what occurs in transfer learning. As in the case of the popular technique of using pre-trained ResNet on Imagenet and vision based tasks, the transferred weights capture many features which would have to have been learned in the secondary network either way \cite{weiss2016survey}. Here, we also want to be able to use learned features to aid in learning a new task, but we want to ensure that these are general features that can be maximally shared across all tasks the system may be expected to learn. Teasing out this difference is aided in investigating how data flows through different parts of the system’s topology. For example, the Fusiform Face Area is shown to activate when distinguishing between faces and also in bird experts for distinguishing between birds \cite{tong2008fusiform}. This multi-use activation of the same subnetwork is one example of bootstrap learning. We might expect similar submodules to be redeployed across different kinds of learning domains.

\subsubsection{Challenges}
We recognize that adding modularity alone will probably not multiplex information in a way that provides human-level general intelligence. The works of Coward, Buszaki, and Turing suggest algorithms employed by the human brain may give rise to this generality based on eliciting sequences of computation across different modules  \cite{coward2005system,gyorgy2019brain,turing1969intelligent}. Furthermore, special message passing algorithms between the regions would benefit from that of publish and subscribe models; that regions in the system can communicate with each other when necessary.

\section{Conclusion}
We have highlighted topology as an important topic for future research, as well as the need for algorithms that allow for learning without rewriting previous learning. We suspect modularity in biological intelligent systems plays a crucial role in achieving this need while conserving resources. This paper also contributes to discussions on AGI agent architectures by suggesting some useful considerations on how to test for generality in regards to CF and bootstrapping. Works from \cite{wang2019paired,rosa2019badger,koutnik2013evolving,stanley2019designing,aich2021elastic} have attempted to tackle these issues separately, while others \cite{silver2021reward,wang2021modular} have come from the angle of single-use systems, all of which overwrite prior learning to a unsatisfactory degree. Still, some works have come from the position of down-scaling theoretical models without resource constraints \cite{goertzel2007artificial,hutter2003gentle,turing2009computing}. While these approaches are useful, we find the proposed approach is more appropriate for the pursuit of general intelligence. We believe this paper provides a framework for investigating computationally tractable general learning systems. Promising directions involve evolutionary and intelligent topology design that utilizes modularity to minimize rewriting prior learning. We hope this paper serves as a tool for those interested in benchmarking early progress in generality. There is room for many future works to provide insight into the approach presented here as well as how to measure generality in terms of bootstrapping and learning new without forgetting old. We have tried to clarify how modularity may offer a method of learning without rewriting and the role it plays in topology design. Much work needs to be done to improve that understanding and find ways to appropriately incorporate modularity in AGI frameworks.

\bibliographystyle{splncs04}


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

The Higgs boson interactions with top quarks and electroweak gauge bosons are one of the most important subjects of studies at the LHC. The strength of the top-Higgs Yukawa coupling can be probed directly without making any assumptions on its nature in the associated production processes involving both a top quark and a Higgs boson in the final state. Alongside the intensively studied associated production with a top-antitop quark pair, $t \bar t H$, the associated process of a single top quark and a Higgs boson ($t H$) production has been attracting a lot of attention lately. Although in the Standard Model (SM) the $tH$ production rate at the LHC is only one tenth of the rate for the $t \bar t H$ production ~\cite{Maltoni:2001hu},  the process is of great interest.  It offers a unique possibility to directly measure not only the value of the Yukawa coupling but also a relative sign between it and the coupling of a Higgs boson to a $W$ boson. The $tH$ cross section receives contributions from diagrams with $H$ radiated by the $t$ quark as well as from diagrams with $H$ radiated by the $W$ boson. The value of the cross section depends strongly on the interference between these diagrams. In fact, it is the interference between diagrams involving the two couplings that is responsible for the small value of the cross section in the SM.  For the same reason the $tH$ process offers a very sensitive probe of new physics interactions. The small cross section notwithstanding, first experimental searches for the process have been already conducted at the LHC~\cite{Khachatryan:2015ota,Sirunyan:2018lzm} and the rate for $H$ boson production in the $tH$ process is included in the measurements together with the $t \bar t H$ process~\cite{Aad:2020mkp, Aad:2020ivc, Sirunyan:2020icl}.

Given on one side the importance of the $tH$ production for new physics searches and, on the other side, much larger statistical samples and expected better understanding of the systematic uncertainties of the future LHC measurements, the precision with which the theoretical predictions are known becomes a central issue. At lowest order in perturbation theory, the production mechanisms for the associated production of the Higgs boson with a single top quark can be classified according to the virtuality of a $W$, which is either exchanged or emitted in the process. Correspondingly, the three mechanisms are the $t$-channel ($ qb \to t H q'$), $s$-channel ($ q \bar q' \to t H \bar b$) and the $tW$ associated channel ($g b \to t H W$). Strictly speaking, this classification is only valid assuming five active quark flavours, i.e. in the so-called five-flavour scheme (5FS) which involves $b$ quarks in the initial state. In the four-flavour scheme (4FS), initial state $b$-quarks can only appear as a result of a $g \to b \bar b$ splitting. First calculations of the next-to-leading order (NLO) QCD correction for the $t$-channel (in 4FS and 5FS) and $s$-channel production in the SM were reported in~\cite{Demartin:2015uha}, whereas an analysis of the NLO QCD corrections to the $t$-channel process in the framework of SM Effective Field Theory (SMEFT) was undertaken in~\cite{Degrande:2018fog}. Recently, NLO QCD+EW precision for the $tHj$ production was reached in the combined analysis of all channels in 5FS ~\cite{Pagani:2020mov}.

With the full NNLO calculations being technically out of reach for the $ 2 \to 3$ processes considered here, a question remains if at least some part of the higher order corrections can be investigated already now. Higher-order corrections due to soft gluon emissions provide a class of corrections which can be identified and accounted for in a systematic manner, as shown in the case of the  $ t \tb H$ or $t \bar t Z/W$ production~\cite{Kulesza:2015vda, Kulesza:2016vnq, Kulesza:2017ukk, Kulesza:2018tqz, Kulesza:2020nfh, Li:2014ula, Broggio:2015lya, Broggio:2016zgg, Broggio:2016lfj, Broggio:2017kzi, Broggio:2019ewu}. In particular, it has been reported in~\cite{Kulesza:2017ukk, Kulesza:2018tqz, Kulesza:2020nfh}, that adding soft gluon corrections calculated at the next-to-next-to-leading-logarithmic (NNLL) accuracy leads to a significant reduction of the scale uncertainty and a correction to the $ t \tb H$ or $t \bar t Z$ total cross section of up to 20\%, depending on the choice of the central scale. Of course, one cannot expect a similar conclusion to automatically apply to the single associated top production, as for example the initial state and final state partons are different in both cases. In particular, in contrast to  the $tH$ production, the $ t \tb H$ process involves the gluon-gluon channel, so it is reasonable to expect smaller impact of soft gluon corrections in the former case.   Nevertheless, since the processes are often analysed together by the experiments, theoretical predictions for the two processes should be known with the same accuracy. In a similar, in this respect, case of the associated heavy gauge boson production, the $t \bar t W$ process, taking place at LO only in the quark channel, is known at the NLO+NNLL accuracy, i.e. the same level of accuracy as reached for the $t \bar t Z$ process.

So far, only the predictions for the $t$-channel $tH$ production have been considered in the context of soft gluon emission.  More specifically,  in~\cite{Forslund:2020lnu,Forslund:2021evo} contributions to the NNLO corrections due to such emission were investigated. That study used the so-called one particle inclusive (1PI) kinematics. Performing full resummed calculations, i.e. obtaining numerical predictions which include soft gluon corrections to all order in perturbation theory in the direct resummation framework (i.e. in Mellin space) is known to be notoriously difficult in 1PI, see e.g~\cite{Hinderer:2018nkb}. For processes involving  jet production it has been achieved only for inclusive enough quantities~\cite{deFlorian:2007fv}. Notably, 1PI threshold resummation for single-inclusive jet production has been carried out in the SCET framework~\cite{Liu:2017pbb}.

In this work, we investigate the application of the direct Mellin resummation techniques to the $s$-channel $tH$ production. While the contribution from the $s$-channel to the $tH$ production is much smaller than the one from the $t$-channel, it proves a useful setting for an exploratory study. The $s$-channel does not involve initial state $b$-quarks and therefore it remains largely insensitive to the FS choice, which is not the case for the $t$-channel predictions. Given the known problems with performing fully resummed calculations in 1PI kinematics, we focus instead on the three-particle invariant mass kinematics, which has successfully been applied to the production of three massive particles, some carrying colour,  in the finals state~\cite{Kulesza:2015vda, Kulesza:2016vnq, Kulesza:2017ukk, Kulesza:2018tqz, Kulesza:2020nfh}. Application of these kinematics to the process in question requires considering kinematical information on a jet $j$ produced  together with a top quark and a Higgs boson.
This, in turn, brings up the question of an appropriate treatment of the jet. In the presence of particles other than massless partons in the final state at LO, 1PI kinematics allows to perfom calculations of observables which do not directly involve the jet by treating it inclusively, i.e. as unobserved. In other words, the threshold variable in 1PI kinematics does not directly involve the momenta of the massless parton. This is not the case in the invariant mass kinematics. Following~\cite{Kidonakis:1998bk}, in this kinematics the jet can be treated as massless or massive at the partonic threshold, leading to different structures of the logarithmic corrections. Massless jet approach results in double logarithmic contributions in the resummed factor describing the jet. Giving up the assumption of a vanishing jet mass at threshold leads to only single logarthmic contributions, albeit with coefficients depending on the jet size parameter. In the following, we consider the impact of threshold corrections on the prediction in both approaches. Our resummed results, take into account leading logarithmic (LL) and  next-to-leading logarithmic (NLL) terms together with ${\cal O} (\alpha_S)$ non-logarithmic contributuions which do not vanish at threshold.

The rest of the paper is structured as follows: in Section 2 we present the theoretical framework and the corresponding analytical expressions. Section 3 focuses on the numerical predictions. First we discuss the relation between the full NLO cross sections and the expansion of the resummed cross section truncated at NLO. Then we present the resummed NLL predictions matched to the NLO results. Section 4 contains a summary  of our results.

\section{NLL resummation in the triple invariant mass kinematics: theoretical framework}
\label{s:theory}

In the following we consider the corrections due to gluon emissons in the limit of the invariant mass $Q$ of the top quark ($t$), the Higgs boson ($H$) and the jet ($j$) becoming close the partonic center-of-mass energy $\hat s$,  {$Q^2 \equiv  (p_t +p_H+ p_j)^2 \to \hat s$. This corresponds to the limit $\hat \rho \to 1$ of the threshold variable $\hat \rho \equiv Q^2/\hat s$. In the $s$-channel process, the final state jet is initated by the $b$-quark. In this work we treat the $b$ quark as massless, i.e. we use a jet function to describe the (quasi-)collinear emission from the final state quark. Further, we follow~\cite{Kidonakis:1998bk} regarding the treatment of the final state jet, and apart from an assumption of a massless jet with $p_j^2=0$, we also consider the case of a non-vanishing jet mass. For massless jets  $Q^2 =(p_t+p_H)^2 + 2 (p_t+p_H) \cdot p_j$.

In general, the construction of the resummation formalism for a $2 \to 3$ process builds up on the earlier work for the $2 \to 2$ processes~\cite{Kidonakis:1998bk, Contopanagos:1996nh, Kidonakis:1998nf, Bonciani:2003nt}.  Here we adapt the expressions for the $t \bar t H (W,Z)$  from~\cite{Kulesza:2017ukk, Kulesza:2018tqz} to account for the presence of a final state jet instead of a top quark. The factorization of the soft emission, necessary to perform resummation, takes place in Mellin space, where the Mellin moments $N$ of the partonic cross sections are taken w.r.t. the variable $\hat \rho$.   At the next-to-leading logarithmic (NLL) accuracy, the resummed expression for the partonic cross section has a schematical form
\begin{eqnarray}
\tilde\hat \sigma^{{\rm (NLL)}}_{i j}=  \mathrm{Tr} \left[\mathbf{H}_{ij \to tHk} \,\, \mathbf{S}_{ij \to tHk} \right] \Delta_i \Delta_{j} J_k \,
\label{eq:res:fact_diag_NLL}
\end{eqnarray}

Since the soft radiation is coherently sensitive to the colour structure of the underlying hard process $i j \to t H k$,  the hard function $\mathbf{H}$, representing hard off-shell dynamics of the process, and the soft function $\mathbf{S}$, describing pure soft radiation, are matrices in color space. The $\Delta_i$ jet functions contain logarithmic contributions due to (soft-)colinear radiation from the incoming quarks or antiquarks, while $J_k$ corresponds to the final state jet function describing emissions of a quark $k$.

In our calculations, following~\cite{Kidonakis:1998bk}, we define the final state jet four-momentum as the sum over the four-momenta of particles flowing into a cone of half-aperture (angular radius) $\delta$ around the jet axis. This approach corresponds to the small cone approximation~\cite{Furman:1981kf, Aversa:1989xw} used in jet calculations.  In the small cone approximation the $\delta$ parameter is related to the jet radius $R$ through $\delta=R/\cosh(\eta)$, where $\eta$ is the pseudorapidity of the jet. The approximation is known to work very well even up to relatively large values of $R \sim 0.7$~\cite{Jager:2004jh,Mukherjee:2012uz}.

The incoming jet functions are independent of the process and read
\begin{equation}
\Delta_i   =  \exp\left[\int_0^1 \text{d}z \, \frac{z^{N-1}-1}{1-z}  \left\{ \int^{(1-z)^2 Q^2}_{\mu_F^2}\frac{dq^2}{q^2} A_i (\alpha_{\rm s}(q^2))\right\} \right]\,,
\end{equation}
where $\mu_F$ is the factorization scale. The outgoing jet function has a different form depending on how the jet is treated~\cite{Kidonakis:1998bk}. Specifically, if the jet is required to be massless it is given by
\begin{equation}
J_k  =  \exp\left[\int_0^1 \text{d}z \, \frac{z^{N-1}-1}{1-z}  \left\{\int_{(1-z)^2 Q^2}^{(1-z) Q^2}\frac{dq^2}{q^2} A_k(\alpha_{\rm s}(q^2))+\frac{1}{2} B_k(\alpha_{\rm s}((1-z) Q^2))\right\}\right]\,.
\label{eq:jet:massless}
\end{equation}
Otherwise, for massive jets we have
\begin{equation}
J_k  = \exp\left[\int_0^1 \text{d}z \, \frac{z^{N-1}-1}{1-z}  C^{\rm M}_k \left(\alpha_{\rm s}((1-z)^2 Q^2)\right)\right]\,.
\end{equation}
As discussed in~\cite{Kidonakis:1998bk}, the outgoing jet function contains double threshold logarithms if the jets are treated as massless, but only single logarithms in the case of a non-vanishing jet mass. Another important difference is that in the former case the jet function $J_k$ does not depend on the $\delta$ parameter, in contrast to the logarithmic dependence on $\delta$ which enters through the coefficient $C^{\rm M}_k$ in the latter case. The coefficients $A_i$, $B_i$, $C^{\rm M}_i$ are all perturbative series in $\alpha_{\rm s}$,
$$
F_i = \frac{\alpha_{\rm s}}{\pi} F_i^{(1)} + \frac{\alpha_{\rm s}^2}{\pi^2} F_i^{(2)} + \dots\,,
$$
with $F_i$ standing in for any of the coefficients. In order to achieve the NLL accuracy, $A_i^{(1)}, A_i^{(2)}, B_i^{(1)}$ and $C_i^{\rm M, (1)}$ are needed. At LO, the $tH$ production involves only quarks in the initial and final state, i.e. $i=q$. The LL $A_q^{(1)}$ coefficient as well as the NLL $A_q^{(2)},\  B_q^{(1)}$ coefficients of the quark jet functions are well known~\cite{KT, CET}
\begin{eqnarray}
A_q^{(1)} &=& C_F, \qquad A_q^{(2)} = {1 \over 2} C_F K\quad   {\rm with} \quad K=C_A\left(  {67 \over 18} -{\pi^2 \over 6}\right) -{5 \over 9} N_f \nonumber \\
B_q^{(1)} &=& \gamma_q = -\frac{3}{2} C_F.
\end{eqnarray}
For jets massive at threshold, our calculations return
\begin{equation}
C^{\rm M, (1)}_q = C_F \log \left( {Q^2 \over \delta^2 E_q^2} \right) \,,
\end{equation}
with the $\delta$ dependence entering logarithmically, as expected.

The function $\mathbf{S}$ sums logarithmic contributions from soft wide-angle emission. In Eq.(\ref{eq:res:fact_diag_NLL}) it is already rescaled by the inverse of the soft function contributions for the incoming and outgoing lines, to avoid double counting of soft contributions included in the $\Delta_i$ and $J_k$ functions used in this work. In the most general case, the solution of the evolution equation  for $\mathbf{S}$~\cite{Kidonakis:1998bk} involves path-order exponentials of the integrals over the soft anomalous dimension $\mathbf{ \Gamma} ={ \alpha_{\rm s} \over \pi} \mathbf{\Gamma}^{(1)} + \dots$. To reach the NLL accuracy, we need to know the one-loop soft anomalous dimension $\mathbf{ \Gamma}^{(1)}$.  For the process $q \bar q^{\prime} \to t H k$ with $q, \bar q^{\prime},k$ corresponding to massless quark lines we find
\begin{eqnarray}
\Gamma^{(1)}_{q \bar q^{\prime} \rightarrow t H k}&=&C_F\mathbbm{1} \left[\log\left(\frac{s_{t k}-m_t^2}{m_t Q}\right)-\frac12\left(1-i\pi\right)\right]\nonumber\\
&&-\frac12 \boldsymbol{C_I}\left[\log\left(\frac{s_{tk}-m_t^2}{m_t Q}\right)-\log\left(\frac{t_{qk}\left(t_{\bar q^{\prime} t}-m_t^2\right)}{m_t Q^3}\right)\right]\nonumber\\
&&-\mathbf{T}_{qt}\log\left(\frac{t_{\bar q^{\prime} k}\left(t_{q t}-m_t^2\right)}{t_{q k}\left(t_{\bar q^{\prime} t}-m_t^2\right)}\right) \,,
\label{eq:SAD}
\end{eqnarray}
where we have $s_{ij}=(p_i+p_j)^2$,  $t_{ij}=(p_i-p_j)^2$, $\boldsymbol{C_I}=\text{Diag}\left(0, C_A\right)$ and
\begin{eqnarray}
\mathbf{T}_{qt}&=&
\begin{pmatrix}
0 & -\frac{C_F}{2C_A} \\
-1 & -2C_F+\frac{C_A}{2}
\end{pmatrix}
\end{eqnarray}
or
\begin{eqnarray}
\mathbf{T}_{q\bar{t}}&=&
\begin{pmatrix}
0 & \frac{C_F}{2C_A} \\
1 & 2C_F-C_A
\end{pmatrix}\,
\end{eqnarray}
for the $q \bar q^{\prime} \to \bar t H k$ process. Our results  for  $\mathbf{ \Gamma}^{(1)}$  are obtained in the singlet-octet $s$-channel basis and agree with the result recently published in the literature~\cite{Forslund:2020lnu}.

Apart from knowing the soft anomalous dimension, evaluation of the soft function at any scale also requires  knowing the boundary condition $\mathbf{\tilde S}$   in the solution of the renormalization group equation at $\mu_R = Q/\bar N$, with $\bar N = N e^{\gamma_E}$ and $\gamma_E$ the Euler constant, cf. Eq.~(\ref{eq:soft:R:evol}) below.  $\mathbf{\tilde S}$ is a purely eikonal function which can be calculated perturbatively
\begin{equation}
\mathbf{\tilde S}= \mathbf{\tilde S}^{(0)} + \frac{\alpha_{\rm s}}{\pi} \mathbf{\tilde S}^{(1)} + \dots
\label{eq:softexp}
\end{equation}
The leading-order soft function is given by
\begin{eqnarray}
\mathbf{\tilde S}^{(0)}=
\begin{pmatrix}
C_A^2 & 0 \\
0 & \frac{C_A C_F}{2} \,
\end{pmatrix}.
\end{eqnarray}
With $\mathbf{\tilde S}^{(0)}$ and $\mathbf{\Gamma}^{(1)}$, NLL precision can be reached. Any higher precision  requires also the expression for  $\mathbf{\tilde S}^{(1)}$. For massless initial state \mbox{$i=q,\bar q^{\prime}$}  and final state $k$ quarks we obtain finite contributions in the form
\begin{eqnarray}
\mathbf{ \tilde S}^{(1)}&= &\mathbf{\tilde S}^{(0)} \left\{ C_F\mathbbm{1}\left[S_{tt}-2S_{tk}\right] +\boldsymbol{C_I}\left[S_{q\bar q^{\prime}}+S_{tk}-S_{qk}-S_{\bar q^{\prime} t}\right] \right. \nonumber\\
&&+\left. 2\mathbf{T}_{qt}\left[S_{qt}+S_{\bar q^{\prime} k}-S_{qk}-S_{\bar q^{\prime} t}\right] \right\}
\end{eqnarray}
with
\begin{eqnarray}
S_{q \bar q^{\prime}}&=&-\frac{\pi^2}{12}\nonumber\\
S_{tt}&=&\frac{1}{2\beta_t}\log\left(\frac{1+\beta_t}{1-\beta_t}\right)\nonumber\\
S_{tk}&=&- \frac{\pi^2}{24}+ \frac{1}{8} \log^2 \frac{1+\beta_t}{1-\beta_t} + \frac{1}{2} \left\{ \text{Li}_2 \left(1- \frac{2E_t E_k (1 -\beta_t)}{s_{tk}-m_t^2}\right) \right. \nonumber \\
&&\left. +\text{Li}_2 \left(1- \frac{2E_t E_k (1 +\beta_t)}{s_{tk}-m_t^2}\right) \right\}\nonumber\\
S_{it}&=&- \frac{\pi^2}{24}+ \frac{1}{8} \log^2 \frac{1+\beta_t}{1-\beta_t}  + \frac{1}{2} \left\{ \text{Li}_2 \left( 1 - \frac{E_t Q(1 +\beta_t)}{ m_t^2 - t_{it}}\right)\ \right. \nonumber \\
&& \left. +\text{Li}_2 \left( 1 - \frac{E_t Q(1 -\beta_t)}{ m_t^2 - t_{it}}\right) \right\}\nonumber\\
S_{ik}&=&- \frac{\pi^2}{12}  + \frac{1}{2}\text{Li}_2 \left( 1 + \frac{2E_k Q}{t_{ik}}\right)
\end{eqnarray}
and  $\beta_t=\sqrt{1-m_t^2/E_t^2}$.

Calculations of the soft function in the dimensional regularization (DR) scheme leading to the above results deliver also information on double and single poles in $\epsilon$. The one-loop soft function corresponding to $\mathbf{\tilde S}^{(1)}$ does not contain $1/\epsilon^2$ poles, while the coefficients of the $1/\epsilon$ terms return the soft anomalous dimension given in Eq.~(\ref{eq:SAD}), thus providing an additional check of our calculations.

In order to simplify the solution of the renormalization group equation for the soft function, and, correspondingly, the expression in eq.~(\ref{eq:res:fact_diag_NLL}), we write it in the color basis in which the one-loop soft anomalous dimension $\mathbf{ \Gamma}^{(1)}$, driving the evolution of the soft-function, is  diagonal. We denote all matrices in this basis with the subscript $R$.  The $IJ$ element of the diagonalised one-loop soft anomalous dimension is thus
$$
\mathbf{\Gamma}^{(1)}_{R, IJ} = \lambda_I^{(1)} \delta_{IJ}
$$
where $\lambda_I^{(1)}$ are the eigenvalues of $\mathbf{\Gamma}^{(1)}$. The soft function has then the form
\begin{equation}
\mathbf{S}_{R, IJ} = \mathbf{\tilde S}_{R, IJ} \exp\left[\frac{\log(1-2\lambda)}{2 \pi b_0}
\left(\left( \lambda^{(1)} _{I}\right)^{*}+\lambda^{(1)} _{J}\right)\right]\, ,
\label{eq:soft:R:evol}
\end{equation}
with $\lambda= \alpha_{\rm s}(\mu_R^2) b_0 \log N$, $b_0= (11C_A-2 nf)/(12\pi)$ and $\mu_R$ the renormalization scale.

The NLL precision of the resummed expression can be increased by including the ${\cal O}(\alpha_{\rm s})$ terms in the expansions of the hard function
$\mathbf{H} = \mathbf{H}^{\mathrm{(0)}} + \frac{\alpha_{\rm s}}{\pi}  \mathbf{H}^{\mathrm{(1)}} +\dots$
and the soft function $\mathbf{\tilde S}$, cf. Eq.~(\ref{eq:softexp}), as well as collinear non-logarithmic contributions. In the  literature on direct QCD resummation such terms are often collectively referred to as the "${\mathbf C^{(1)}}$" coefficient.
We will therefore refer to the precision obtained by adding the information on ${\mathbf C^{(1)}}$ to the NLL result as "NLLwC"~\footnote{This accuracy is formally equivalent to the "NLL'" notation used in the SCET literature on resummation.} and give the corresponding expressions below. In the $R-$basis we have
\begin{eqnarray}
 \frac{d \tilde\hat \sigma^{{\rm (NLLwC)}}_{ij }}{dQ^2}&&\hspace{-0.9cm}(N,Q^2,\{m^2\},\mu_{F}^2, \mu_{R}^2) =  \label{eq:res:fact_diag_NLLwC_R} \\
 &&\hspace{-0.9cm} \ (\mathbf{H}_{R,IJ} \mathbf{\tilde S}_{R,JI})(Q^2, \{m^2\},\mu_{F}^2, \mu_{R}^2, \delta) \nonumber \\
&\times&
\Delta_i(N+1, Q^2,\mu_{F}^2, \mu_{R}^2) \Delta_j(N+1, Q^2,\mu_{F}^2, \mu_{R}^2) \nonumber \\
&\times& J_k(N+1, \delta, Q^2, \mu_{R}^2 )\, \exp\left[\frac{\log(1-2\lambda)}{2 \pi b_0}
\left(\left( \lambda^{(1)} _{J}\right)^{*}+\lambda^{(1)} _{I}\right)\right], \nonumber
\end{eqnarray}
where $\{m^2\}$ stands for a set of mass parameters characteristic for the process.
The NLLwC precision is achieved with
\begin{eqnarray}
\mathbf{H}_{R}\, \mathbf{\tilde S}_{R}& = &\mathbf{H}^{\mathrm{(0)}}_{R} \mathbf{\tilde S}^{\mathrm{(0)}}_{R} \left( \mathbbm{1} + \frac{\alpha_{\rm s}}{\pi} \ \mathbf{C}^{(1)}_{R} \right) \label{eq:hardXsoft} \\
&=& \mathbf{H}^{\mathrm{(0)}}_{R} \mathbf{\tilde S}^{\mathrm{(0)}}_{R} + \frac{\alpha_{\rm s}}{\pi}\left[ \mathbf{H}^{\mathrm{(0)}}_{R} \mathbf{\tilde S}^{\mathrm{(1)}}_{R} + \mathbf{V}^{\mathrm{(1)}}_{R} \mathbf{\tilde S}^{\mathrm{(0)}}_{R}+ \mathbf{H}^{\mathrm{(0)}}_{R} \mathbf{\tilde S}^{\mathrm{(0)}}_{R} \left( 2 \tilde J_{in}^{(1)} + \tilde J_{out}^{(1)}  \right) \right]  \nonumber
\end{eqnarray}
If only the  $\mathbf{H}^{\mathrm{(0)}}_{R} \mathbf{\tilde S}^{\mathrm{(0)}}_{R}$ term, with its trace corresponding to the  leading order partonic cross section, is kept in the expression above, the precision of the resummed cross section in Eq.~(\ref{eq:res:fact_diag_NLLwC_R}) reduces to NLL.
Eq.~(\ref{eq:hardXsoft}) defines the NLLwC precision as the NLL precision supplemented with additional information included in the coefficient $ \mathbf{C}^{(1)}$. It consists of ${\cal O(\alpha_{\rm s})}$ contributions from the virtual corrections  $\mathbf{V}^{\mathrm{(1)}}$, collinear non-logarithmic contributions ${\tilde J}_{in}^{(1)}$ from incoming quark jets, as well as ${\cal O(\alpha_{\rm s})}$ non-logarithmic contributions ${\tilde J}_{out}^{(1)}$ from the outgoing jets. The virtual corrections $\mathbf{V}^{\mathrm{(1)}}$ are calculated numerically using the  aMC@NLO code~\cite{Alwall:2014hca}, whereas the collinear contributions are calculated analytically. In particular, for the outgoing quark jet with zero invariant mass we find the finite non-logarithmic contributions of the form
\begin{eqnarray}
{\tilde J}_{out}^{(1)} &=&C_F \left\{\frac{1}{4}\log^2 \frac{\mu_R^2}{Q^2}+\frac{3}{4}\log \frac{\mu_R^2}{Q^2}+ \frac{7}{4} + \frac{3}{4}\gamma_E - \frac{1}{2} \gamma_E^2 - \frac{ 5\pi^2}{12}
  - \log^2 \frac{2}{\delta} \right. \nonumber \\
& -& \log^2 \frac{E_k \, \delta}{Q}
+    \frac{1}{4} \Bigg[ \frac{Q^2}{Q^2 + E_k^2 \bar{N} \delta^2} + 3 \log \left(1+ \frac{Q^2}{E_k^2 \bar{N} \delta^2}\right) \nonumber \\
&+ &2 \log^2 \bigg(1+ \frac{Q^2}{E_k^2 \bar{N} \delta^2}\bigg) + 4 \text{Li}_2 \left( \frac{Q^2}{Q^2+E_k^2 \bar{N} \delta^2} \right) \nonumber \\
&+ &  \left. \left. \frac{\pi^2}{3}\frac{2\,Q^6 + 5 \,  E_k^2 \, \bar{N} \, Q^4 \, \delta^2 + 4 \, E_k^4 \, \bar{N}^2 \, Q^2\, \delta^4}{2\, (Q^2+E_k^2 \bar{N} \delta^2)^3} \right] \right\} \,,
\label{eq:jet1massless}
\end{eqnarray}
whereas in the massive case we have
\begin{eqnarray}
{\tilde J}_{out}^{{\rm M},(1)} =  C_F  \hspace{-1.5em} && \left(  2 -  \frac{\pi^2}{6} + \frac{1}{4}\log^2 \frac{\mu_R^2 }{Q^2} + \frac{3}{2}\log \frac{\mu_R}{\delta E_k} -\log^2 \frac{2}{\delta}  + \log^2 \frac{E_k \delta} {Q} \right. \nonumber \\
&&- \left.  2\gamma_E \log \frac{Q}{\delta E_k} \right)
\label{eq:jet1massive}
\end{eqnarray}
where $k$ corresponds to a massless final state quark.
In Eq.~(\ref{eq:jet1massless}), apart from terms of ${\cal O} (1)$ in $N$, we show additional contributions formally subleading in $N$ that, as we discuss in the next section, turn out to be relevant numerically. These contributions given in the square bracket arise for $\bar N \sim Q^2/(E_k^2 \delta^2)$ and  vanish in the limit  $\bar N \gg Q^2/ (E_k^2 \delta^2)$.

Our calculations of jet functions, leading to Eqs.~(\ref{eq:jet1massless}) and~(\ref{eq:jet1massive}), return coefficients of  the $1/\epsilon^2$ and $1/\epsilon$ poles corresponding to cusp and non-cusp contributions to the jet  anomalous dimension. They also reproduce the same coefficients of the double and single logarithms of $N$ as given by the first-order expansion of a massless jet function, cf. Eq.~(\ref{eq:jet:massless}). For a massive jet function they return only a single logarithm of $N$, accompanied by a coefficient logarithmic in $\delta$, in agreement with literature~\cite{Kidonakis:1998bk}.

The numerical predictions in the next section for the resummation-improved hadronic cross sections for the $s$-channel process $h_1 h_2 \to t H +X $  are obtained through matching the NLO cross section with the hadronic equivalent of Eq.~(\ref{eq:res:fact_diag_NLLwC_R}). In this way, the full information from the NLO calculations is preserved, while  double counting of the LO and NLO contributions included in the NLLwC expression avoided. More specifically,
\begin{eqnarray}
\label{hires}
\frac{d\sigma^{\rm (NLO+NLLwC)}_{h_1 h_2 }}{d Q^2}(Q^2,\{m^2\},\mu_{F}^2, \mu_{R}^2) &=&
\frac{d\sigma^{\rm (NLO)}_{h_1 h_2 }}{d Q^2}({Q^2},\{m^2\},\mu_{F}^2, \mu_{R}^2) \\ \nonumber
&+&   \frac{d \sigma^{\rm
		(res-exp)}_
	{h_1 h_2 }}{d{Q^2}}( Q^2,\{m^2\},\mu_{F}^2, \mu_{R}^2)
\end{eqnarray}
with
\begin{eqnarray}
\label{invmel}
&& \!\!\!\!\!\! \frac{d \sigma^{\rm
		(res-exp)}_{h_1 h_2 }}{d{Q^2}} ({Q^2},\{m^2\},\mu_{F}^2, \mu_{R}^2) \! = \!\!\! \sum_{i,j=\{q, \bar q^{\prime}\}}\,
\int_{\sf C}\,\frac{dN}{2\pi
	i} \; \rho^{-N} f^{(N+1)} _{i/h{_1}} (\mu_{F}^2) \, f^{(N+1)} _{j/h_{2}} (\mu_{F}^2) \nonumber \\
&& \!\!\!\!\!\! \! \times\! \left[
\frac{d \tilde\hat \sigma^{\rm (NLLwC)}_{q \bar q^{\prime}}}{d Q^2} (N,Q^2,\{m^2\},\mu_{F}^2, \mu_{R}^2)
-  \frac{d \tilde\hat \sigma^{\rm (NLLwC)}_{q \bar q^{\prime}}}{d Q^2} (N, Q^2,\{m^2\},\mu_{F}^2, \mu_{R}^2)
{ \left. \right|}_{\scriptscriptstyle({\rm NLO})}\, \! \right], \nonumber \\
\end{eqnarray}
where $ d \tilde\hat \sigma^{\rm (NLLwC)}_{q \bar q^{\prime}} / dQ^2 |_{\scriptscriptstyle({\rm NLO})}$ represents the perturbative expansion of $ d \tilde\hat \sigma^{\rm (NLLwC)}_{q \bar q^{\prime}} / dQ^2$ truncated at the same order of $(\alpha_{\rm s})$ as in the NLO result.
The Mellin moments of the parton distribution functions are defined in  the usual way as
$$
f^{(N)}_{i/h} (\mu_{F}^2) = \int_0^1 dx \, x^{N-1} f_{i/h}(x, \mu_{F}^2)\,,
$$
and the inverse Mellin transform is evaluated using a contour ${\sf C}$ in the complex $N$ space using the minimal prescription method developed in~\cite{Catani:1996yz}. The total cross sections are calculated by integrating the invariant mass distribution in Eq.~(\ref{hires}) over $Q^2$.

\section{Numerical results }
\label{s:results}
In this section we discuss the results for NLO+NLL resummation of soft gluon  correction for the $s$-channel $tH$ production at the LHC collision energy $\sqrt S$=13 TeV. The NLO cross section is calculated using the aMC@NLO code~\cite{Alwall:2014hca}. We use  $m_t = 173$ GeV and $m_H = 125$ GeV and employ PDF4LHC15 parton distribution sets~\cite{Butterworth:2015oua, Dulat:2015mca, Harland-Lang:2014zoa, Ball:2014uwa, Gao:2013bia, Carrazza:2015aoa}, apart from LO results which are calculated with the LO MMHT2014 set~\cite{Harland-Lang:2014zoa}. Unless otherwise stated, the results shown in the plots are obtained for the central scale choice $\mu_{R,0}=\mu_{F,0}=\mu_0=Q$ and for the jet radius $R=0.6$.

We begin by investigating how well an expansion of the resummed cross section approximates the NLO cross section. More precisely, we compare the expansion with the contribution to the NLO cross section originating from the $q \bar q'$ channel, i.e. the NLO cross section with the $qg$ channel removed, which we call "NLO (no qg)". Since the $qg$ channel is subleading in powers of $N$ w.r.t. the $ q \bar q'$ channel, it cannot be resummed within the current formalism and therefore should not be taken into consideration in the comparison. Given that the threshold definition we use here depends on the variable $Q$, in the first step we study differential distributions in $Q$ for a phenomenologically viable value of the $R$ parameter, $R=0.6$. In Fig.~\ref{fig1r06} we compare the NLO (no qg) distribution with the expansion of the NLL cross section up to NLO, called "NLLwC$|_{\rm NLO}$".
Additionally, in the massless case we also show distributions calculated without power suppressed terms in $N$, cf. Eq.~(\ref{eq:jet1massless}).
It is immediately clear from Fig.~\ref{fig1r06}(b) that the power suppressed corrections play a vital role. Therefore  from now on in all results where the final state jet  is treated as massless we will include these contributions. For most of the range of $Q$ shown in the plots, the expansion overestimates the NLO (no qg) distribution by around a few percent. Only at the lower end of the spectrum, the expansion underestimates the NLO result but this region is relatively small and does not provide a dominant contribution to the total cross section. We observe that the expanded results are quantitatively very similar independent of the assumption on the jet mass. It turns out that this is also a manifestation of the numerical importance of the power suppressed contributions. Namely, we find that the analytical expressions for  NLLwC$|_{\rm NLO}$ cross sections agree between the two approaches  in the limit $\delta^2 \ll 1/\bar N$, provided the power suppressed  contributions are included. Correspondingly, this behaviour also implies a relatively low relevance of large-$N$ terms for the considered process.

The lower panels of the plots in Fig.~\ref{fig1r06} display the ratio of the invariant mass distributions to the LO distribution. Firstly we notice that the NLO (no qg) corrections are substantial for all values of $Q$ and can be well above 50\% for the $Q$ values corresponding to the peak of the distribution. In comparison with the ratios for the NLO (no qg) distributions, the ratios of the NLLwC$|_{\rm NLO}$  to the LO distributions are flatter. This corresponds to the NLLwC$|_{\rm NLO}$ distributions closer in shape to the LO distributions than the NLO (no qg) distributions, in accordance with the former ones being calculated in the 3-particle kinematics limit.

\begin{figure}[t!]
\centering
\includegraphics[width=0.49\textwidth]{fig2ia_massive_PDF4LHC15_R06}
\includegraphics[width=0.49\textwidth]{fig2ib_massless_PDF4LHC15_R06}
\caption{Differential distributions in $Q$ for massive (left plot) and massless (right plot) jets  for $R=0.6$: (top) absolute values of the NLO(no qg), NLLwC|NLO distributions and (bottom) ratios of NLLwC|NLO and NLO(no qg) results to the LO predictions.}
\label{fig1r06}
\end{figure}

In Fig.~\ref{fig1r01} we study the invariant mass distributions obtained for a smaller value of $R$, $R=0.1$. For smaller $R$, the NLO (no qg) distribution in $Q$ naturally becomes narrower since narrower jets are less massive and the invariant mass of the full system decreases correspondingly. This can be also observed by comparing lower panels of Figs.~\ref{fig1r06} and~\ref{fig1r01}. Additionally, we see that the quality of the approximation of the NLO (no qg) result decreases slightly in comparison to the $R=0.6$ case. The dependence of the NLLwC$|_{\rm NLO}$ results on $\delta$ manifests itself in the higher absolute values for these distributions at smaller $R$. This dependence can be seen as an artifact of our setup, where by choosing the threshold variable dependent on the invariant mass $Q$, we are forced to define a jet, thus introducing the dependence on $\delta$ in the  NLLwC$|_{\rm NLO}$ cross sections. The dependence must cancel against $\delta$ dependence of the contributions from hard radiation, not included in the soft approximation, in the inclusive NLO cross section. At this very small value of $R$ some of the  the $\delta$-dependent terms can be very big individually, as can be seen in the right plot in Fig.~\ref{fig1r01} by comparing results with and without power subleading terms in $N$ in the massless jet case.

\begin{figure}[t!]
\centering
\includegraphics[width=0.49\textwidth]{fig1a_massive_PDF4LHC15_R01}
\includegraphics[width=0.49\textwidth]{fig1b_massless_PDF4LHC15_R01}
\caption{Differential distributions in $Q$ for massive (left plot) and massless (right plot) jets for $R=0.1$: (top) absolute values  of the NLO(no qg), NLLwC|NLO distributions and (bottom) ratios of NLLwC|NLO and NLO(no qg) results to the LO predictions).}
\label{fig1r01}
\end{figure}

In Fig.~\ref{fig2} we show scale variation for $\mu=\mu_F=\mu_R$ around the central scale $\mu_0=Q$ of a set of $K$ factors, defined as the ratios of the NLO, NLO (no qg) and NLLwC$|_{\rm NLO}$ total cross sections to the LO cross section. First we observe that the full NLO corrections to the total cross section are substantial, above 35\% within the range of $\mu$ considered here. Fig.~\ref{fig2} also shows that most of the NLO cross section is provided by the quark channel, while the NLO $qg$ terms at $\mu/\mu_0=1$ provide contributions of around 5\% to the full NLO result or of around 20\% to the NLO corrections. These contributions are negative for the most of the scale range shown, specifically  where the difference between NLO (no qg) and NLO is the most pronounced. As already discussed, in our approach resummation is performed for the dominant $q\bar q^{\prime}$ channel, while the $qg$ channel is included at the NLO accuracy.

The overestimation of the NLO (no qg) results by the  NLLwC$|_{\rm NLO}$  approximation for most values of $Q$ observed for differential distributions carries over to the total cross sections. Consistently with differential results, we also observe that the numerical values of the NLLwC$|_{\rm NLO}$ expansions agree very well between results for massless and massive jets. The expansions overestimate the NLO (no qg) result by around 3\% at $\mu=\mu_0=Q$. In terms of $K$-factors, this corresponds to $K_{\rm NLO\;(no\;qg)}$ of around $1.45$ at this value of $\mu$, while the expansions provide $K$-factors of around 1.5, see Fig.~\ref{fig2}. The agreement among the $K$-factors improves significantly as the scale increases, indicating that the scale-dependent terms provide considerable contributions to the NLO (no qg) cross section.

\begin{figure}
\centering
\includegraphics[width=0.55\textwidth]{fig3_mudep_Kfact_expansions}
\caption{Total $K$-factors as a function of $\mu/\mu_0$ for $R=0.6$: ratios of the total cross sections calculated at the NLO, NLO(no qg) and NLLwC|NLO accuracy to the LO total cross section with the final state jet treated as massive (M) or massless (m).}
\label{fig2}
\end{figure}

The dependence of the NLO expansions of the resummed cross section on the jet parameter $R$ is plotted in Fig.~\ref{fig3} for the same $K$-factors as shown in Fig.~\ref{fig2}. As $R$ gets smaller, the terms logarithmic in $\delta$ become more and more relevant, driving the difference between the expansion of the resummed cross section and the NLO result further apart. However, the dependence on $R$ is relatively mild for moderate $R$ values. Overall, for the central scale choices $\mu_0=Q$ and a wide range of $R$ values, the NLO(no qg) result differs from the NLLwC$|_{\rm NLO}$ expansion by up to a few percent.

\begin{figure}[t!]
\centering
\includegraphics[width=0.55\textwidth]{fig4_Rdep_Kfact_expansions}
\caption{Total $K$-factors as a function of $R$ at $\mu=\mu_0=Q$: ratios of the total cross sections calculated at the NLO, NLO(no qg) and NLLwC|NLO accuracy with massive (M) and massless (m) jets to the LO total cross section.}
\label{fig3}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[width=0.49\textwidth]{fig5a_massive_PDF4LHC15_R06}
\includegraphics[width=0.49\textwidth]{fig5b_massless_PDF4LHC15_R06}
\caption{Differential distributions in $Q$ for massive (left plot) and massless (right plot) jets for $R=0.6$: (top) absolute values of the NLO, NLO+NLLwC distributions and (bottom) ratio of the NLO+NLLwC result to the NLO prediction.}
\label{fig4}
\end{figure}

Next, we discuss the resummed NLLwC results matched to the NLO predictions. Fig.~\ref{fig4} shows differential distributions in $Q$ for massive and massless jets. We see that the corrections beyond NLO taken into account in the resummed formula amount to 1-4\%, depending on the value of $Q$ and the treatment of the jet. Inspecting the ratios of the NLO+NLLwC to NLO distributions, we also observe that the resummed corrections increase with growing $Q$, especially for massive jets. This behaviour is expected, as soft gluon corrections should be more pronounced closer to the threshold, approched with higher $Q$.

Fig.~\ref{fig5} shows the NLO+NLLwC and NLO total cross sections as functions of the scale ratio $\mu/\mu_0$ for $\mu_0=Q$.  We see that resummation provides a reduction of the scale dependence for both massless and massive jets, with a more noticable effect for massless jets. The higher-order NLLwC corrections modify the absolute value of the total cross section at $\mu=\mu_0=Q$  by up to 3\% percent. The value of these corrections depends on the treatment of the final state jet since $\tilde {J}_{out}^{(1)}$ functions and the resummed final state jet factors $J_k$ differ between the two approaches. Consequently, even if the ${\cal O}(\alpha_{\rm s})$  NLLwC$|_{\rm NLO} $ expanded results are very close  to each other numerically, cf. Fig.~\ref{fig2},  the resummed cross sections can differ visibly, as they take into account different terms starting from  ${\cal O}(\alpha_{\rm s}^2)$. For example, treating jets as massive leads to resummation of terms with single logarithms of $N$ in the final state jet function. They come together with a $\log \delta$ coefficient, effectively partially resumming also the $\log \delta$ terms. In the massless jet approach, final state jet function does not contain any dependence on $\delta$.

On the other hand, as demonstrated in Fig.~\ref{fig6},  if the scale variation of results obtained with various central scale choices is considered, the NLO+NLLwC results show much smaller spread in absolute values than the NLO results, even after accounting for the spread due to possible various treatments of the final state jet. It clearly indicates a great potential of resummation methods to increase the stability of the theoretical predictions and motivates further study of the subject.

\begin{figure}
\centering
\includegraphics[width=0.55\textwidth]{fig6_mudep_Q_R06_PDF4LHC15_matching_LO}
\caption{Total cross section as a function of $\mu/\mu_0$ for $R=0.6$ at the LO, NLO and the NLO+NLLwC accuracy with the final state jet treated as massive (M) or massless (m).}
\label{fig5}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{fig7a_mu_matching_scales_massive_PDF4LHC15}
\includegraphics[width=0.49\textwidth]{fig7b_mu_matching_scales_massless_PDF4LHC15}
\caption{Total cross section as a function of $\mu/\mu_0$ for $R=0.6$ at the NLO and NLO+NLLwC accuracy for massive (left plot) and massless (right plot) jets. Results are shown for three central scale choices $\mu_0=Q$ (solid lines), $\mu_0=H_T/2$ (dashed lines), and $\mu_0=H_T/6$ (dashed-dotted lines). }
\label{fig6}
\end{figure}

In our studies, we use a moderate value of the jet radius, $R=0.6$. It can be then expected that logarithmic terms in $R$, or equivalently $\delta$, should not provide large contributions to the cross section.
While the dependence of the matched resummed results on the assumption regarding the jet mass is not very big, and amounts to 3\% of the total NLO cross section at $\mu_0=Q$, cf. Fig.~\ref{fig5}, it is only slightly smaller than the overall scale uncertainty estimated for $\mu_R=\mu_F$. As we have noted above, the difference between the two results is driven by the $\delta$-dependent terms. Their importance is  heightened by the relative smallness of the pure soft gluon corrections arising in the invariant mass threshold limit for the $tH$ production.  This is particularly visible in the massless jet case, where, as we have discussed, in order to numerically approach the NLO result,  the corrections of ${\cal O}(1/N)$, originating from the region $\delta^2 \bar N \sim 1$, cannot be neglected. While we have included these terms in our calculations, in general they can only generate a part of power-subleading terms at higher order, what may be seen as a drawback of the massless jet approach. In contrast, the NLLwC formula we use for a massive final-state jet does not involve terms formally subleading in $N$. However, this approach does not treat the logarithmic terms in $\delta$ on the same footing, resumming them only partially. At the heart of these intricacies and differences between the two formulations lies a discord between the inclusive treatment of the jet, embodied in the corresponding factorization formula and the notion of the invariant mass threshold, which as we have seen involves some specification of the jet momentum.  While for processes with bigger contributions from soft gluon corrections, i.e. those involving gluons in the initial state, the approach presented here might lead to a very small dependence on the treatment of the jet, our studies show that the case of the $tH$ process requires much more care.

Apart from considering a different threshold variable in a more inclusive kinematics, one could also consider a more exclusive process, i.e. $t$, $H$ and jet in the final state. This would imply a different factorization formula, including non-global logs and separating soft-collinear modes from the global soft and collinear modes, which offers a possibility to sum logarithms of $R$, as demonstrated within the soft-collinear field theory framework~\cite{Liu:2017pbb, Becher:2016mmh, Larkoski:2015zka, Larkoski:2016zzc}. Such an analysis is, however,  beyond the scope of the present paper.

\section{Summary}

The NLO corrections to the $s$-channel $tH$ production are significant, at the level of 35\%. The main contribution to them comes from the quark-initiated production channel, which indicates the importance of higher-order corrections in this channel. The work presented here focuses on exploring ways in which resummation methods in direct QCD resummation can be used to provide an estimate of the higher-order corrections. Of course, as the $tH$ production takes place only in the quark-antiquark channel at LO, the resummation corrections cannot to be expected to be as important as for  the process with gluons in the initial state. However, it is a formally well defined set of higher-order corrections that can be systematically taken into account to all orders and studied numerically. The relevance of the measurement of the $tH$ process for gaining information on the top Yukawa coupling, as well as the increasing experimentall precision make such studies more and more important with time.

In this exploratory work, we investigated a method of threshold resummation in direct QCD, set in the invariant mass kinematics.  The definition of the threshold variable through the invariant mass of the final state system inevitably involves jet's momentum, leading to a residual dependence on a jet-size parameter $R$.  Our numerical results carry typical characteristic of the resummed predictions such as reduced scale dependence in comparison with the fixed-order results. However, it turns out that the pure soft gluon corrections at the invariant mass threshold, i.e. terms logarithmic in $N$, play a relatively small role for the process at hand. This is reflected in the fact that in the two cases we consider here, i.e. treating the final state jet as either massless or massive, quantitatively very similar $\cal{O}(\alpha_{\rm s})$ corrections  are obtained at different level of accuracy in powers of $N$. Naturally, in such a situation non-logarithmic terms, including also those  dependent on $R$, carry more relevance. In consequence, the resummed results in the two considered cases differ numerically from each other, as they take into account different subsets of higher-order terms.
To the best of our knowledge, though substantial numerical differences resulting from the treatment of the final state jet have been observed before, e.g. for the single-inclusive jet production~\cite{Kumar:2013hia, deFlorian:2013qia}, this is the first time where this difference is attributed, at least in part, to the terms subleading in powers of $N$. Consequently, our studies indicate the importance of the power-subleading terms in $N$,
underlining a general need for more developments aimed at their systematic all-order treatment. Similarly, our results also indicate that the calculations of the full higher-order corrections are very much needed in order to improve the description of the $tH$ production process.\\

\pagebreak
\noindent
{\Large \bf Acknowledgements}\\
This work has been supported in part by the Deutsche Forschungsgemeinschaft (DFG) grant KU3103/2. V.T. acknowledges funding from the European Union's Horizon 2020 research and innovation programme as part of the Marie Sk\l odowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104), while L.M.V. acknowledges support from the DFG Research Training Group ``GRK 2149: Strong and Weak Interactions - from Hadrons to Dark Matter''.

\bibliographystyle{JHEP}
\providecommand{\href}[2]{#2}\begingroup\raggedright