N 3 LO soft-gluon corrections in single-particle-inclusive kinematics and H + H − production

We calculate the complete soft-gluon corrections for the production of colorless final states through N 3 LO in single-particle-inclusive kinematics. We present explicit analytical results and use them to study higher-order QCD corrections for the production of a heavy charged Higgs pair ( H + H − ) via quark-antiquark annihilation in the Two-Higgs-Doublet Model at LHC energies. We calculate the NNLO soft-gluon and virtual QCD corrections as well as the N 3 LO soft-gluon corrections to the total cross section and the charged-Higgs rapidity distribution. This is the first calculation of complete N 3 LO soft-gluon corrections for a process in single-particle-inclusive kinematics, and the results can be applied to other processes with colorless final states.


Introduction
Precision in calculations of cross sections for hard-scattering processes requires the inclusion of higher-order QCD corrections.A particular class of these corrections arises from the emission of soft gluons.These soft radiative corrections can be formally resummed [1][2][3][4][5] to all orders in the strong coupling after taking Mellin moments or Laplace transforms of the cross section and using its factorization properties, and the resummed cross section can be used as a generator of fixed-order corrections.Soft-gluon resummation has been successfully developed and applied to processes with massive final states, such as top-quark pair production [3][4][5][6][7][8], because the cross section receives large contributions from soft-gluon emission near partonic threshold due to the large mass.In all these processes, the soft-gluon corrections are numerically dominant and account for the overwhelming majority of the complete corrections at higher orders.
For processes with complex color flows, such as top production, one does not yet have all the ingredients necessary to calculate the complete set of soft-gluon corrections at next-to-next-to-nextto-leading order (N 3 LO).For processes with trivial color flow, and in particular with colorless (e.g.electroweak) final states, the N 3 LO soft-gluon corrections are known in fully inclusive kinematics [9][10][11].However, complete N 3 LO soft-gluon contributions have not yet been calculated for a partonic process in single-particle-inclusive (1PI) kinematics.In this paper we provide explicit analytical results for the complete soft-gluon corrections through N 3 LO for 1PI partonic processes with no color in the final state.As a concrete application, we employ these results to compute higher-order QCD corrections for charged-Higgs pair production in the Two-Higgs-Doublet Model (2HDM) using 1PI kinematics.
In this work we present the most precise predictions for quark-induced production of H + H − at the LHC, that occurs through diagrams of Drell-Yan type, in the 2HDM by computing the total cross section and the charged-Higgs rapidity distribution for this process at approximate next-tonext-to-leading order (aNNLO) and approximate N 3 LO (aN 3 LO) in QCD.Results at aNNLO are derived by adding the second-order soft-plus-virtual corrections to the full next-to-leading order (NLO) QCD calculation, while aN 3 LO results are obtained by further adding the third-order softgluon corrections.
The 2HDM represents one of the simplest extensions of the Standard Model Higgs sector, introducing an additional SU (2) L Higgs doublet.This model gives rise to a very rich phenomenology due to its extended Higgs sector (for a review see [12] and references therein).The model predicts the existence of three neutral scalar bosons (h, H, and A) and a pair of charged Higgs bosons (H + and H − ).
The dominant production channel for the charged Higgs at the Large Hadron Collider (LHC) is in association with a top quark p p → t H − , for which soft-gluon corrections were calculated in Refs.[13][14][15], followed by the associated production with a W boson p p → H + W − (see [16] for soft-gluon corrections), and pair production p p → H + H − .In this work we focus on the last process, which can occur at leading order (LO) through quark-antiquark annihilation (q q → H + H − ) [17][18][19] and gluon fusion (gg → H + H − ) [20][21][22][23][24]. In the case of light quarks in the initial state, the pair production of charged Higgs bosons is of Drell-Yan type, with the exchange of a photon and a Z boson in the s-channel.NLO QCD corrections for this channel have been computed in the 2HDM and the Minimal Supersymmetric Standard Model (MSSM) in [17,19].In the case of bottom quarks in the initial states, there are additional diagrams that account for the exchange of a neutral Higgs boson in the s-channel and a top quark in the t-channel.The NLO QCD corrections to b b → H + H − have been computed in [18,19].The charged-Higgs pair production through gluon fusion is a one-loop process already at LO and it has been computed in the 2HDM and the MSSM in [20][21][22][23][24].
In models like the type II 2HDM, the quark-antiquark annihilation process q q → H + H − is the dominant channel for low values of tan β, which is the ratio of the vacuum expectation values for the two doublets.On the other hand, for large values of tan β, the gg and b b channels have cross sections which are enhanced by a tan 4 β factor, and this can compensate for the loop suppression or the small bottom-quark parton densities and can make these channels dominant, as shown in [19].However, recent studies showed that, for type II 2DHM with softly broken Z 2 symmetry, the size of cos(β − α) has to be lower than O(10 −2 ) and that low values of tan β are favored [25][26][27][28].Therefore, the calculations presented in this work are particularly relevant for those scenarios.
The paper is organized as follows.In Section 2, we describe the soft-gluon resummation formalism and provide explicit analytical results for the soft-gluon corrections through N 3 LO in 1PI kinematics for processes with colorless final states.In Section 3 we introduce the Two-Higgs-Doublet model and apply the resummation formalism of Section 2 to charged-Higgs pair production.In Section 4, we present numerical results for the total cross section through aN 3 LO at LHC energies.In Section 5, we present results for the charged-Higgs rapidity distribution for two representative masses of the charged Higgs boson.We conclude in Section 6.

Soft-gluon corrections at N 3 LO in 1PI kinematics
In this section, we discuss the resummation formalism that we use for the calculation of softgluon corrections in 1PI kinematics for partonic processes with colorless final states, e.g. with electroweak final-state particles, such as H + H − production.The origin of the soft-gluon corrections is from the emission of soft (i.e.low-energy) gluons from the initial state partons, which result in partial cancellations of infrared divergences between real-emission and virtual diagrams close to partonic threshold.We derive the resummed cross section via Laplace transforms, factorization, and renormalization-group evolution.We then provide explicit analytical results for the soft-gluon corrections through N 3 LO.
We consider first leading-order partonic processes of the form where q and q are quarks and antiquarks in the protons, and A and B are a pair of colorless particles, each of mass m (e.g.H + H − ), with particle A the observed particle in 1PI kinematics.We will also discuss later the changes needed in the analytical expressions if we instead have gluons in the initial state, i.e. processes gg → AB, or if the two final-state particles have different masses, m A and m B .
We The resummation of soft-gluon corrections is a consequence of the factorization properties of the (in general, differential) cross section under Laplace transforms in 1PI kinematics.We first write the differential hadronic cross section, dσ pp→AB , as a convolution of the differential partonic cross section, dσ q q→AB , with the parton distribution functions (pdf), ϕ q/p and ϕ q/p , as where µ F is the factorization scale, and x a , x b are momentum fractions of partons q, q, respectively, in the colliding protons.We also define the hadron-level variables S = (P , where P a and P b denote the momenta of the colliding protons, and We then consider the parton-parton cross section dσ q q→AB , which is of the same form as Eq.(2.2) but with incoming partons instead of hadrons [3][4][5], dσ q q→AB (S 4 ) = dx a dx b ϕ q/q (x a ) ϕ q/q (x b ) dσ q q→AB (s 4 ) , and we define its Laplace transform as Using the expression for S 4 /S in Eq. ( 2.3), we can rewrite the Laplace transform, Eq. (2.5), of the parton-parton cross section in Eq. (2.4) as where N a = N (−u 1 /s) and N b = N (−t 1 /s).
Next, we introduce a refactorization of the cross section via new functions H q q→AB , S q q→AB , ψ q/q , ψ q/q [3][4][5].The hard function H q q→AB is purely short-distance and infrared safe while the soft function S q q→AB describes the emission of noncollinear soft gluons.The coupling of the soft gluons to the partons in the hard-scattering process is described by eikonal (Wilson) lines, i.e. ordered exponentials of the gauge field.The functions ψ q/q and ψ q/q differ from the pdf ϕ q/q and ϕ q/q , and they describe collinear emission from the incoming partons [1,[3][4][5].The refactorized form of the cross section [5,29] is then dσ q q→AB = dw a dw b dw S ψ q/q (w a ) ψ q/q (w b ) where the w's are dimensionless weights, with w a and w b for ψ q/q and ψ q/q , respectively, and w S for S q q→AB .The argument in the delta function of Eq. (2.7) arises from rewriting Eq. (2.3) in terms of the new weights [5], as After taking a Laplace transform of Eq. (2.7), we have dw a e −Nawa ψ q/q (w a ) All N -dependence has now been absorbed into the functions ψq/q , ψq/q , and S. Comparing Eqs.(2.6) and (2.9), we get an expression for the hard-scattering partonic cross section in Laplace transform space, d σqq→AB (N ) = ψq/q (N a ) ψq/q (N b ) φq/q (N a ) φq/q (N b ) (2.10) The resummed differential cross section in Laplace-transform space is derived from the renormalizationgroup evolution of ψq/q / φq/q , ψq/q / φq/q , and Sqq→AB in Eq. (2.10), and it is given by the expression [1][2][3][4][5] In the first exponential of Eq. (2.11), we have with an analogous expression for E q(N b ) [1].The first integrand in Eq. (2.12) has the perturbative expansion A q = (α s /π)A (1)  q + (α The first term in this expansion is A (1)  q = C F [1], where C F = (N 2 c − 1)/(2N c ) and N c = 3 is the number of colors.The second term, A (2)  q , is given by [2, 30] where C A = N c and n f is the number of light-quark flavors (in the next section we will set n f = 5 for H + H − production).Finally, the third term, A (3)  q , is given by [31] Here and below The corresponding expressions for processes with gluons in the initial state (i.e.gg → AB) are given at these perturbative orders by q , where n = 1, 2, 3.The second integrand in Eq. (2.12) can be expanded as D q = (α s /π)D (1)  q + (α s /π) 2 D (2)  q + (α s /π) 3 D (3)  q + • • •.In Feynman gauge, the one-loop term D (1)  q = 0, while at two loops [32] and at three loops [9] The corresponding expressions for processes with gluons in the initial state are given at these perturbative orders by q , where n = 1, 2, 3.In the second exponent of Eq. (2.11), the quantity γ q/q is the moment-space anomalous dimension of the MS density ϕ q/q [33-37] which can be expressed as γ q/q (N ) = −A q ln Ñ + γ q , and similarly for γ q/q .The parton anomalous dimension can be expanded as γ q = (α s /π)γ (1)  q + (α s /π) 2 γ (2)  q + • • • with γ (1)  q = 3C F /4 and (2.17) For gluons, we have γ (1)  g = β 0 /4, where β 0 = 11C A /3 − 2n f /3 is the coefficient of the first term in the perturbative expansion of the β function [38,39], and By expanding Eq. (2.11) to fixed order and inverting back to momentum space, we can calculate the soft-gluon corrections completely through N 3 LO.The Born differential cross section at partonic level can be written as where the specific form of F B q q→AB depends on the process (see the next section for H + H − production).If we perturbatively expand H q q→AB = H (0) The NLO soft-plus-virtual corrections are given by where µ R is the renormalization scale and The coefficients of the D 1 (s 4 ) and D 0 (s 4 ) terms in Eq. (2.20) are, respectively, The coefficient of the δ(s 4 ) term is where q q→AB S (0) q q→AB /F B q q→AB is the contribution from the one-loop virtual corrections.For Drell-Yan type processes with quark-antiquark annihilation, such as [40][41][42].For gluon-initiated processes, we replace C F by C A in Eqs.(2.22) and (2.23), and we use the corresponding V 1 (for example, see [43,44] for the NLO virtual corrections for gg → H).
The NNLO soft-plus-virtual corrections are given by where q q→AB S (0) q q→AB + H (1) q q→AB /F B q q→AB is the contribution from the virtual two-loop corrections.For Drell-Yan type processes with quark-antiquark annihilation, such as H + H − production, V 2 is given by the expression [45] For processes gg → AB, we replace in Eq. ( 2.24) all explicit appearances of C F by C A , and we also use the corresponding expressions for A g , D g , γ g , c 3 , c 2 , and c 1 as detailed above, as well as the appropriate expressions for V 1 and V 2 (for example, see [46] for the NNLO virtual corrections for gg → H).
The N 3 LO soft-gluon corrections are given by where and with [47][48][49] the coefficient of the second term in the perturbative expansion of the β function.For processes gg → AB, we replace in Eqs.(2.26), (2.27), (2.28), and (2.29) all explicit appearances of C F by C A , and we use the corresponding expressions for A g , D g , γ g , c 3 , c 2 , and c 1 as detailed above, as well as the appropriate expressions for V 1 and V 2 .
Finally, we consider the more general case where the observed particle A has mass m A , which is different from the mass m B of the other final-state particle B. In that case, we simply replace m by m A everywhere it appears explicitly in Eqs.(2.21) through (2.29), and we also redefine and use B and u 1 = u − m 2 B everywhere.With these simple changes all the expressions above can be used for the more general case of unequal masses.
It is worth emphasizing that the expressions for the higher-order soft-gluon corrections presented in this section are universal and can be applied to any electroweak production process that occurs through quark-antiquark annihilation or gluon-gluon fusion.The only process-dependent parts are encoded into the virtual correction terms V 1 and V 2 , which depend on the specific process.For example, explicit expressions of those terms for various electroweak processes can be found in [11].

Charged Higgs pair production in the 2HDM
In this section we discuss the production of a charged Higgs pair via quark-antiquark annihilation in the 2HDM that occurs through diagrams of Drell-Yan type.For concreteness, we consider a CP-conserving type II 2HDM with a softly-broken Z 2 symmetry defined by the transformations Φ 1 → Φ 1 and Φ 2 → −Φ 2 , where Φ 1 and Φ 2 are the two Higgs doublets.In this model, the Higgs doublet Φ 1 is responsible for the Yukawa interactions that generate masses for the down-type quarks and charged leptons, while Φ 2 is responsible for the Yukawa interactions that generate masses for the up-type quarks.The Higgs potential is given by: where all the masses (m 1 , m 2 , m 12 ) and quartic couplings (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ) are taken to be real in order to avoid CP violation.After electroweak symmetry breaking, this model predicts the existence of three neutral scalar fields (two CP-even states h and H, and a CP-odd state A) and a pair of charged Higgs bosons H ± .The eight parameters of the scalar Higgs potential in Eq. (3.1) are usually rewritten in terms of more physical ones: a common choice is to consider the electroweak vacuum expectation value υ = 246 GeV, the four Higgs masses m h , m H , m A , m H + , the tangent tan β, the cosine cos(β − α), and the discrete symmetry soft-breaking term m 12 (for more details see Section III of [50]).For this particular model, recent studies [25][26][27][28] showed that the size of cos(β − α) has to be lower than O(10 −2 ) and that low values of tan β are favored.As already mentioned, charged-Higgs pair production is a process that can occur through quarkantiquark annihilation and gluon-fusion channels.In this work we consider only quark-induced production which is the dominant production mechanism at low tan β values [19].Leading-order Feynman diagrams for quark-induced production are shown in Fig. 1.In the case of light quarks in the initial state, we have diagrams of Drell-Yan type where the production of a pair of charged Higgs bosons occurs through the exchange of a photon and a Z boson in the s-channel (left diagram of Fig. 1).In the case of bottom quarks, in addition to Drell-Yan type diagrams, the production of charged Higgs bosons also occurs through the exchange of the CP-even neutral Higgs bosons (h, H) in the s-channel (middle diagram of Fig. 1) and the top-quark in the t-channel (right diagram of Fig. 1).These extra diagrams might interfere with Drell-Yan type ones and introduce some dependence on the scalar potential parameters tan β, cos(β − α) and m 2 12 .In this work, we consider only Drell-Yan type production of a pair of charged Higgs bosons in the 5-flavor scheme where we include also a massless bottom quark in the initial state.In this way, our H + H − production process depends only on the mass of the final-state charged Higgs and is independent of the other 2HDM scalar potential parameters, like tan β, cos(β − α), and m 2 12 .In considering the total contribution of both light and heavy quarks in the initial state, we are neglecting some bottom-quark contributions (middle and right diagrams of Fig. 1), which might give a sizable contribution, depending on the specific values of the Higgs potential parameters.
However, for small values of tan β favored by recent fits, i.e. 2 < tan β < 4, these contributions are of the order of 2-3% at LO and NLO, as also shown in [19], which are smaller than the aNNLO corrections for the Drell-Yan type processes and within the theoretical uncertainty of our calculation obtained from scale and pdf variation.We would like to emphasize that the goal of this section is to provide the most precise prediction for processes that occur through Drell-Yan type diagrams, which for light quarks represent their complete contribution.Providing the total b b contribution, which depends on the values of several Higgs potential parameters, is beyond the scope of this paper.
Therefore, the quantity F B q q→AB entering in the Born differential cross section in Eq. (2.19) computed using Drell-Yan type diagrams is given by for up-type quarks, while for down-type quarks we have These expressions are used in Eq. (2.19) to compute the Born cross section, and in Eqs.(2.20), (2.24), and (2.26) to derive the higher-order soft-gluon corrections through N 3 LO.
4 Total cross sections for H + H − production at the LHC In this section we present results for the total cross section of charged-Higgs pair production via quark-antiquark annihilation q q → H + H − at the LHC, for two center-of-mass energies, namely 13 TeV and 13.6 TeV.As already discussed in the previous section, this cross section depends only on m H + and is independent of the other Higgs potential parameters.In our calculation we work in the 5-flavor scheme and we include the bottom quark in the initial state, which is taken to be massless.Lower bounds on the mass of the charged Higgs of the order of 570-800 GeV have been obtained in Ref. [51].In light of these constraints, the charged-Higgs mass is taken here to vary between 500 and 1500 GeV.Lower values of m H + are excluded, while higher values, although allowed, are not considered because they give very small cross sections.The complete LO and NLO QCD results are calculated using MadGraph5 aMC@NLO [52] and an ad hoc UFO model for the CP-conserving type II 2HDM, which can be found at [53].We take m W = 80.377 GeV, m Z = 91.1876GeV, and α −1 = 127.9.We use the same MSHT20 aN 3 LO pdf set [54] for computing results at every perturbative order, in order to show how each order in the series contributes to the aN 3 LO cross section.
For the total cross section, the central results are obtained by setting the factorization and renormalization scales both equal to m H + .Scale uncertainties are obtained by varying independently µ R and µ F between m H + /2 and 2m H + , and pdf uncertainties are also computed.The cross sections at aNNLO are computed by adding the second-order soft-plus-virtual QCD corrections to the complete NLO result.The third-order soft-gluon corrections are further added to derive the result at aN 3   In Fig. 2 and Fig. 3 we show the total cross section for H + H − production at the LHC with 13 TeV and 13.6 TeV energy, respectively.These plots show the central value, obtained with µ = m H + (where µ = µ R = µ F ), of the LO, NLO, aNNLO, and aN 3 LO cross sections as a function of the mass of the charged Higgs boson.The inset plots display the K-factors of the higher-order cross sections relative to the LO result.At each perturbative order, we find that the K-factors slowly increase with increasing charged-Higgs mass.For instance, at both 13 TeV and 13.6 TeV, the NLO K-factor varies from 1.15 for m H + = 500 GeV to 1.20 for m H + = 1500 GeV.The NNLO soft-plus-virtual corrections provide an additional 5% to 6% increase, and the N 3 LO soft-gluon corrections give a further 0.5% to 0.6% increase.Thus, the aN 3 LO K-factor varies from 1.21 for m H + = 500 GeV to 1.26 for m H + = 1500 GeV.
In Table 1 and Table 2 we present some numerical values for the total cross section of H + H − production in proton-proton collisions, with a center-of-mass energy of 13 TeV and 13.6 TeV, respectively.These cross sections are shown for select charged-Higgs masses, namely 600, 800, 1000, and 1200 GeV.They are calculated at LO, NLO, aNNLO, and aN 3 LO and they are shown together with scale and pdf uncertainties.We want to stress that the soft-gluon corrections are numerically dominant.The result derived with only NLO soft-gluon corrections differs from the exact NLO cross section by only around one percent or less for the masses and energies considered here.Furthermore, the aNNLO result would change by only 2 per mille if we did not include the δ(s 4 ) terms (i.e. the virtual corrections) at NNLO; the soft-gluon corrections are numerically by far the dominant contribution to the higher-order QCD corrections.As shown in Table 1 for the results at 13 TeV energy, the scale uncertainty at LO is large and is reduced significantly at NLO.The higher-order soft-plus-virtual corrections at NNLO decrease the uncertainty even more.Finally, the addition of the N 3 LO soft-gluon corrections further decreases this uncertainty for larger charged-Higgs masses but not for smaller ones.The pdf uncertainty varies from about ±3% for a mass of 600 GeV to ±7% for a 1200 GeV mass, at each perturbative order.The results in Table 2 at 13.6 TeV energy show similar behavior for the scale uncertainty at different perturbative orders.Also in this case, the pdf uncertainty varies at each perturbative order from about ±3% for m H + = 600 GeV to ±6% for m H + = 1200 GeV.
Regarding pdf uncertainties, we can see that, for both center-of-mass energies, they increase with increasing mass of the charged Higgs boson.Moreover, they become bigger than the scale ones already at NLO and therefore, starting at that order, they should be considered the dominant pp → H + H − cross section via quark-antiquark annihilation at 13  source of theoretical uncertainty.As argued in [54], with the use of aN 3 LO pdf, the factorisation scale variation is already contained within the predicted pdf uncertainties.
It is also interesting to study the relative contributions of the different logarithmic powers to the higher-order corrections.As an example, we consider a charged-Higgs mass of 600 GeV and an energy of 13.6 TeV.Setting µ F = µ R = m H + and using the notation of Eq. (2.21), the NNLO corrections from the D 3 terms alone are 0.0134 fb, from the D 3 terms plus the D 2 terms they are 0.0189 fb, from D 3 + D 2 + D 1 they are 0.0150 fb, and from D 3 + D 2 + D 1 + D 0 they are 0.0118 fb.By further adding the NNLO δ(s 4 ) terms, we get the final NNLO soft-plus-virtual correction of 0.0123 fb.Thus, we see that the virtual contribution to the NNLO soft-plus-virtual corrections is small, and the soft-gluon corrections alone are dominant.We also see that the leading-logarithmic contribution (i.e. from D 3 ) alone is quite close to the complete soft-plus-virtual contribution at NNLO.
For the N 3 LO soft-gluon contributions, again with a 600 GeV charged-Higgs mass at 13.6 TeV energy, we find a contribution of 0.00993 fb from the Thus, at N 3 LO, the leading-logarithmic contribution (i.e. from D 5 ) is very different from the complete soft-gluon correction, in contrast to what we found above for the NNLO corrections.We also see that the N 3 LO soft-gluon corrections are an order of magnitude smaller than those at NNLO.
As mentioned in the previous section, for small values of tan β, the additional contributions from b b → H + H − are rather small and less than the aNNLO corrections of the Drell-Yan type q q → H + H − diagrams.One can calculate these b b contributions at NLO for any choice of the Higgs-potential parameters, and then add them to our aN 3 LO result for the Drell-Yan type q q contributions.For example, for tan β = 3, m H + = 600 GeV, m H = 400 GeV, cos(β − α) = 0, and m 12 = 0 GeV, the missing b b contribution at 13 TeV energy is 0.0035 fb.Looking at the first line of Table 1, we see that this is much smaller than the aNNLO corrections of 0.011 fb for q q as well as the total uncertainty of the overall q q cross section.Adding the b b contribution to the aN 3 LO q q cross section gives a total of 0.232 fb.Clearly, in this case it is not necessary to calculate the aNNLO and aN 3 LO corrections for the b b contribution since they would be much smaller than the precision of the numbers displayed.

Charged-Higgs rapidity distributions
As the formalism presented in Section 2 indicates, soft-gluon resummation is derived not only for total cross sections but also for differential distributions.In this section we present results for the charged-Higgs rapidity distribution in H + H − production at the LHC with 13.6 TeV center-of-mass energy, for two specific charged-Higgs masses, namely m H + = 600 GeV and m H + = 800 GeV.Again, we use MadGraph5 aMC@NLO [52] for the complete NLO results to which we add the second-order soft-plus-virtual corrections to derive aNNLO distributions.The third-order soft-gluon corrections are also added in order to derive results at aN 3 LO.In Fig. 4, we plot the central (µ = m H + ) results at LO, NLO, aNNLO, and aN 3 LO for the H + rapidity distribution in H + H − production at 13.6 TeV LHC energy with MSHT20 aN 3 LO pdf, for m H + equal to 600 GeV.The K-factors relative to the LO results are shown in the inset plot.The enhancement from the aNNLO soft-plus-virtual corrections is significant while the aN 3 LO soft-gluon contribution is smaller than 1% in the rapidity range shown in the figure, similar to the total cross section case.We see that the K-factors are relatively flat for central values of the charged-Higgs rapidity, but they begin to grow for rapidities larger than 1.We also note that the distribution becomes very small for rapidity values above 2.
Regarding the theoretical uncertainties from scale variation, we note that they are essentially the same as those for the total cross section in the rapidity range plotted, and only become bigger at very large rapidities, where the numerical value of the distribution is very small.As for the total cross section, the higher-order corrections reduce the scale dependence.The pdf uncertainties are also similar to those of the total cross section, around 3%, in the rapidity range shown in the plot,  but they increase at larger rapidities.In Fig. 5, we plot the central (µ = m H + ) results at LO, NLO, aNNLO, and aN 3 LO for the H + rapidity distribution in H + H − production at 13.6 TeV LHC energy with MSHT20 aN 3 LO pdf, for m H + = 800 GeV.Although the numerical values for the distribution are much smaller than those found with the smaller mass in Fig. 4, the K-factors relative to the LO results in the inset plot show a similar pattern, i.e. they are relatively flat for central values of the charged-Higgs rapidity but begin to grow for rapidities larger than 1.Also, the theoretical uncertainties from scale variation are again essentially the same as those for the total cross section in the rapidity range plotted, and only become bigger at very large rapidities, where the distribution is negligible.Finally, the pdf uncertainties are also similar to those of the total cross section for an 800 GeV mass, around 4%, in the rapidity range shown in the plot, but they increase at larger rapidities.

Conclusions
In this paper we have presented the first calculation of the complete soft-gluon corrections in the production of colorless final states through N 3 LO in single-particle-inclusive kinematics.We have derived the detailed analytical expressions for these corrections and, as a concrete application, we have used our results to study higher-order QCD corrections for the production of a charged Higgs pair (H + H − production) at the LHC via quark-antiquark annihilation in the 2HDM.This calculation is particularly relevant for type II 2HDM, where the value of tan β is constrained by experiments to be small and the quark-antiquark annihilation production mechanism is indeed the dominant channel for such small values of tan β.
We have computed cross sections at aNNLO by adding the second-order soft-plus-virtual QCD corrections to the complete NLO result, and at aN 3 LO by further adding the third-order soft-gluon corrections.Our results show that the NLO corrections increase the LO cross section by 15% to 20% for charged-Higgs masses in the range from 500 to 1500 GeV, with the exact numbers depending on the collision energy and the charged-Higgs mass, while the NNLO soft-plus-virtual corrections provide an additional 5% to 6% increase.The N 3 LO soft-gluon corrections provide a further increase of 0.5% to 0.6%.In general, the soft-gluon corrections are by far the dominant contribution to the higher-order QCD corrections.For instance, the exact NLO cross section differs by only around one percent or less from the one with NLO soft-gluon corrections alone, and the aNNLO result would change by only 2 per mille if we did not include the virtual corrections at NNLO.For the total cross section, we find that the theoretical uncertainties from scale variation get substantially reduced by going to higher orders.
We have also computed aNNLO and aN 3 LO charged-Higgs rapidity distributions at the LHC with 13.6 TeV center-of-mass energy, for two representative masses, namely m H + = 600 GeV and m H + = 800 GeV.We found that the enhancement from the aNNLO soft-plus-virtual corrections is significant while the aN 3 LO soft-gluon contribution is smaller, similar to the total cross section case.The K-factors are relatively flat for central values of the charged-Higgs rapidity, but begin to grow for rapidities larger than one.The theoretical uncertainties from scale variation and from the pdf are basically the same as those of the total cross section for rapidities less than two, where the rapidity distribution is not negligible.
Finally, we want to emphasize that this work opens the way for using soft-gluon corrections through N 3 LO in single-particle-inclusive kinematics to calculate more precise predictions of production rates of colorless final states in hadron-hadron collisions, both in the Standard Model and beyond.

Figure 1 :
Figure 1: LO Feynman diagrams for quark-induced H + H − production at the LHC.

Figure 2 :
Figure2: The total cross sections at LO, NLO, aNNLO, and aN 3 LO for q q → H + H − production in pp collisions at13 TeV energy with µ F = µ R = m H + using MSHT20 aN 3 LO pdf.The inset plot displays the NLO/LO, aNNLO/LO, and aN 3 LO/LO K-factors.

Figure 3 : 3 LOTable 1 :
Figure3: The total cross sections at LO, NLO, aNNLO, and aN 3 LO for q q → H + H − production in pp collisions at 13.6 TeV energy with µ F = µ R = m H + using MSHT20 aN 3 LO pdf.The inset plot displays the NLO/LO, aNNLO/LO, and aN 3 LO/LO K-factors.

Figure 4 :
Figure4: The rapidity distribution of the H + with a mass of 600 GeV at LO, NLO, aNNLO, and aN 3 LO for q q → H + H − production in pp collisions at 13.6 TeV energy with µ F = µ R = m H + using MSHT20 aN 3 LO pdf.The inset plot displays the NLO/LO, aNNLO/LO, and aN 3 LO/LO K-factors.

Figure 5 :
Figure5: The rapidity distribution of the H + with a mass of 800 GeV at LO, NLO, aNNLO, and aN 3 LO for q q → H + H − production in pp collisions at 13.6 TeV energy with µ F = µ R = m H + using MSHT20 aN 3 LO pdf.The inset plot displays the NLO/LO, aNNLO/LO, and aN 3 LO/LO K-factors.
as well as a partonic threshold variable s 4 = s + t 1 + u 1 = (p 2 + p g ) 2 − m 2 , with p g the momentum of an additional gluon in the final state.Near partonic threshold p g → 0 and, thus, s 4 → 0. The soft-gluon corrections appear in the perturbative series as logarithms of s 4 , i.e. [(ln k (s 4 /m 2 ))/s 4 ] + , with 0 ≤ k ≤ 2n − 1 at nth order in the strong coupling, α s . LO.

Table 2 :
The total cross section for pp → H + H − via quark-antiquark annihilation at LO, NLO, aNNLO, and aN 3 LO at the LHC with √ S = 13.6 TeV and MSHT20 aN3LO pdf.The central results are for µ F = µ R = m H + and are shown together with scale and pdf uncertainties.