# Catalytic role of HI in the interstellar synthesis of complex organic
molecule

Shuming Yang  Laboratory for Relativistic Astrophysics, Department of Physics,
Guangxi University, Nanning 530004, China  Peng Xie  School of Chemistry and
Chemical Engineering, Guangxi University, Nanning 530004, China  Enwei Liang
Laboratory for Relativistic Astrophysics, Department of Physics, Guangxi
University, Nanning 530004, China  Zhao Wang<EMAIL_ADDRESS>Laboratory for
Relativistic Astrophysics, Department of Physics, Guangxi University, Nanning
530004, China

###### Abstract

Using quantum chemical calculations, we model the pathways for synthesizing
two purine nucleobases, adenine and guanine, in the gas-phase interstellar
environment, surrounded by neutral atomic hydrogen (HI). HI is found active in
facilitating a series of fundamental proton transfer processes of organic
synthesis, including bond formation, cyclization, dehydrogenation, and H
migration. The reactive potential barriers were significantly reduced in the
alternative pathways created by HI, leading to a remarkable increase in the
reaction rate. The presence of HI also lowered the reactive activation
temperature from 757.8 K to 131.5-147.0 K, indicating the thermodynamic
feasibility of these pathways in star-forming regions where some of the
reactants have been astronomically detected. Our findings suggest that HI may
serve as an effective catalyst for interstellar organic synthesis.

Astrochemistry, Astrobiology, ISM: molecules

##  I Introduction

Finding the synthesis pathways of interstellar complex organic molecules
(COMs) of biological interest is paramount for explaining the origin of the
elementary building blocks of life detected on meteorites, comets and
asteroids (Burton et al., 2012; Morse & Chan, 2019; Singh et al., 2013; Misra
et al., 2017). The formation of biomolecules through direct collisions between
abiotic molecules is generally challenging in the cold environment of the
interstellar medium (ISM). Photochemical and radiation chemical processes,
particularly on ice/dust grains, are commonly regarded as the principal steps
towards synthesizing interstellar COMs (Materese et al., 2020; Öberg, 2016;
Herbst & Van Dishoeck, 2009; Mifsud et al., 2021; Carrascosa et al., 2021).
While the ion-neutral or radical routes are the commonly considered mechanisms
for the gas-phase synthesis of interstellar COMs, catalyzed reactions are
emerging as alternative processes since the majority of the detected
interstellar molecules are neutral and have a closed shell.

Catalysts, often in the form of free radicals, ions, or metal clusters, are
frequently utilized in laboratory experiments of organic chemistry to convert
a challenging molecule-molecule reaction into a more straightforward molecule-
radical or radical-radical reaction (Qin et al., 2017; Zhao et al., 2020). For
the synthesis of interstellar organics, previous works suggest that ice/dust
grains, hydronium, and carbamic acid could function as catalysts (Potapov et
al., 2019; Roy et al., 2007; da Silva & de Araujo, 2017; Saladino et al.,
2013; Qi et al., 2018; Hanine et al., 2020; Mendoza et al., 2004; Rimola et
al., 2012; Vinogradoff et al., 2015). For example, radicals such as ·NH2 and
·CN are considered to participate the reactions for synthesizing nucleobases
from HCN, and are capable of reducing the barrier and may even render some
reactions barrierless (Jeilani et al., 2013, 2016). Nonetheless, the low
abundance of the previously suggested catalysts often renders the catalyzed
reaction seemingly impossible in ISM, as it requires the collision of three
bodies. Fortunately, this is not always the case. If the concentration of the
catalyst is in great excess of that of the reactants, the reaction can occur
just as easily as a two-body collision. Neutral atomic hydrogen (HI) is the
simplest and the most abundant free radical in the universe. It dominates the
chemical composition of the ISM, and is virtually ubiquitous in every
interstellar environment. In hot molecular cores, the concentration of HI is
about $7-10$ orders of magnitude higher than that of COMs Herbst & Van
Dishoeck (2009).

Although HI is a common reactant or product in simple interstellar reactions,
its catalytic activity in interstellar organic synthesis remains unclear in
most astrochemistry models (Bergantini et al., 2014; Carrascosa et al., 2020;
Chuang et al., 2016; Fedoseev et al., 2015, 2016; Miksch et al., 2021; Slate
et al., 2020; Singh et al., 2014; Misra et al., 2017). In comparison to its
molecular form, atomic hydrogen is rare on Earth and is difficult to prepare
and control in laboratory experiments. As a result, quantum chemical
calculations are expected to shed light on whether HI plays a crucial role in
the interstellar formation of COMs (Zamirri et al., 2019). Aiming at this
question, we conducted a new series of ab initio quantum chemical calculations
to explore the potential catalytic activities of HI in the formation pathways
of adenine (Ade) and guanine (Gua).

Ade and Gua are two purine nucleobases constituting the genetic code in the
nucleic acids. After having been detected in meteorites, their
extraterrestrial origin and prebiotic formation are a long-debated focus of
attention in attempts to answer the abiogenesis question (Burton et al., 2012;
Peeters et al., 2003). Intensive efforts are therefore devoted to find the
synthesis pathways of nucleobases in laboratory experiments simulating the
chemical evolution of abiotic molecules in interstellar environments (Sandford
et al., 2020; Chuang et al., 2016; Fedoseev et al., 2016; Parker & Kaiser,
2017; Merz Jr et al., 2014; Bera et al., 2017; Glaser et al., 2007; Choe,
2021; da Silva & de Araujo, 2017; Parker et al., 2015; Gupta et al., 2011;
Fulvio et al., 2021; Parker et al., 2015). In particular, Ade, Gua and purine
derivatives have been synthesized through photochemical reactions of purine in
an interstellar ice analogue (Materese et al., 2017, 2018). 2-aminopurine and
isoguanine are assumed to be two key intermediates for these reactions.
Nevertheless, the HI-catalyzed synthesis of nucleobases has yet to be
investigated.

##  II Methods

We consider novel synthesis pathways of Ade and Gua starting from
1h-pyrimidine-2-one C4H4N2O, a direct intermediate toward nucleobases
producted by two interstellar molecules: formamide H2NCHO and vinyl cyanide
H2CCHCN (Lu et al., 2021). Our proposed pathways involve two steps,
illustrated in Figure 1. In the 1st step (via the paths A-D), C4H4N2O reacts
with either cyanamide (NH2CN) or carbodiimide (HNCNH) (both detected in Sgr B2
and Orion-KL (Turner et al., 1975; McGuire et al., 2012; Belloche et al.,
2013; White et al., 2003) to form two direct purine-base precursors,
2-hydroxypurine (C5H4N4O) and 2-oxopurine (C5H4N4O), via the paths A-D. In the
2nd step, these two precursors produce Ade or Gua by separately reacting with
·NH2 via the paths E-H. We focus on the catalytic activities of HI in the both
two steps.

Figure 1: Formation scheme of Ade and Gua.

Multiple possible collision trajectories were simulated for each reaction
pathway. The transitional and final molecular structures in single molecular
collision events were determined by structural optimization using first
principles calculations (Maeda et al., 2015; Barone et al., 2015; Puzzarini,
2022; Zamirri et al., 2019; Pan et al., 2021), in which the system’s total
potential is minimized in the framework of the density functional theory (DFT)
at M06-2X/6-31+G(d,p) level as implemented in Gaussian 16 B.01 (Frisch et al.,
2016). We chose this level of theory in view of its proven records of accuracy
in predicting physical chemical properties of relevant compounds (Mardirossian
& Head-Gordon, 2017; de Jesus et al., 2021). In our previous work Lu et al.
(2021), this method was compared with calculations based on Møller Plesset
perturbation theory in the second order (MP2), and a good agreement was
achieved.

Vibration frequency calculations were here carried out to confirm the
intermediate (IM) and the transition state (TS). Intrinsic reaction coordinate
(IRC) calculations were performed to ensure that the TS reasonably connect the
reactant (RC) and the product (PD). To incorporate HI in our calculations, we
initially positioned a HI in a transition state with the reactant. We then
performed IRC coordinate scanning to ensure that the HI molecule’s placement
was appropriate. Subsequently, we simulated multiple trajectories of
collisions between the HI and the reactant from various positions. A static
calculation is used for this step, unlike what occurs in molecular dynamics
Wang (2020, 2019). Finally, we selected trajectories with the lowest energy
cost to facilitate the reaction. The atomic coordinate data of all optimized
molecular structures, including RC, IM, TS, and PD, are provided in the
standard XYZ chemical file format in the Supplementary Information’s Data
section.

The Gibbs free energy $G$ was calculated for every optimized structure of RC,
IM, TS or PD at a thermodynamic temperature of $150$ K under the assumption of
non-interacting particles (McQuarrie & Simon, 1999) (p321-325). This
representative temperature corresponds to a typical temperature of hot
molecular cores. The potential barrier $\Delta G$ of a reaction pathway was
determined as the difference in $G$ between the highest-energy states
(including IM, TS, or PD) and the RC in each of the two steps described in
Figure 1. This definition of the reactive barrier assumes that the energy does
not leave the system by irradiation on a timescales typical of the atomic
rearrangements, which is in general much faster than the needed vibration
cascade. As multiple collision trajectories were simulated for each path, only
the energetically favorable trajectories with the lowest $\Delta G$ are
discussed in the following.

##  III Results and discussions

Figure 2: Minimum reactive potential barrier for different reaction paths.

the reactive barrier $\Delta G$ is a crucial parameter for assessing the
feasibility and significance of a reaction in the ISM, as it strongly
correlates with the reaction rate and activation temperature. Figure 2
displays the minimum energy costs of producing Ade and Gua from C4H4N2O
through various paths, and compares $\Delta G$ for reactions with (red bars)
and without (black bars) the involvement of HI. The results demonstrate that
HI greatly promotes the reactions.

In the absence of HI, the reaction pathway (Path A-D, left half) shows a high
barrier with $\Delta G>83$ kcal/mol, which is expected since the fused
pyrimidine-imidazole rings in Ade and Gua are complex to synthesize without a
catalyst. However, the addition of HI significantly facilitates the reactions,
reducing $\Delta G$ to $11.6-19.4$ kcal/mol. Path D, involving HNCNH, shows
the lowest barrier among the proposed paths, indicating that it is more active
than its isomer NH2CN. Considering steps 1 and 2 consecutively, the optimal
pathway is either Paths B+G/D+G or D+H, with HI as a catalyst. The
corresponding optimal barrier for synthesizing Ade or Gua is around $15.8$ or
$14.1$ kcal/mol, respectively.

Figure 3: Reaction rate constant versus temperature for the optimal paths
without (black curve) and with HI (red and blue curves). The ternary reactions
with HI are assumed to be pseudo binary reactions.

As the temperature of a gas rises, its constituent molecules collide more
frequently and with greater kinetic energy. As a result, the fraction of
collisions that can surmount the barrier $\Delta G$ increases, following the
Arrhenius relation. In chemical kinetics, the Eyring-Polanyi equation (Eyring,
1935) describes the relationship between the reaction rate constant $k$,
$\Delta G$, and $T$,

$k=\frac{C\kappa_{B}T}{h}e^{-\frac{\Delta G}{RT}}.$ (1)

where $C$ is the transmission coefficient, $\kappa_{B}$ is the Boltzmann
constant, $R$ is the molar gas constant, and $h$ is the Planck constant. Using
this equation, we estimate $k$ as a function of $T$ as plotted in Figure 3,
for the optimal reactive paths of Ade and Gua with and without HI. Generally
speaking, most significant reactions in the ISM typically proceed with a rate
constant of 10-11 cm3s-1 or larger (above the dashed horizontal orange line).
Without catalyst, such a rate constant would roughly require a reactive
temperature $>757.8$ K with $\Delta G=83.9$ kcal/mol, as indicated by the
dashed line intersecting with the black curve in Figure 3. The feasibility of
such a reaction would be low in the generally cold, gas-phase environment of
dense molecular clouds, as the reactants would have to coincide with radiation
having an appropriate amount of energy at the time of collision. In contrast,
the presence of HI significantly accelerates the reaction, leading to a
substantial increase in the reaction rate constant by approximately $20-30$
orders of magnitude (as seen in Figure 3). As a result, the required
activation temperature is respectively reduced to only $131.5$ or $147.0$ K
for Gua or Ade, as indicated by the dashed lines intersecting with the red and
blue curves. To verify our computational model, we compute the half-life time
based on the reaction rate and compare it to the results of Zamirri et al.
(2019), with agreement obtained as shown in Part I of the Supplementary
Information.

Three-body collision processes are often deemed impossible in the ISM, but
this does not hold true for the HI-catalyzed reaction. This is because the
concentration of HI vastly exceeds that of COMs, reaching 107-1010 times in
hot cores Herbst & Van Dishoeck (2009). Thus, ternary reactions with HI are
essentially pseudo-binary reactions. In chemistry models for star-forming
regions, the HI abundance is estimated to be comparable to that of H2 (Garrod
et al., 2008; Sun & Du, 2022). In addition to the remaining background HI from
large-scale molecular clouds, HI also stems from the photo-dissociation of H2
by the radiation emitted from YSOs. This in particular points to the inner
parts of protoplanetary disks, in which the ultraviolet radiation can be
dramatically enhanced and resulting in a chemical structure resembling that of
a photon-dominated region (PDR) (van Dishoeck et al., 2006). In such a region,
photons are usually not energetic enough to ionize hydrogen but can dissociate
molecular hydrogen (Hollenbach & Tielens, 1999). For example, the total HI
mass in a protoplanetary disk is estimated to be $11.4$ $M_{Earth}$ with a
typical vertical column density of $1.7\times 10^{21}$ or $3.3\times 10^{20}$
cm-2 at $5$ or $110$ AU, respectively (Kamp et al., 2008). Moreover, HI is a
common product of many important reactions in ISM (Herbst & Klemperer, 1973),
as proton transfer is a fundamental process in organic reaction (Mayer, 2011).

The involvement of HI makes the proposed pathways thermally feasible in
certain interstellar environments with low activation temperatures. For
instance, the formation sites of young stellar objects (YSO) are often
associated with warm ambient gas, known as hot cores or corinos, with average
temperatures up to $300$ K, where a rich organic chemistry is observed (Herbst
& Van Dishoeck, 2009; Jørgensen et al., 2020). Interestingly, most of here-
used reactants, including formamide (NH2CHO), vinyl cyanide (CH2CHCN), NH2CN,
and HNCNH, are detected to co-exist in SgrB2 and Orion KL hot cores (Turner et
al., 1975; McGuire et al., 2012; Belloche et al., 2013; White et al., 2003;
Scibelli et al., 2021). Note that the reactant C4H4N2O is supposed to be
synthesized from NH2CHO and CH2CHCN (Lu et al., 2021), which are also both
detected in those regions as “hot” molecules with characteristic rotation
temperature $>100$ K (Lopez-Sepulcre et al., 2019; Belloche et al., 2013;
López et al., 2014; Suzuki et al., 2018; Fu & Lin, 2016). Those regions
contain the initial materials as well as enough thermal energy for the
proposed reactions to proceed. It is noteworthy that protoplanetary disks are
also of particular interest for abiogenesis, as the planets formed in them
inherit their chemical compositions. Although only simple organics have been
detected in protoplanetary disks so far (Öberg & Bergin, 2021), the activation
temperature reported in this study for nucleobase synthesis is lower than that
in the inner part of disks, as the temperature gradually increases from $10$ K
to several thousand K from the outer layer to the central forming star
(Akimkin et al., 2013).

Figure 4: Potential energy surface for forming 2-hydroxypurine (C5H4N4O) from
C4H4N2O and NH2CN in via Path A without (upper panel) or with HI (lower
panel). The arrows indicate the direction of molecular collision.

One might naturally inquire about the specific effects of HI in promoting the
reaction. To showcase the catalytic capabilities of HI, Figure 4 provides a
comparison of the reaction pathways in the absence and presence of HI, using
path A as an example. HI is found to facilitate several critical processes,
including cyclization, dehydrogenation, and the bonding of the reactants. In
the absence of HI, the reaction is impeded by a significant barrier of $86.8$
kcal/mol against the formation of an initial bond between the two reactants
(see the step in the red line, upper panel). However, when HI is present in
the vicinity, it can extract a hydrogen atom from the amino group with a
relatively low cost of only $14.4$ kcal/mol, converting NH2CN into the highly-
reactive HNCN radical. This generates an alternative pathway in which the
barrier for bonding the two reactants is considerably lowered to $10.8$
kcal/mol, as indicated by the red line in the lower panel.

In path A, the rate determining step involves two conventional organic
synthesis processes, namely cyclization and dehydrogenation. Initially,
cyclization exhibits a high barrier of around $49.0$ kcal/mol, as indicated by
the blue circle in the upper panel of Figure 4. This barrier arises because
the N end required for forming the imidazole ring is completely saturated.
However, this process is expedited in the alternative pathway created by HI,
where the (1h-pyrimidine-2-one-6-yl)-cyanamide radical (IM21) possesses an
unsaturated N end, as depicted by the blue rectangle in the lower panel of
Figure 4. Although cyclization still requires approximately $30.6$ kcal/mol,
the energy retained in the preceding sub-steps helps to surmount this barrier,
owing to the low energy of IM21.

The next step is dehydrogenation, in which a C-bonded and an N-bonded hydrogen
atom in the imidazole ring are converted into a H2 molecule, resulting in the
formation of 2-hydroxypurine C5H4N4O, as shown in Figure 4. The original
barrier for this process is roughly $53.2$ kcal/mol, which is substantially
reduced in the alternative pathway created by HI in the lower panel (green
line). One of the H atoms in IM23 is taken by the free HI, at a small cost of
only $5.7$ kcal/mol, to generate the stable IM24, 1,6-dihydro-purin-2-one
radical. This radical ultimately yields 2-hydroxypurine (C5H4N4O) by directly
removing another H atom.

Figure 5: HOMO (left) and LUMO (right) of C4H4N2O (a), NH2CN (b) ,HNCN (c).

The catalytic activity of HI is demonstrated using the initial step of Path A
in Figure 4 as an example within the framework of the frontier molecular
orbital theory. This step involves the formation of a bond between the
reactants, C4H4N2O and NH2CN, to initiate the reaction. However, because the
highest occupied molecular orbital (HOMO) of C4H4N2O (located on the left in
Figure 5 (a)) and the lowest unoccupied molecular orbital (LUMO) of NH2CN
(outlined in red) have mismatched symmetry ($\pi$ vs $\sigma$ orbital), a high
energy is required to form the bond. The feasibility of a reaction process is
closely related to the conservation of orbital symmetry, and a lack of
symmetry can impede the reaction (Woodward & Hoffmann, 1965, 1969).
Consequently, a significant HOMO-LUMO gap of $7.1$ kcal/mol arises, making the
reaction difficult to proceed. Fortunately, HI has the ability to extract a
hydrogen atom from NH2CN, thereby transforming it into a radical HNCN, which
can readily initiate the reaction because its LUMO symmetry (outlined in
green) is consistent with the HOMO symmetry of C4H4N2O ($\pi$ vs $\pi$
orbital). As a result, the HOMO-LUMO gap is significantly reduced to $3.6$
kcal/mol, allowing for the reaction to proceed smoothly.

Figure 6: Potential energy surface to form Gua from 2-oxopurine (C5H4N4O) and
·NH2 via Path H without (upper) or with HI (lower). The processes in the two
paths become the same after IM152 and IM163.

HI has also been found to aid in hydrogen migration, as demonstrated in path
H, where 2-oxopurine (C5H4N4O) reacts with ·NH2 to form Gua, as depicted in
Figure 6. To initiate this reaction, a proton must move from an N site to an
adjacent O site, which requires approximately $31.0$ kcal/mol (see yellow line
in the upper panel). This barrier can be reduced by the presence of HI, which
binds to the relevant O atom, creating a highly stable intermediate, IM161,
2-hydroxy-3H-Purine radical, as shown in the lower panel. This intermediate
sits in a deep potential well, and the energy it contains aids in overcoming
the subsequent processes of dehydrogenation and bonding between ·NH2 and a
$sp^{2}$-hybridized C.

The catalytic effects of HI are observed in all of the simulated reaction
paths (B-G), as provided in Part II of the Supplementary Information.
Additionally, several stable intermediates are identified in the proposed
reaction pathways, including 2-aminopurine, purinol, 2-hydroxypurine,
2-oxopurine, and isoguanine. Some of these intermediates align with the
findings of prior laboratory experiments. For instance, 2-aminopurine (IM153
or IM164 in Figure 6) was discovered by Materese et al. (2017, 2018) in a
photochemistry experiment on an interstellar ice analogue and was believed to
be a critical intermediate in the synthesis of Ade, Gua, and purine
derivatives. The results of this experiment also suggest the same function for
isoguanine and purinol. The present study elucidates the pathways from these
intermediates to the purine bases. Finally, the IR and electronic adsorption
spectra of these intermediates are calculated using DFT calculations (Kovács
et al., 2020; Meng et al., 2021; Wu et al., 2022) and are presented in Part
III of the Supplementary Information to facilitate potential astronomical
observations.

##  IV Conclusions

The proposed new pathways for the synthesis of Ade and Gua demonstrate the
efficiency of HI as a catalyst in promoting various fundamental organic
synthesis processes, including bond formation, cyclization, dehydrogenation,
and H migration. These are typical proton transfer processes in organic
synthesis, where HI can effectively lower the potential barrier and
significantly increase the reaction rate constant by up to $20-30$ orders of
magnitude, depending on the temperature. As a result, the activation
temperature required for the proposed reaction is lowered from 757.8 K to
131.5-147.0 K when HI is present. This suggests that the thermal feasibility
of gas-phase purine base synthesis with HI could be possible in environments
such as hot molecular cores or protoplanetary disks. From a broader
perspective, given the prevalence of HI as a major component of interstellar
gas and the fundamental role of proton transfer in organic synthesis, the HI-
catalyst mechanism should have wide implications for the formation of organic
molecules in the ISM.

## References

  * Akimkin et al. (2013) Akimkin, V., Zhukovska, S., Wiebe, D., et al. 2013, ApJ, 766, 8
  * Barone et al. (2015) Barone, V., Biczysko, M., & Puzzarini, C. 2015, AcChR, 48, 1413
  * Belloche et al. (2013) Belloche, A., Müller, H. S., Menten, K. M., Schilke, P., & Comito, C. 2013, A&A, 559, A47
  * Bera et al. (2017) Bera, P. P., Stein, T., Head-Gordon, M., & Lee, T. J. 2017, AsBio, 17, 771
  * Bergantini et al. (2014) Bergantini, A., Pilling, S., Rothard, H., Boduch, P., & Andrade, D. 2014, MNRAS, 437, 2720
  * Burton et al. (2012) Burton, A. S., Stern, J. C., Elsila, J. E., Glavin, D. P., & Dworkin, J. P. 2012, Chem Soc Rev, 41, 5459
  * Carrascosa et al. (2020) Carrascosa, H., Cruz-Díaz, G., Muñoz Caro, G. M., Dartois, E., & Chen, Y. 2020, MNRAS, 493, 821
  * Carrascosa et al. (2021) Carrascosa, H., González Díaz, C., Muñoz Caro, G., Gómez, P., & Sanz, M. L. 2021, MNRAS, 506, 791
  * Choe (2021) Choe, J. C. 2021, ApJ, 914, 136
  * Chuang et al. (2016) Chuang, K.-J., Fedoseev, G., Ioppolo, S., Van Dishoeck, E., & Linnartz, H. 2016, MNRAS, 455, 1702
  * da Silva & de Araujo (2017) da Silva, J. B. P., & de Araujo, A. P. M. 2017, ESC, 1, 376
  * de Jesus et al. (2021) de Jesus, D. N., da Silva, J. M., Tejero, T. N., et al. 2021, MNRAS, 501, 1202
  * Eyring (1935) Eyring, H. 1935, JChPh, 3, 107
  * Fedoseev et al. (2016) Fedoseev, G., Chuang, K.-J., van Dishoeck, E. F., Ioppolo, S., & Linnartz, H. 2016, MNRAS, 460, 4297
  * Fedoseev et al. (2015) Fedoseev, G., Cuppen, H. M., Ioppolo, S., Lamberts, T., & Linnartz, H. 2015, MNRAS, 448, 1288
  * Frisch et al. (2016) Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 2016, Gaussian 16 vB.01 (Wallingford, CT: Gaussian Inc)
  * Fu & Lin (2016) Fu, L., & Lin, G.-M. 2016, RAA, 16, 182
  * Fulvio et al. (2021) Fulvio, D., Potapov, A., He, J., & Henning, T. 2021, Life, 11, 568
  * Garrod et al. (2008) Garrod, R. T., Weaver, S. L. W., & Herbst, E. 2008, ApJ, 682, 283
  * Glaser et al. (2007) Glaser, R., Hodgen, B., Farrelly, D., & McKee, E. 2007, AsBio, 7, 455
  * Gupta et al. (2011) Gupta, V., Tandon, P., Rawat, P., Singh, R., & Singh, A. 2011, A&A, 528, A129
  * Hanine et al. (2020) Hanine, M., Meng, Z., Lu, S., et al. 2020, ApJ, 900, 188
  * Herbst & Klemperer (1973) Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505
  * Herbst & Van Dishoeck (2009) Herbst, E., & Van Dishoeck, E. F. 2009, ARA&A, 47, 427
  * Hollenbach & Tielens (1999) Hollenbach, D. J., & Tielens, A. 1999, RvMP, 71, 173
  * Jeilani et al. (2013) Jeilani, Y. A., Nguyen, H. T., Newallo, D., Dimandja, J.-M. D., & Nguyen, M. T. 2013, PCCP, 15, 21084
  * Jeilani et al. (2016) Jeilani, Y. A., Williams, P. N., Walton, S., & Nguyen, M. T. 2016, PCCP, 18, 20177
  * Jørgensen et al. (2020) Jørgensen, J. K., Belloche, A., & Garrod, R. T. 2020, ARA&A, 58, 727
  * Kamp et al. (2008) Kamp, I., Freudling, W., Robberto, M., Chengalur, J., & Keto, E. 2008, PhyS, 2008, 014013
  * Kovács et al. (2020) Kovács, P., Zhu, X., Carrete, J., Madsen, G. K., & Wang, Z. 2020, ApJ, 902, 100
  * Lopez-Sepulcre et al. (2019) Lopez-Sepulcre, A., Balucani, N., Ceccarelli, C., et al. 2019, ESC, 3, 2122
  * López et al. (2014) López, A., Tercero1, B., Kisiel, Z., et al. 2014, A&A, 572, A44
  * Lu et al. (2021) Lu, S., Meng, Z., Xie, P., Liang, E., & Wang, Z. 2021, A&A, 656, A84
  * Maeda et al. (2015) Maeda, S., Harabuchi, Y., Ono, Y., Taketsugu, T., & Morokuma, K. 2015, IJQC, 115, 258
  * Mardirossian & Head-Gordon (2017) Mardirossian, N., & Head-Gordon, M. 2017, MolPh, 115, 2315
  * Materese et al. (2020) Materese, C. K., Gerakines, P. A., & Hudson, R. L. 2020, AcChR, 54, 280
  * Materese et al. (2018) Materese, C. K., Nuevo, M., McDowell, B. L., Buffo, C. E., & Sandford, S. A. 2018, ApJ, 864, 44
  * Materese et al. (2017) Materese, C. K., Nuevo, M., & Sandford, S. A. 2017, AsBio, 17, 761
  * Mayer (2011) Mayer, J. M. 2011, AcChR, 44, 36
  * McGuire et al. (2012) McGuire, B. A., Loomis, R. A., Charness, C. M., et al. 2012, ApJL, 758, L33
  * McQuarrie & Simon (1999) McQuarrie, D. A., & Simon, J. D. 1999, Molecular thermodynamics (Mill Valley, CA: Univ. Science Books) (Sterling Publishing Company)
  * Mendoza et al. (2004) Mendoza, C., Ruette, F., Martorell, G., & Rodriguez, L. 2004, ApJ, 601, L59
  * Meng et al. (2021) Meng, Z., Zhu, X., Kovács, P., Liang, E., & Wang, Z. 2021, ApJ, 922, 101
  * Merz Jr et al. (2014) Merz Jr, K. M., Aguiar, E. C., & da Silva, J. B. P. 2014, JPCA, 118, 3637
  * Mifsud et al. (2021) Mifsud, D. V., Hailey, P. A., Muiña, A. T., et al. 2021, FrASS, 8, 757619
  * Miksch et al. (2021) Miksch, A. M., Riffelt, A., Oliveira, R., Kästner, J., & Molpeceres, G. 2021, MNRAS, 505, 3157
  * Misra et al. (2017) Misra, A., Tandon, P., et al. 2017, RAA, 17, 1
  * Morse & Chan (2019) Morse, A. D., & Chan, Q. H. 2019, ESC, 3, 1773
  * Öberg (2016) Öberg, K. I. 2016, ChRv, 116, 9631–9663
  * Öberg & Bergin (2021) Öberg, K. I., & Bergin, E. A. 2021, PhR, 893, 1
  * Pan et al. (2021) Pan, L., Carrete, J., & Wang, Z. 2021, JPCM, 34, 015303
  * Parker & Kaiser (2017) Parker, D. S., & Kaiser, R. I. 2017, ChSRv, 46, 452
  * Parker et al. (2015) Parker, D. S., Kaiser, R. I., Kostko, O., et al. 2015, ApJ, 803, 53
  * Peeters et al. (2003) Peeters, Z., Botta, O., Charnley, S., Ruiterkamp, R., & Ehrenfreund, P. 2003, ApJ, 593, L129
  * Potapov et al. (2019) Potapov, A., Theulé, P., Jäger, C., & Henning, T. 2019, ApJL, 878, L20
  * Puzzarini (2022) Puzzarini, C. 2022, FrASS, 8, 252
  * Qi et al. (2018) Qi, H., Picaud, S., Devel, M., Liang, E., & Wang, Z. 2018, ApJ, 867, 133
  * Qin et al. (2017) Qin, Y., Zhu, L., & Luo, S. 2017, ChRv, 117, 9433
  * Rimola et al. (2012) Rimola, A., Sodupe, M., & Ugliengo, P. 2012, ApJ, 754, 24
  * Roy et al. (2007) Roy, D., Najafian, K., & von Ragué Schleyer, P. 2007, PNAS, 104, 17272
  * Saladino et al. (2013) Saladino, R., Botta, G., Delfino, M., & Di Mauro, E. 2013, ChEurJ, 19, 16916
  * Sandford et al. (2020) Sandford, S. A., Nuevo, M., Bera, P. P., & Lee, T. J. 2020, ChRv, 120, 4616
  * Scibelli et al. (2021) Scibelli, S., Shirley, Y., Vasyunin, A., & Launhardt, R. 2021, MNRAS, 504, 5754
  * Singh et al. (2013) Singh, A., Misra, A., Tandon, P., et al. 2013, RAA, 13, 912
  * Singh et al. (2014) Singh, A., Misra, A., Tandon, P., et al. 2014, RAA, 14, 275
  * Slate et al. (2020) Slate, E. C., Barker, R., Euesden, R. T., Revels, M. R., & Meijer, A. J. 2020, MNRAS, 497, 5413
  * Sun & Du (2022) Sun, J., & Du, F. 2022, RAA, 22, 065022
  * Suzuki et al. (2018) Suzuki, T., Ohishi, M., Saito, M., et al. 2018, ApJS, 237, 3
  * Turner et al. (1975) Turner, B., Liszt, H., Kaifu, N., & Kisliakov, A. 1975, ApJ, 201, L149
  * van Dishoeck et al. (2006) van Dishoeck, E. F., Jonkheid, B., & van Hemert, M. C. 2006, FaDi, 133, 231
  * Vinogradoff et al. (2015) Vinogradoff, V., Duvernay, F., Fray, N., et al. 2015, ApJL, 809, L18
  * Wang (2019) Wang, Z. 2019, JPCC, 123, 15166
  * Wang (2020) Wang, Z. 2020, JPCC, 124, 3851
  * White et al. (2003) White, G. J., Araki, M., Greaves, J., Ohishi, M., & Higginbottom, N. 2003, A&A, 407, 589
  * Woodward & Hoffmann (1965) Woodward, R. B., & Hoffmann, R. 1965, JAChS, 87, 395
  * Woodward & Hoffmann (1969) Woodward, R. B., & Hoffmann, R. 1969, AngCh Int, 8, 781
  * Wu et al. (2022) Wu, R.-Q., Long, M.-Q., Zhang, X.-J., et al. 2022, RAA, 22, 035014
  * Zamirri et al. (2019) Zamirri, L., Ugliengo, P., Ceccarelli, C., & Rimola, A. 2019, ESC, 3, 1499
  * Zhao et al. (2020) Zhao, L., Kaiser, R. I., Lu, W., et al. 2020, AngChIE, 59, 11334


Gluon Fusion Production at NLO: Merging the Transverse Momentum and the High-
Energy Expansions

Luigi Bellafrontea***email<EMAIL_ADDRESS>Giuseppe
Degrassib†††email<EMAIL_ADDRESS>Pier Paolo
Giardinoa‡‡‡email<EMAIL_ADDRESS>Ramona Gröberc§§§email:
<EMAIL_ADDRESS>Marco Vittib¶¶¶email<EMAIL_ADDRESS>
(a) Instituto Galego de Física de Altas Enerxías, Universidade de Santiago de
Compostela, 15782 Santiago de Compostela, Galicia-Spain
(b) Dipartimento di Matematica e Fisica, Università di Roma Tre and
INFN, sezione di Roma Tre, I-00146 Rome, Italy
(c) Dipartimento di Fisica e Astronomia ’G. Galilei’, Università di Padova and
INFN, sezione di Padova, I-35131 Padova, Italy

The virtual corrections to $gg\to HH$ and $gg\to ZH$ are analytically
evaluated combining an expansion in the small transverse momentum of the final
particles with an expansion valid at high energies. The two expansion methods
describe complementary regions of the phase space and we merge their results,
extending the range of validity of both expansions using Padé approximants. We
show that this approach can reproduce the available numerical results
retaining the exact top quark mass dependence with an accuracy well below the
1% level. Our results allow a fast and flexible evaluation of the virtual
corrections of the considered processes. Furthermore, they are available in
different renormalisation schemes of the top quark mass.

##  1 Introduction

At the Large Hadron Collider (LHC) gluon fusion is the most relevant
production mechanism for Higgs physics in single-Higgs production, $gg\to H$,
and pair production, $gg\to HH$, and plays an important role in the production
of a Higgs boson in association with a $Z$ boson, $gg\to ZH$. Precise
predictions for such processes are necessary in order to measure the
properties of the Higgs boson accurately. Since the Higgs boson is a colorless
particle, these processes are all loop-induced and mediated by a heavy
particle, mainly the top quark. As known since more than twenty years, QCD
corrections to the Born result are very large and the perturbative expansion
converges slowly (see for example [1, 2]). This implies that the knowledge of
the higher-order QCD corrections is very important. However, the evaluation of
these corrections is extremely challenging already at next-to-leading order
(NLO), since it requires the computation of two-loop diagrams.

In this paper we are going to consider the virtual corrections to Higgs
production via gluon fusion at the NLO level. In general, the degree of
difficulty in the evaluation of loop diagrams grows with the number of energy
scales present in the diagram. In the case of single-Higgs production the
relevant diagrams feature a triangular topology and, consequently, depend upon
only two scales, namely the Higgs mass, $m_{\scriptscriptstyle H}$, and the
top mass111 All the quarks but the top are assumed to be massless., $m_{t}$.
In this case, the functional dependence of the result upon the top mass can be
expressed in terms of one single variable, $m_{\scriptscriptstyle
H}^{2}/m_{t}^{2}$. Due to this “simplified” one-scale situation, exact
analytic results for the NLO corrections are available since many years [3, 4,
5, 6].

In the case of processes with two particles in the final state the situation
is more complicated. Indeed these processes receive contributions not only
from triangle diagrams, that can be calculated adapting the exact analytic
results obtained for single-Higgs production, but also from box-topology
diagrams. In pair production, $gg\to HH$, the box diagrams depend upon four
scales, namely $\hat{s},\,\hat{t},\,m_{t},m_{\scriptscriptstyle H}$, where
$\hat{s},\,\hat{t}$, and $\hat{u}$ are the Mandelstam variables which satisfy
the condition

$\hat{s}+\hat{t}+\hat{u}=2\,m_{\scriptscriptstyle H}^{2}~{}.$ (1)

Concerning associated production, $gg\to ZH$, a fifth energy scale is present,
i.e. the mass of the $Z$ vector boson, $m_{Z}$.

Exact analytic results for two-loop box diagrams with several energy scales
cannot be derived with the present computational technology. Instead, usually
two different strategies are followed in order to evaluate the two-loop box
contribution in Higgs production via gluon fusion. a) A fully numerical exact
evaluation [7, 8, 9, 10, 11]. b) An approximate analytic evaluation that takes
advantage of hierarchies among the various energy scales present in the
diagrams, in order to reduce the number of scales in the problem. Thus, its
validity is restricted to specific regions of the phase space. The method used
is based on the expansion of the diagrams in terms of ratios of small energy
scales vs. large energy scales, in order to obtain a result that retains an
exact dependence upon the large energy scales. Concerning the small ones, in
order to simplify further the evaluation, expansions in term of ratios between
small energy scales is often used.

The former strategy, although accurate, is very demanding from a computational
point of view, requiring a high degree of optimization in order to obtain a
result in a “reasonable”, although usually quite long, computer time.
Furthermore, this approach is not very flexible with respect to the
modification of the input parameters.

Strategy b) provides accurate results valid in specific regions of the phase
space without requiring heavy computational work, i.e. in a short computer
time. Examples of this approach of evaluating the two-loop box contribution
are:

  1. i)

The infinite-top-mass limit [12, 13] refined by the inclusion of powers in the
large top-mass expansion (LME) [14, 15, 16, 17]. Here, $m_{t}$ is assumed to
be the large energy scale while $\hat{s},\>\hat{t},\>m_{\scriptscriptstyle
H}$, and in associated production also $m_{\scriptscriptstyle Z}$, are
considered to be the small ones. Thus, the validity of this approach is
restricted to phase-space regions where $\hat{s}/(4m_{t}^{2})\leq 1$.

  2. ii)

The evaluation via an expansion in the transverse momentun,
$p_{\scriptscriptstyle T}$, of the final-state particles [18, 19]. Here,
$\hat{s}$ and $m_{t}$ are assumed to be the large energy scale while
$m_{\scriptscriptstyle H},\>m_{\scriptscriptstyle Z}$ and
$p_{\scriptscriptstyle T}$, that can be traded for $\hat{t}$, are considered
to be the small ones. The validity of this approach is restricted to phase-
space regions where $|\hat{t}|/(4m_{t}^{2})\lesssim 1$.

  3. iii)

The evaluation via a high-energy (HE) expansion [20, 21, 22]. Here
$\hat{s},\>\hat{t}$ are assumed to be the large energy scale while
$m_{t},\>m_{\scriptscriptstyle H}$ and $m_{\scriptscriptstyle Z}$, with
$m_{t}\gg m_{\scriptscriptstyle H},\,m_{\scriptscriptstyle Z}$, are considered
to be the small ones. The validity of this approach is restricted to phase-
space regions where $|\hat{t}|/(4m_{t}^{2})\gtrsim 1$.

  4. iv)

The evaluation via an expansion in terms of small external masses [23, 24,
25]. Here $\hat{s},\>\hat{t},\>m_{t}$ are assumed to be the large energy scale
while $m_{\scriptscriptstyle H}$ and $m_{\scriptscriptstyle Z}$ are considered
to be the small ones. This approach basically covers the entire phase space of
the considered processes. However, since the reduction of scales in this
approach is minimal, one ends up with the evaluation of Master Integrals (MIs)
that are much more complicated than those appearing in the i)–iii) cases. As a
consequence the evaluation of the box contribution in any point of the phase
space requires a longer computer time than in the approaches i)–iii).

As an alternative approach, refs. [26, 27] proposed to reconstruct the full
result from its LME version, supplemented by the non-analytic part of the
diagrams near the top threshold, via a conformal mapping and Padé
approximants.

In this paper we propose an alternative way to derive the full top-mass
dependence in Higgs production via gluon fusion, based on the merging of the
$p_{\scriptscriptstyle T}$ expansion in ii) with the HE expansion in iii) that
individually are valid in complementary regions of the phase space. Since the
numerical evaluations of the two expansions are quite fast from a
computational point of view, our proposal allows a fast evaluation of the
virtual corrections to Higgs production via gluon fusion that is accurate in
the entire phase space.

The key point of our analysis is to extend the fixed-order results both in the
$p_{\scriptscriptstyle T}$ expansion [18, 19] and in the HE expansion [21, 22]
up to or beyond their border of validity, i.e. $\hat{t}\simeq 4m_{t}^{2}$, in
order to merge the two analytic approximations. This is done by constructing a
[1/1] Padé approximant for the $p_{\scriptscriptstyle T}$-result and a [6/6]
Padé approximant for the HE-result. We point out that the extension of the HE
expansion via Padé approximants has been already considered in refs. [21, 22].

The paper is organized as follows. In section 2 we introduce the different
expansions as well as the method for combining them. In section 3 we validate
the method at LO both at the level of form factors and at the level of the
partonic cross section, with a focus on the $gg\to HH$ process. In section 4
we present the merging of the two expansions at NLO. In order to show the
flexibility of our approach with respect to the modification of the input
parameters, in the same section we present the result for the two-loop virtual
contribution in $gg\to HH$ and $gg\to ZH$ for two different choices of the top
quark mass, namely the on–shell and the $\overline{\text{MS}}$ mass. Finally,
we conclude in section 5.

##  2 Method

We start by considering the process $g(p_{1})\,g(p_{2})\rightarrow
3(p_{3})\,4(p_{4})$, where 3 and 4 are two neutral222Dealing with neutral
particles in the final state implies the absence of mixed top-bottom diagrams,
as they would appear e.g. in $gg\to W^{+}W^{-}$ production, for which the
proposed method cannot be straightforwardly applied. particles with masses
$m_{3}$ and $m_{4}$, respectively. Taking all momenta to be incoming, the
partonic Mandelstam variables are

$\hat{s}=(p_{1}+p_{2})^{2},\qquad\hat{t}=(p_{1}+p_{3})^{2},\qquad\hat{u}=(p_{2}+p_{3})^{2},$
(2)

and the transverse momentum $p_{\scriptscriptstyle T}$ of the final-state
particles can be written as

$p_{\scriptscriptstyle
T}^{2}=\frac{\hat{t}\hat{u}-m_{3}^{2}m_{4}^{2}}{\hat{s}}.$ (3)

As suggested in refs. [18, 19], if the amplitude of the process is written in
terms of (anti)symmetric form factors with respect to the exchange
$\hat{t}\leftrightarrow\hat{u}$, then it is sufficient to discuss only the
forward contribution to the cross section. Therefore, in the following we will
always assume that $|\hat{t}|\leq|\hat{u}|$ and that

$\hat{t}=-\frac{1}{2}\left(\hat{s}-m_{3}^{2}-m_{4}^{2}-\sqrt{\lambda(\hat{s},m_{3}^{2},m_{4}^{2})-4\hat{s}\,p_{\scriptscriptstyle
T}^{2}}\right),$ (4)

where $\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$ is the Källén function.

In the forward regime, the validity of both the $p_{\scriptscriptstyle T}$ and
HE expansions is limited by the condition

$|\hat{t}|\simeq 4m_{t}^{2},$ (5)

i.e. for any fixed value of $\hat{s}$, the $p_{\scriptscriptstyle T}$
expansion provides reliable results when $|\hat{t}|\lesssim 4m_{t}^{2}$ while
the HE expansion is accurate for $|\hat{t}|\gtrsim 4m_{t}^{2}$, if the fixed
$\hat{s}>4m_{t}^{2}$. However, we find that in the vicinity of the point
$|\hat{t}|=4m_{t}^{2}$ the fixed-order results in the $p_{\scriptscriptstyle
T}$ expansion and in the HE expansion are both divergent (see fig. 1). As a
consequence, a straightforward combination of the $p_{\scriptscriptstyle
T}$-expanded and the HE-expanded results cannot allow for an accurate
description of the above region, and this fact prevents a full coverage of the
phase space. We point out that this situation does not change substantially
when higher orders in both the expansions are computed.

Alternatively, the convergence of the expanded results can be improved by
considering the respective Padé approximants. Indeed, starting from a given
Taylor expansion of an exact function $f(x)$ around $x=0$ up to the first $r$
terms

$f(x)\simeq\sum_{k=0}^{r-1}c_{k}x^{k},$ (6)

it is possible to construct the associated Padé approximant, defined as

$[m/n](x)=\frac{p_{0}+p_{1}x+\dots+p_{m}x^{m}}{1+q_{1}x+\dots q_{n}x^{n}},$
(7)

provided that $m+n+1=r$. Specifically, by Taylor-expanding the r.h.s. of eq.
(7), the $\\{p_{i},q_{j}\\}$ coefficients of the Padé approximant can be
written in terms of the $c_{k}$ ones known from eq. (6), by solving a system
of linear equations. Usually, $[m/n]$ Padé approximants such that $m=n$ give
the best improvement in the convergence of the original Taylor expansion, and
we consider only these combinations in our study. In the
$p_{\scriptscriptstyle T}$-expanded results, at NLO, only the first three
terms in eq. (6) are known and therefore we are limited to construct a $[1/1]$
Padé approximant (we will refer to this as the $p_{\scriptscriptstyle
T}$-Padé). Instead the availability of many terms in the HE-expansion results
allows to consider several $[n/n]$ approximants (defined as HE-Padé).

When calculating the $p_{\scriptscriptstyle T}$-Padé, care is to be taken in
the treatment of the expansion parameters. As discussed in refs.[18, 19], not
only the $p_{\scriptscriptstyle T}$ but also the masses of the external
particles are understood as small parameters. Since these are all treated on
the same footing with respect to the large scales set by $\hat{s}$ and
$m_{t}$, we can write the general expression for a $p_{\scriptscriptstyle
T}$-expanded form factor $F$ in the amplitude in terms of a scaling parameter
$x$

$F(x)=\sum_{N=0}^{2}x^{N}\sum_{i+j+k=N}c_{ijk}~{}(p_{\scriptscriptstyle
T}^{2})^{i}(m^{2})^{j}(\Delta_{m})^{k}\equiv\sum_{N=0}^{2}x^{N}c_{N}$ (8)

where $m$ is interpreted as $m_{H}$ and $m_{Z}$ for $gg\rightarrow HH$ and
$gg\rightarrow ZH$, respectively, and $\Delta_{m}=(m_{H}^{2}-m_{Z}^{2})/2$ is
included only for the $ZH$ case (see ref. [19]). Starting from eq. (8) we can
then obtain the corresponding $[1/1]$ Padé approximant with respect to the
limit $x\rightarrow 0$

$[1/1](x)=\frac{p_{0}+p_{1}x}{1+q_{1}x},$ (9)

with

$p_{0}=c_{0}\qquad p_{1}=c_{1}-\frac{c_{0}c_{2}}{c_{1}}\qquad
q_{1}=-\frac{c_{2}}{c_{1}},$

and subsequently set $x=1$ in eq. (9).

We want to clarify a possible source of ambiguity concerning the limit of
validity of the $p_{T}$ expansion. Indeed, while in the previous works we
suggested that this expansion is valid for $p_{T}^{2}\lesssim 4m_{t}^{2}$, as
the comparison at LO between the $p_{\scriptscriptstyle T}$-expanded and exact
result seems to indicate, in this paper we follow a more conservative approach
and we consider as limit of validity for the $p_{\scriptscriptstyle T}$
expansion $|\hat{t}|\lesssim 4m_{t}^{2}$. Additionally, we checked that the
same complementarity for the $p_{\scriptscriptstyle T}$ and HE expansions can
be observed when choosing $p_{T}^{2}=4m_{t}^{2}$ as limit of validity.

We now discuss the procedure adopted to construct the Padé approximants from
the HE expansion. Following the prescription of ref. [28] (see also [29, 30]),
we initially arrange the various orders $F^{(i)}$ of the HE expansion for a
given form factor as follows

$F(x)=F^{(0)}+\sum_{l=1}^{L}\left(F^{(2l-1)}m_{t}^{(2l-1)}+F^{(2l)}m_{t}^{(2l)}\right)x^{l}=\sum_{l=0}^{L}d_{l}x^{l},$
(10)

where orders related to odd powers of $m_{t}$ are grouped with the orders
related to the next even power. Then, we construct $[n/n]$ approximants in $x$
with $2n=L$ from eq. (10) using the analytic expressions available in [31,
32], and setting $x=1$. We remark that our Padé approximants are obtained in a
fully symbolic way, whereas in refs. [28, 29] all the kinematical quantities
are fixed to the respective numerical values before the Padés are constructed
in $x$. Furthermore, in comparison to refs.[28, 22], we only studied $[n/n]$
HE-Padés up to $n=6$. In those references Padés with $n>6$ were also
considered in order to extrapolate the results in the region
$|\hat{t}|<4m_{t}^{2}$, for a fixed $\hat{s}$, and characterize the relative
uncertainties of different $[m/n]$ Padés. In our case, because the region
$|\hat{t}|<4m_{t}^{2}$, is more accurately described by the results of the
$p_{\scriptscriptstyle T}$ expansion, we find that a [6/6] HE-Padé is more
than enough to perform the merging with the $p_{\scriptscriptstyle T}$-result
and, at the same time, to describe accurately the high-energy region.

The $p_{\scriptscriptstyle T}$-Padé and the HE-Padé extend the range of
validity of each expansion beyond its limit. As discussed in the next section,
the $p_{\scriptscriptstyle T}$ and the HE Padés are accurate enough to bridge
the gap around the phase-space region $|\hat{t}|\simeq 4m_{t}^{2}$. Then, an
accurate approximation of the exact result for any phase-space point
($\hat{s},\hat{t}$) can be obtained by choosing as switching point between the
Padé-improved expansions any point in the region $|\hat{t}|\sim 4m_{t}^{2}$.
For simplicity we choose to use the $p_{\scriptscriptstyle T}$-Padé when
$|\hat{t}|<4m_{t}^{2}$ and the HE-Padé when $|\hat{t}|\geq 4m_{t}^{2}$, for
any fixed value of $\hat{s}$. We recall that in our discussion we just
consider the forward region $|\hat{t}|\leq|\hat{u}|$, while the result in the
complementary phase-space region is obtained using the symmetry of our form
factors under $\hat{t}\leftrightarrow\hat{u}$. Noticing that, when
$|\hat{t}|\leq|\hat{u}|$, the maximum absolute value of $\hat{t}$ as a
function of $\hat{s}$ is given by
$|\hat{t}|_{\text{max}}=1/2(\hat{s}-m_{3}^{2}-m_{4}^{2})$ our choice
corresponds to using the $p_{\scriptscriptstyle T}$-Padé up to the partonic
energy $\hat{s}_{c}=8m_{t}^{2}+m_{3}^{2}+m_{4}^{2}$. In this energy region
($\sqrt{\hat{s}_{c}}\simeq 500$ GeV for $gg\rightarrow HH$ and $gg\rightarrow
ZH$) at the LHC more than 2/3 of the hadronic cross section is concentrated.

##  3 The $p_{\scriptscriptstyle T}$ and HE expansions vs the exact results at
LO

In this section we assess the reliability of our merging procedure by studying
how well the combination of the $p_{\scriptscriptstyle T}$-Padé and the HE-
Padé can reproduce the exact LO results for $HH$ and $ZH$ production via gluon
fusion. For the sake of simplicity, we discuss in detail only the
$gg\rightarrow HH$ process, but we verified that similar conclusions can be
drawn for $gg\rightarrow ZH$.

Using the same notation of ref. [16], we recall that the amplitude for
$gg\rightarrow HH$ can be expressed as

$A^{\mu\nu}=\frac{G_{\mu}}{\sqrt{2}}\frac{\alpha_{S}(\mu_{R})}{2\pi}\delta_{ab}T_{F}\hat{s}[A_{1}^{\mu\nu}F_{1}+A_{2}^{\mu\nu}F_{2}],$
(11)

where $F_{1}$ and $F_{2}$ are the form factors associated to the spin-0 and
spin-2 projectors, respectively. Both triangle and box diagrams contribute to
$F_{1}$

$F_{1}=F_{\triangle}\frac{3m_{\scriptscriptstyle
H}^{2}}{\hat{s}-m_{\scriptscriptstyle H}^{2}}+F_{\square},$ (12)

whereas the $F_{2}$ form factor receives contribution only from boxes. Our
goal is to improve the evaluation of the box contributions, therefore we focus
on the discussion of $F_{\square}$ and $F_{2}$.

(a)

(b)

(c)

(d)

Figure 1: Modulus of the box form factors contributing to $gg\rightarrow HH$
at LO, for a fixed value of (a,c) $\sqrt{\hat{s}}=0.9$ TeV and (b,d)
$\sqrt{\hat{s}}=2$ TeV. In the upper part of each plot, the exact prediction
(solid black line) is shown together with the $p_{\scriptscriptstyle T}$ and
HE expansions (solid blue and purple lines, respectively) and with the [1/1]
$p_{\scriptscriptstyle T}$ and [6/6] HE Padé approximants (dashed light blue
and pink lines, respectively). The bottom part of each plot shows the ratio of
the above results to the exact prediction.

The LO results for these form factors, denoted as $F_{\square}^{LO}$ and
$F_{2}^{LO}$, are shown in fig. 1, for fixed values of the partonic center-of-
mass energy. Only large values of $\hat{s}$ are shown in fig. 1 because for
small $\hat{s}$ values the $p_{\scriptscriptstyle T}$-expanded results are
very accurate [18]. The $p_{\scriptscriptstyle T}$-expanded and HE-expanded
results are represented by the blue and purple solid lines, respectively, and
they deviate from the exact result, shown as a solid black line, at
$|\hat{t}|/4m_{t}^{2}\simeq 1$, as anticipated in the previous section. The
light blue dashed line stands for the [1/1] $p_{\scriptscriptstyle T}$-Padé,
while the pink dashed line represents the [6/6] HE-Padé. One can see that the
Padé results show an improved convergence with respect to the fixed-order
expansions. The bottom part of the plots in fig. 1 shows the ratio of the
expanded and Padé results to the exact one. Indeed, fig. 1(a,b) shows that in
the case of $F_{\square}^{LO}$ for $|\hat{t}|/4m_{t}^{2}=1$ the differences of
the Padé results with respect to the exact prediction are negligible. For the
$F_{2}$ form factor, whose contribution to the cross-section is much smaller
than the one of the $F_{1}$ form factor, the difference is always below 5%,
see fig. 1(c,d). We notice that, when comparing the accuracies of the Padé
approximants, larger discrepancies can be attributed to the
$p_{\scriptscriptstyle T}$-Padé. Indeed, being the latter a [1/1] Padé, it is
expected to be a less refined approximation than the [6/6] HE-Padé. Still, in
the case of $F_{\square}^{LO}$ the differences between the two Padé near
$|\hat{t}|/4m_{t}^{2}=1$ are negligible. We also notice that, as $\hat{s}$
increases, larger values of $|\hat{t}|$ are allowed by the kinematics, and the
relative importance of the HE expansion increases.

The improvement in convergence provided by the Padé approximants is such that
the merging of the two results discussed in the previous section can reproduce
the exact prediction with good accuracy for every value of $\hat{t}$, for any
$\hat{s}$. While we refrain from showing more examples here, we note that we
studied the behaviour of all the box contributions to $gg\rightarrow HH$ and
$gg\rightarrow ZH$ at several values of $\hat{s}$. We explicitly checked that,
among the various possibilites, a [6/6] HE-Padé is more than enough for an
accurate merging. Furthermore, we observed that the value
$|\hat{t}|=4m_{t}^{2}$ is a good choice as a merging point for the
$p_{\scriptscriptstyle T}$ and HE Padé approximants. The high level of
accuracy of our merging method can be observed in fig. 2, where the partonic
cross section at LO is shown for $gg\rightarrow HH$ and $gg\rightarrow ZH$.
One can see that deviations of the combination of the $p_{\scriptscriptstyle
T}$\- and HE-Padé with respect to the exact prediction never exceed 1%.

(a)

(b)

Figure 2: Partonic cross section at LO for (a) $gg\rightarrow HH$ and (b)
$gg\rightarrow ZH$. The upper part of each plot shows the exact prediction
(solid line) together with the merging of the $p_{\scriptscriptstyle T}$ and
HE Padé approximants (dashed line). The bottom part of each plot shows the
ratio of the merged result to the exact prediction.

##  4 Merging the $p_{\scriptscriptstyle T}$ and HE expansions at NLO

In the previous section we showed that the merging of the
$p_{\scriptscriptstyle T}$\- and HE-Padé can accurately reproduce the exact LO
prediciton. In this section we present the merging of the NLO
$p_{\scriptscriptstyle T}$-expanded and HE-expanded results improved by the
respective Padé approximants.

(a)

(b)

(c)

(d)

Figure 3: Modulus of the box form factors contributing to $gg\rightarrow HH$
at NLO, for a fixed value of (a,c) $\sqrt{\hat{s}}=0.9$ TeV and (b,d)
$\sqrt{\hat{s}}=2$ TeV. The $p_{\scriptscriptstyle T}$ and HE expansions are
shown as solid blue and purple lines, respectively, while the [1,1]
$p_{\scriptscriptstyle T}$\- and [6,6] HE-Padé are shown as dashed light blue
and pink lines, respectively.

In fig. 3 the NLO contributions to $F_{\square}$ and $F_{2}$ are
shown333$F_{\square}^{NLO}$ and $F_{2}^{NLO}$ are the form factors as defined
in eq. (11) but do not contain the double triangle diagrams that can be
expressed in terms of products of one-loop integrals and as such are computed
analytically in exact top-mass dependence [16]. . The relative behaviour of
the various approximations is analogous to what we observed at LO. For low
values of $|\hat{t}|$ the fixed-order and Padé-improved $p_{\scriptscriptstyle
T}$-expanded results agree well. Increasing the value of $|\hat{t}|$ up to the
merging region, $|\hat{t}|\sim 4m_{t}^{2}$, the $p_{\scriptscriptstyle
T}$-Padé becomes close to the Padé-improved HE expansion. For values above
$|\hat{t}|=4m_{t}^{2}$ the $p_{\scriptscriptstyle T}$ and HE Padé approximants
show small deviations as expected. The NLO study shows the same qualitative
behaviour as the LO one. This makes us confident that the proposed merging
procedure works well also at NLO.

We now compare our evaluation of the virtual corrections for the di-Higgs
production process444We note that for $ZH$ production no public code including
the results of the full computation [11] is currently available. Hence we
refrain from making any comparisons for $ZH$ production. with the numerical
result provided as a grid in ref. [33]. This reference summarizes the work of
ref. [28], where the numerical calculation in exact top-mass dependence of
ref. [8] was supplemented by the result in the HE expansion of ref. [21]. The
comparison is done on the quantity

$\Delta\hat{\sigma}_{virt}=\int_{\hat{t}_{-}}^{\hat{t}_{+}}\frac{\alpha_{s}}{32\pi^{2}}\frac{1}{\hat{s}^{2}}\mathcal{V}_{fin}d\hat{t},$
(13)

where the finite part of the virtual corrections $\mathcal{V}_{fin}$ is
defined as in ref. [26]. The results are shown in fig. 4. The grid of ref.
[33] shows very good agreement with our results at every invariant mass,
except for the first few bins at low $M_{HH}$. The reason is a large
uncertainty of the numerical grid on the low $M_{HH}$ bins, that are described
by only a few points in the numerical grid due to their small contribution to
the total cross section. For moderate and large $M_{HH}$ we observe
differences below 1%.

Figure 4: Upper pannel: $\Delta\hat{\sigma}_{virt}$ using the numerical grid
provided in ref. [28] (dark blue line) and our combination of HE expansion and
small $p_{\scriptscriptstyle T}$ expansion (light blue dashed line). The lower
panel shows the ratio of the two lines.

Finally, we show that our merging approach is flexible with respect to the
modification of the input parameters by computing the virtual corrections for
various renormalisation schemes of the top quark mass. It was noted in refs.
[9, 10] that the di-Higgs production process suffers from a large uncertainty
associated to the renormalisation scheme of the top quark mass. In particular,
an uncertainty on the NLO cross section between $+4\%$ and $-18\%$ [34] is
related to the change from the on-shell renormalisation scheme to the
$\overline{\text{MS}}$ scheme for the top mass, with the latter evaluated at
different values of the renormalisation scale ($m_{t}$, $M_{HH}$ and
$M_{HH}/4$). The results presented so far have been calculated using the on-
shell scheme for the top mass, however the form factors in the
$\overline{\text{MS}}$ scheme can be obtained by simply shifting our result
according to:

$F_{i}^{NLO,\overline{\text{MS}}}=F_{i}^{NLO,\text{OS}}-\frac{1}{4}\frac{\partial
F_{i}^{LO}}{\partial m_{t}^{2}}\Delta_{m_{t}^{2}}$ (14)

with $i=\Delta,\Box,2$ and

$\Delta_{m_{t}^{2}}=2m_{t}^{2}C_{F}\left[-4+3\log\left(\frac{m_{t}^{2}}{\mu^{2}}\right)\right]\,.$
(15)

Notice that the numerical values for the top quark mass have to be adjusted to
their $\overline{\text{MS}}$ values, which we evaluate following refs. [35,
36]. Since the LO results are available analytically, we can first build the
Padé approximants of the $p_{\scriptscriptstyle T}$\- or HE-expanded
$F_{i}^{NLO,\text{OS}}$ form factors and then calculate the shift to the the
$\overline{\text{MS}}$ scheme using the full LO form factors in eq. (14).
Alternatively, we can use the expanded LO form factors to perform the shift in
eq. (14) order by order in the $p_{\scriptscriptstyle T}$-expanded and the HE-
expanded results, and then build Padé approximants on the
$F_{i}^{NLO,\overline{\text{MS}}}$ form factors. The difference between these
two approaches turns out to be well below the 0.5 permille level everywhere
except near the top-mass threshold where the difference is at the percent
level. For the $gg\to ZH$ process, the shift to the $\overline{\text{MS}}$
scheme can be applied in analogy to eq. (14) on the associated form
factors555The $\Delta\hat{\sigma}_{virt}$ for $gg\rightarrow ZH$ was defined
as in ref. [19].. For the triangle contributions we always use the results
available in full top-mass dependence.

We present our results in fig. 5. We observe that indeed the
$\Delta\hat{\sigma}_{virt}$ show a non-negligible dependence on the
renormalisation scheme for the top quark mass. This holds true both for $gg\to
HH$ and for $gg\to ZH$. In particular, for $gg\to ZH$ we see a shift of the
maximal value of the differential $\Delta\hat{\sigma}_{virt}$, which is
associated to the lower top quark mass value. This may point to a significant
uncertainty related to the top-mass scheme also in the case of $gg\to ZH$. We
notice, however, that the effects observed in fig. 5 are likely to be
partially compensated by similar modifications in the LO cross section and by
the inclusion of the real-emission corrections at NLO. Therefore, we use fig.
5 only as an illustration of the flexibility of the merging method discussed
in this paper, and we leave a full assessment of the effects due the change of
the top-mass renormalization scheme to a future work. For the same reason we
refrain from providing a full comparison with ref. [34] for $gg\to HH$.

(a)

(b)

Figure 5: The integrated virtual corrections in different renormalisation
schemes for the top quark mass for (a) $gg\rightarrow HH$ and (b)
$gg\rightarrow ZH$. The dashed lines show $\Delta\hat{\sigma}_{virt}$ in the
$\overline{\text{MS}}$ scheme using the top quak mass evaluated at different
choices of the renormalisation scale, while the black solid line shows the on-
shell result.

##  5 Conclusion

In this paper, we combined an expansion in small $p_{\scriptscriptstyle T}$
with a HE expansion for the processes $gg\to HH$ and $gg\to ZH$ and we showed
that this combination leads to results that describe the whole phase space of
the considered processes with high accuracy.

The expansion in small $p_{\scriptscriptstyle T}$ is valid for
$|\hat{t}|\lesssim 4m_{t}^{2}$, while the HE expansion is valid for
$|\hat{t}|\gtrsim 4m_{t}^{2}$ and $\hat{s}>4m_{t}^{2}$. For a successful
combination of the two expansion methods, we extended the validity range of
both expansions making use of simple Padé approximants. This allowed us to
describe also the region $|\hat{t}|\simeq 4m_{t}^{2}$ very accurately.

We verified our method first at LO, where exact analytic results are
available. At the level of the differential partonic cross section, the
difference with the results in exact top quark mass dependence never exceeded
1% for any of the $M_{HH}$ bins. In a second step, we compared our results
with a numerical grid containing the results of a computation in exact top
quark mass dependence for the $gg\to HH$ process [28]. Also at NLO we could
verify the good numerical accuracy of our approach, observing differences
below 1% for the invariant mass bins where the numerical accuracy of the grid
is assumed to be sufficiently high.

We have also shown that in our analytic calculation the top quark mass
renormalisation scheme can be easily changed. This shows the great flexibility
of our analytic approach with respect to the numerical one. With a running
time of well below 1 s per phase-space point our results can be well
implemented in a fast and versatile Monte Carlo program.

## Acknowledgements

We thank Lina Alasfar and Xiaoran Zhao for useful discussions. We would also
like to thank the authors of ref. [22] for a careful reading of the manuscript
and for their useful comments. G.D. would like to thank the Department of
Physics of the University of Rome ”La Sapienza” for the kind hospitality
during part of this project. The work of G.D. and M.V. was partially supported
by the Italian Ministry of Research (MUR) under grant PRIN 20172LNEEZ. The
work of L.B. and P.P.G. has received financial support from Xunta de Galicia
(Centro singular de investigación de Galicia accreditation 2019-2022), by
European Union ERDF, and by “María de Maeztu” Units of Excellence program
MDM-2016-0692 and the Spanish Research State Agency.

## References

  * [1] S. Catani, D. de Florian, M. Grazzini, and P. Nason, “Soft gluon resummation for Higgs boson production at hadron colliders,” JHEP 07 (2003) 028, arXiv:hep-ph/0306211.
  * [2] V. Ahrens, T. Becher, M. Neubert, and L. L. Yang, “Origin of the Large Perturbative Corrections to Higgs Production at Hadron Colliders,” Phys. Rev. D 79 (2009) 033013, arXiv:0808.3008 [hep-ph].
  * [3] M. Spira, A. Djouadi, D. Graudenz, and P. M. Zerwas, “Higgs boson production at the LHC,” Nucl. Phys. B453 (1995) 17–82, arXiv:hep-ph/9504378 [hep-ph].
  * [4] R. Harlander and P. Kant, “Higgs production and decay: Analytic results at next-to-leading order QCD,” JHEP 0512 (2005) 015, arXiv:hep-ph/0509189 [hep-ph].
  * [5] U. Aglietti, R. Bonciani, G. Degrassi, and A. Vicini, “Analytic Results for Virtual QCD Corrections to Higgs Production and Decay,” JHEP 01 (2007) 021, arXiv:hep-ph/0611266 [hep-ph].
  * [6] C. Anastasiou, S. Beerli, S. Bucherer, A. Daleo, and Z. Kunszt, “Two-loop amplitudes and master integrals for the production of a Higgs boson via a massive quark and a scalar-quark loop,” JHEP 01 (2007) 082, arXiv:hep-ph/0611236 [hep-ph].
  * [7] S. Borowka, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk, and T. Zirke, “SecDec-3.0: numerical evaluation of multi-scale integrals beyond one loop,” Comput. Phys. Commun. 196 (2015) 470–491, arXiv:1502.06595 [hep-ph].
  * [8] S. Borowka, N. Greiner, G. Heinrich, S. Jones, M. Kerner, J. Schlenk, U. Schubert, and T. Zirke, “Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence,” Phys. Rev. Lett. 117 no. 1, (2016) 012001, arXiv:1604.06447 [hep-ph]. [Erratum: Phys. Rev. Lett.117,no.7,079901(2016)].
  * [9] J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, M. Spira, and J. Streicher, “Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme,” Eur. Phys. J. C 79 no. 6, (2019) 459, arXiv:1811.05692 [hep-ph].
  * [10] J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca, M. Spira, and J. Streicher, “Higgs-Pair Production via Gluon Fusion at Hadron Colliders: NLO QCD Corrections,” JHEP 04 (2020) 181, arXiv:2003.03227 [hep-ph].
  * [11] L. Chen, G. Heinrich, S. P. Jones, M. Kerner, J. Klappert, and J. Schlenk, “$ZH$ production in gluon fusion: two-loop amplitudes with full top quark mass dependence,” JHEP 03 (2021) 125, arXiv:2011.12325 [hep-ph].
  * [12] S. Dawson, S. Dittmaier, and M. Spira, “Neutral Higgs boson pair production at hadron colliders: QCD corrections,” Phys. Rev. D58 (1998) 115012, arXiv:hep-ph/9805244 [hep-ph].
  * [13] L. Altenkamp, S. Dittmaier, R. V. Harlander, H. Rzehak, and T. J. E. Zirke, “Gluon-induced Higgs-strahlung at next-to-leading order QCD,” JHEP 02 (2013) 078, arXiv:1211.5015 [hep-ph].
  * [14] J. Grigo, J. Hoff, K. Melnikov, and M. Steinhauser, “On the Higgs boson pair production at the LHC,” Nucl. Phys. B875 (2013) 1–17, arXiv:1305.7340 [hep-ph].
  * [15] J. Grigo, J. Hoff, and M. Steinhauser, “Higgs boson pair production: top quark mass effects at NLO and NNLO,” Nucl. Phys. B900 (2015) 412–430, arXiv:1508.00909 [hep-ph].
  * [16] G. Degrassi, P. P. Giardino, and R. Gröber, “On the two-loop virtual QCD corrections to Higgs boson pair production in the Standard Model,” Eur. Phys. J. C76 no. 7, (2016) 411, arXiv:1603.00385 [hep-ph].
  * [17] A. Hasselhuhn, T. Luthe, and M. Steinhauser, “On top quark mass effects to $gg\to ZH$ at NLO,” JHEP 01 (2017) 073, arXiv:1611.05881 [hep-ph].
  * [18] R. Bonciani, G. Degrassi, P. P. Giardino, and R. Gröber, “Analytical Method for Next-to-Leading-Order QCD Corrections to Double-Higgs Production,” Phys. Rev. Lett. 121 no. 16, (2018) 162003, arXiv:1806.11564 [hep-ph].
  * [19] L. Alasfar, G. Degrassi, P. P. Giardino, R. Gröber, and M. Vitti, “Virtual corrections to $gg\to ZH$ via a transverse momentum expansion,” JHEP 05 (2021) 168, arXiv:2103.06225 [hep-ph].
  * [20] J. Davies, G. Mishima, M. Steinhauser, and D. Wellmann, “Double-Higgs boson production in the high-energy limit: planar master integrals,” JHEP 03 (2018) 048, arXiv:1801.09696 [hep-ph].
  * [21] J. Davies, G. Mishima, M. Steinhauser, and D. Wellmann, “Double Higgs boson production at NLO in the high-energy limit: complete analytic results,” JHEP 01 (2019) 176, arXiv:1811.05489 [hep-ph].
  * [22] J. Davies, G. Mishima, and M. Steinhauser, “Virtual corrections to $gg\to ZH$ in the high-energy and large-$m_{t}$ limits,” JHEP 03 (2021) 034, arXiv:2011.12314 [hep-ph].
  * [23] X. Xu and L. L. Yang, “Towards a new approximation for pair-production and associated-production of the Higgs boson,” JHEP 01 (2019) 211, arXiv:1810.12002 [hep-ph].
  * [24] G. Wang, Y. Wang, X. Xu, Y. Xu, and L. L. Yang, “Efficient computation of two-loop amplitudes for Higgs boson pair production,” Phys. Rev. D 104 no. 5, (2021) L051901, arXiv:2010.15649 [hep-ph].
  * [25] G. Wang, X. Xu, Y. Xu, and L. L. Yang, “Next-to-leading order corrections for $gg\to ZH$ with top quark mass dependence,” arXiv:2107.08206 [hep-ph].
  * [26] R. Gröber, A. Maier, and T. Rauh, “Reconstruction of top-quark mass effects in Higgs pair production and other gluon-fusion processes,” JHEP 03 (2018) 020, arXiv:1709.07799 [hep-ph].
  * [27] J. Davies, R. Gröber, A. Maier, T. Rauh, and M. Steinhauser, “Padé approach to top-quark mass effects in gluon fusion amplitudes,” PoS RADCOR2019 (2019) 079, arXiv:1912.04097 [hep-ph].
  * [28] J. Davies, G. Heinrich, S. P. Jones, M. Kerner, G. Mishima, M. Steinhauser, and D. Wellmann, “Double Higgs boson production at NLO: combining the exact numerical result and high-energy expansion,” JHEP 11 (2019) 024, arXiv:1907.06408 [hep-ph].
  * [29] J. Davies, G. Mishima, M. Steinhauser, and D. Wellmann, “$gg\to ZZ$: analytic two-loop results for the low- and high-energy regions,” JHEP 04 (2020) 024, arXiv:2002.05558 [hep-ph].
  * [30] D. Wellmann, Top Quark Mass Effects in Higgs and Z Boson Pair Production and Higgs Boson Decays.  PhD thesis, KIT, Karlsruhe, 2020.
  * [31] https://www.ttp.kit.edu/preprints/2019/ttp19-018/.
  * [32] https://www.ttp.kit.edu/preprints/2020/ttp20-041/.
  * [33] https://github.com/mppmu/hhgrid.
  * [34] J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca, and M. Spira, “$gg\to HH$ : Combined uncertainties,” Phys. Rev. D 103 no. 5, (2021) 056002, arXiv:2008.11626 [hep-ph].
  * [35] K. Melnikov and T. v. Ritbergen, “The Three loop relation between the MS-bar and the pole quark masses,” Phys. Lett. B 482 (2000) 99–108, arXiv:hep-ph/9912391.
  * [36] M. Carena, D. Garcia, U. Nierste, and C. E. M. Wagner, “Effective Lagrangian for the $\bar{t}bH^{+}$ interaction in the MSSM and charged Higgs phenomenology,” Nucl. Phys. B 577 (2000) 88–120, arXiv:hep-ph/9912516.