Higgs and Z boson associated production via gluon fusion in the SM and the 2HDM

We analyse the associated production of Higgs and $Z$ boson via heavy-quark loops at the LHC in the Standard Model and beyond. We first review the main features of the Born $2\to 2$ production, and in particular discuss the high-energy behaviour, angular distributions and $Z$ boson polarisation. We then consider the effects of extra QCD radiation as described by the $2 \to 3$ loop matrix elements, and find that they dominate at high Higgs transverse momentum. We show how merged samples of 0-- and 1--jet multiplicities, matched to a parton shower can provide a reliable description of differential distributions in $ZH$ production. In addition to the Standard Model study, results in a generic two-Higgs-doublet-model are obtained and presented for a set of representative and experimentally viable benchmarks for $Zh^0$, $ZH^0$ and $ZA^0$ production. We observe that various interesting features appear either due to the resonant enhancement of the cross-section or to interference patterns between resonant and non-resonant contributions.


Introduction
Run I of the LHC has brought the discovery [1, 2] of a scalar particle, so far consistent with the Higgs boson predicted by the Standard Model (SM) [3], and it has given the first evidence of the Brout-Englert-Higgs mechanism [4,5] in particle physics. The increased energy and luminosity that will be achieved in Run II at the LHC will allow us to pin down the properties and in particular the strength and form of the interactions of such a boson with all other SM particles. To this aim a vast campaign of measurements of rates and distributions in various production and decay channels is being planned.
Among such processes is the associated production of a Higgs boson together with a vector boson V , either a W ± or a Z, also known as Higgs-strahlung, i.e., at the leading order in QCD, the Drell-Yan production of an off-shell vector boson qq → V * with its subsequent decay V * → V H. While suppressed in the SM with respect to the leading gluon-gluon and vector boson fusion channels, V H production is of phenomenological interest mostly because the presence of the vector boson (and possibly of leptons coming from its decay) in the final state can help to access the large yet challenging H → bb decay mode. For instance, Higgs-strahlung has been the dominant Higgs search mode at the Tevatron [6]. At the LHC, the ATLAS [7] and the CMS [8] collaborations have investigated V H production, with the Higgs boson decaying to a b−quark pair, both reporting small excesses above the background only hypothesis. Searches for Higgs decaying to W + W − [9,10] and to invisible states [11,12] have also been performed by both ATLAS and CMS.
On the theory side, predictions for ZH production are known at NNLO in QCD and at NLO electroweak in EW theory. The NNLO QCD cross section includes the Drell-Yan type terms of O(g 4 α 2 s ) first computed in [13,14]. In addition to Higgs-strahlung, it has been noted that contributions from quark-anti-quark initiated diagrams where the Higgs is emitted from a top quark loop arise at the same order. These diagrams interfere with the LO and NLO Drell-Yan amplitudes and have been computed in [15], where they were found to contribute to the inclusive NNLO cross section at the percent level. Implementations of the NNLO QCD calculations are publicly available in vh@nnlo [16] and HVNNLO [17]. Fully differential NLO QCD and EW results can be obtained with the program HAWK [18,19], while event generation accurate at NLO in QCD (inclusively and for higher jet-multiplicitites), can be nowadays obtained (automatically or semiautomatically) in several frameworks, i.e., MadGraph aMC@NLO [20] / POWHEG [21] + Pythia 8 [22] /Herwig++ [23] and Sherpa [24]. At NNLO, a purely virtual gluon fusion contribution emerges, through the gg → ZH amplitude squared, which at the LHC can be enhanced by the large gluon-gluon luminosity at small Bjorken x. Its contribution to the total cross section has been known for a long time [25,26] and it has been included in the implementations of the NNLO calculations [16,17]. The gluon fusion component is separately gauge invariant, IR and UV finite and accounts for about 10% of the total NNLO cross section at 14 TeV. Being essentially a leading-order contribution, gg → ZH introduces a rather strong scale dependence to the NNLO result, which in turn is known quite precisely. In order to reduce the associated theoretical uncertainty, recently, NLO corrections for the gluon fusion contribution have been estimated by computing them in the infinite top-quark mass limit [27]. The NLO corrections to this process, O(α 3 s ), while formally part of the N 3 LO ZH cross section, are expected to be large, similarly to other gluon fusion processes such as Higgs single or pair production. The computation of the approximate NLO corrections in the infinite top mass limit has confirmed this expectation. The NLO computation in the infinite top mass limit reduces the scale uncertainty by a factor of two, yet the size of the finite top-quark mass effects remains unknown: the exact NLO result requires two-loop multiscale amplitudes whose analytic form is beyond the current advances in the multi-loop technology. In an effort to further reduce the theoretical uncertainties in this process, a soft gluon resummation for the gluon-gluon contribution has been performed in [28] promoting the previous results to NLO+NLL accuracy. We should note here that in contrast with single Higgs production, where the infinite top-quark mass limit provides a good description of the process, and allows the computation of higher order corrections, here, similarly to (yet with even less control than) Higgs pair production, the much higher scales involved make the effective field-theory approach unreliable, especially so at the differential level.
In addition to the SM production mechanism and characteristics, interesting features can be expected from Higgs production in association with a Z boson in beyond the SM scenarios. The Two-Higgs-Doublet-Model (2HDM) is an attractive framework in which Higgs-strahlung can lead to interesting features. First the range of channels is richer: in addition to the production of the light (125 GeV) Higgs boson in association with a Z boson (Zh 0 ), Z associated production of the heavy scalar (ZH 0 ) and pseudoscalar boson (ZA 0 ) are also possible [29]. Experimental searches are already underway to look for signals of these processes, especially in the case where the cross-sections can be enhanced by the resonant production of an intermediate scalar (H 0 or A 0 ) with subsequent decays into Z and a lighter scalar. In particular, CMS has searched for signals of the decay of the pseudoscalar A 0 into a Zh 0 pair [30] and that of the heavy scalar H 0 into a ZA 0 pair [31] and the results have been used to set constraints on the 2HDM parameter space.
So far considerable effort has been devoted to provide accurate total rates in this channel for both the SM and the 2HDM, but accuracy and precision in the differential distributions is also of vital importance. This need becomes more important for experimental analyses which make use of exclusive observables, in order to tame the typically very large QCD backgrounds. Moving in that direction, it has been noted in the literature [29,32] that the gluon induced component can play an important role. The gluon fusion Higgs p T distribution peaks at higher values than the corresponding Drell-Yan one, and therefore its relative contribution becomes more important in boosted Higgs searches, which are preferred experimentally to reduce the backgrounds. The prospects of such searches have improved recently due to progress in jet substructure techniques, after the seminal suggestion in [33].
The aim of this work is to contribute to the understanding of gluon induced ZH production and to improve the predictions for the differential distributions. We consider the 2 → 2 and 2 → 3 matrix elements entering the gluon fusion contribution to Z Higgsstrahlung. We first review the main features of the 2 → 2 ones and then examine the importance of the 2 → 3 contributions. Given the lack of an exact and fully differential NLO computation for this process, we provide a better description of the kinematics for this component by combining the 2 → 2 and 2 → 3 matrix elements in a merged sample, matched to a Parton Shower (PS). This provides a fully exclusive control at the hadron level. A similar approach has been followed for other loop induced processes, such as single Higgs production [34] and Higgs pair production [35,36]. In general, this method provides a better description of the kinematics, yet as the formal accuracy for total rates remains at LO, it is often combined with a normalisation obtained from higher-order computations, when available.
This merging-matching approach makes use of the fact that while tree level fixedorder amplitudes describe reliably the region of hard and well separated jets, the parton shower provides a better description of the soft and collinear regions. Combining the two requires of course a consistent treatment to avoid double-counting, which is achieved by various merging algorithms. Methods that are widely used for tree level merging are CKKW [37,38], CKKW-L [39,40] (and their later improvements [41,42]), and MLM [43]. More recently new methods have been developed to perform the merging at NLO, see for example FxFx [44] and UNLOPS [45], yet not directly applicable to 2 → 2 loop-induced processes at the Born level yet, mainly due the absence of analytic results for the two-loop 2 → 2 matrix elements.
In this work we study gluon induced Higgs-strahlung at the LHC, presenting the first merged-matched results for gg → Zφ, with φ being a generic scalar, by employing the 0 and 1-jet matrix elements for the SM and the 2HDM. This paper is organised as follows. In Section 2, we discuss the process within the SM, first by reviewing the important features coming from the gg → ZH matrix elements. We also consider the behaviour of the 2 → 3 matrix elements, which we then combine with the 2 → 2 ones. We describe our methodology, and present results both at the parton level and after merging and matching to a parton shower. In Section 3, we explore the results of various 2HDM scenarios using the same calculation setup. We draw our conclusions in the final section.

Gluon induced ZH production in the SM
Representative Feynman diagrams contributing to the gg → ZH process in the SM are shown in Fig. 1. Massive fermions, t and b−quarks, run in the box, while all flavours run in the triangle. The contribution of the two light generations to the triangle vanishes as required by the anomaly cancellation. In practice, it is only the axial vector part of the heavy-quark-Z coupling that contributes to the amplitude. The amplitude for this process was first computed in [25,26]. In what follows, we will first review the main features of the 2 → 2 process for gluon induced ZH production before discussing the implications of the 2 → 3 one. A sample of the relevant diagrams contributing to ZHj is shown in Fig. 2. In addition to the gg initial state amplitudes, the qg and qq channels also open up, when an additional jet is allowed. The gg → ZHg amplitudes were used in [46] to calculate the gg part of the ZHj cross-section at the LHC for various jet transverse momentum cuts. In what follows, we will consider these along with the qg and qq diagrams to discuss the behaviour of the 2 → 3 amplitudes and subsequently to obtain a merged sample of 0 and 1-jet multiplicitities.

Calculation setup
In this work, we employ the MadGraph5 aMC@NLO framework [20]. The one-loop amplitudes squared for ZH and ZHj can be obtained with the help of MadLoop [47], which computes one-loop matrix elements using the OPP integrand-reduction method [48] (as implemented in CutTools [49]). A reweighting procedure is then employed to overcome the present limitations concerning event generation for loop-induced processes 1 . A reweighting method has been employed already for a series of processes within the Mad-Graph5 aMC@NLO framework [34,51,52] both at LO and NLO accuracy. This procedure involves generating events through the implementation of a tree-level effective field theory (EFT), in this case obtained by taking the limit of infinite top-quark mass with all other quarks being massless. In practice, a UFO model [53,54] including the effective theory interactions is imported in the simulation framework. After event generation, event weights Loop |/|M 2 EF T |, where |M 2 Loop | represents the numerical amplitude as obtained from MadLoop. In our case, reweighting proves to be efficient in terms of the computational speed, as the loop amplitudes have to be calculated for significantly fewer phase-space points than what would be needed to integrate them directly. Moreover the EFT leads to distributions that are in general harder in the tails, and therefore the EFT events populate regions that are later suppressed by the exact loop matrix elements, resulting to no significant degradation of the statistical uncertainty.

Parton level results
Before proceeding to the technical setup and presenting results of the merging-matching, we consider the salient aspects as observed at the parton level. The findings of this study will reveal some previously unnoticed features of gg → ZH and will act as a motivation to employ a merging-matching procedure in the following section.
In our computation the heavy quark masses are set to: m t =173 GeV and m b =4.75 GeV, while the Higgs mass to m H =125 GeV and the heavy quark Yukawas are given by y q / √ 2 = m q /v. We note here that finite width effects in the propagators of the loops can be taken consistently into account within MadGraph5 aMC@NLO via the implementation of the complex mass scheme [55,56]. The effect of a non-zero top width is shown in Fig. 3, where the matrix element squared for gg → ZH, for 90 0 scattering, is shown as a function of the invariant mass of the ZH system. The correction is more important at the tt threshold, where it reaches 20%. Finally, when integrated over all centre-of-mass energies and scattering angles, we find the top-quark width to modify the gg → ZH cross-section by ∼2% at 14TeV, an effect similar to that observed for single and double Higgs production in [57] and [52], respectively. For the rest of the results presented in this work the width of the top quark is set to zero. An interesting aspect of the gg → ZH matrix element is its angular dependence. While in Fig. 3 we have fixed the scattering angle to 90 0 , in Fig. 4, we show the dependence of the amplitude squared on the centre-of-mass scattering angle, for various values ofŝ. The matrix element starts with no angular dependence at low energies, but varies significantly with the angle at high energies. This angular dependence of the matrix element implies that at high energies, very forward or backward scattering is favoured over 90 0 scattering. This behaviour originates from the interplay between the triangle and the box diagrams, and their respective angular behaviour. As we will also discuss later, box and triangle interfere destructively, with the triangle contribution dominating at low energies. The cancellation becomes nearly exact at high energies, mostly leaving a remainder from the box contribution that is strongly dependent on the scattering angle.
We now proceed to discuss results for the LHC. For these, parton distribution functions (PDFs) are evaluated using the MSTW2008LO set [58] and the central renormalisation and factorisation scales are set to the invariant mass of the ZH system: MadGraph5 aMC@NLO allows the automatic computation of the scale and PDF uncertainties by a reweighting procedure [59]. In our results, scale variations are obtained by varying the scales in the range of µ 0 /2 < µ R,F < 2µ 0 . Table 1 summarises results for the gg → ZH cross section and the "loop-induced" ZHj contribution, originating from the square of the amplitudes shown in Fig. 2. For reference we also include the total pp → ZH cross section at NNLO obtained with vh@nnlo [16]. The NNLO cross section includes the gg → ZH result of the first row, for which excellent agreement has been found between our computation and the result of vh@nnlo. We note that the results shown in the second row of Table 1 for ZHj are obtained using the loop amplitudes shown in Fig. 2. These qg and qq amplitudes can interfere with the Drell-Yan type real emission amplitudes. This interference contribution to the cross section has been  computed in [15] and found to be at the per-mille level. In our computation we use these amplitudes squared, i.e., at O(α 3 s ). It is clear that at this order, other qg loop-induced contributions can enter squared, for example, the set of diagrams where the Z couples to a light quark and the H to a top-quark loop. We have not included these diagrams here, as we consider them of a different origin, but we have checked that their amplitude squared contribution to the cross section is small, below the femtobarn level, and therefore at least one order of magnitude smaller than those in Fig. 2. The interference of this type of diagrams with the Drell-Yan amplitude was computed in [15] and also found to be small.
Given that the ZHj amplitude is divergent in the limit of a collinear or soft jet, we apply a cut on the p T of the jet to obtain finite results. We have set this cut to 30 GeV in Table 1. The 2 → 3 contribution comes mainly from the gg initiated diagrams, with qg giving about 20% of the ZHj cross section. The ZHj contribution is not as suppressed as expected from the extra power of α s , leading to results comparable in size to the gg → ZH cross section. Of course these results are extremely sensitive to the chosen cut for the transverse momentum of the additional parton, as the cross section diverges in the IR limit. Such a problem would not arise in the case of a NLO computation matched to a parton shower, for example with the MC@NLO method [60], in which all divergences are regularised and cancelled for inclusive observables.
The results in Table 1 also demonstrate the problem of large scale uncertainties for the LO gg cross section, that contribute significantly to the total NNLO scale uncertainty. The problem persists also for the loop-induced ZHj contribution. A significant reduction of these intrinsic QCD uncertainties can only be achieved by a complete NLO computation, as discussed in [27].
In addition to the total cross-section results presented in Table 1, interesting observations can be made by studying the differential distributions for the gluon fusion process. We start by presenting distributions for the invariant mass of the ZH system and the Table 1: Cross sections (in fb) for ZH associated production at the LHC at √ s = 8, 13 and 14 TeV. The uncertainties (in percent) refer to scale variations. No cuts are applied to final state particles apart from the jet p T cut in the second row (p j T > 30 GeV) and no Higgs or Z branching ratios are included. The ZHj contribution shown here is obtained from the loop diagrams shown in Fig. 2   transverse momentum of the Higgs in Fig. 5 for the LHC at 14 TeV. These distributions have been shown elsewhere in the literature, for example in [32], yet we consider them here again for completeness. In addition to the gluon fusion results, we also include the NLO Drell-Yan like distributions obtained automatically with MadGraph5 aMC@NLO. We note here that differential NNLO results for the Drell-Yan like contribution can be provided by the code HVNNLO. As is evident from [17], the NNLO computation leads to an overall 20% decrease of the Drell-Yan component but not to any significant shape difference compared to the corresponding NLO one.
The first observation regards the clear presence of the 2m t threshold, at which the gluon fusion amplitude acquires an absorptive part, related to the on-shell gg → tt , tt → ZH scattering, leading to a characteristic rise in the invariant mass distribution. It is evident from Fig. 5 that the gluon fusion component leads to distributions of fundamentally different shape from the Drell-Yan ones and therefore it should be considered in all relevant studies, in particular in the boosted region of p H T > 100 GeV, where its relative importance increases.
In the plots, we decompose the gluon fusion result into the triangle and box components. The two interfere destructively over the whole range of centre-of-mass energies, with the cancellation between the two being nearly exact at high energies. Such a cancellation is due to unitarity: while each of the two diagrams grows with energy the cancellation leads to a well-behaved amplitude at high energies. We stress here that this behaviour of the amplitude is not present in the infinite top mass limit. In this limit, the amplitude for the box diagram vanishes and therefore only the triangle contributes to the amplitude, giving a rather bad approximation of the one-loop amplitude at high energies. In addition, we note that this is a process highly sensitive to the relative phase between the HZZ and Htt couplings. To demonstrate this, in Fig. 5 we show the result obtained by changing the relative sign between the top Yukawa and the HZZ coupling. In pair with other processes where such unitarity cancellations take place, such as H → γγ or pp → tHj [61][62][63], flipping the sign results in an increase in the gluon fusion induced contribution by a factor of five, and much harder distributions as the interference between triangle and box becomes constructive, see Fig. 5. We conclude that, given the difference in the shape as well as the size of the cross section above 2m t , the ZH invariant mass or transverse momentum of the Higgs or the Z could also be used to bound the relative phase between the Higgs couplings to fermions and to vector bosons.
The difference in the p T shape between the Drell-Yan and gluon fusion production persists also in the distribution of the lepton p T coming from the Z decay. Besides, another interesting aspect of the gluon fusion process is that it leads to different angular distributions for the resulting leptons compared to the Drell-Yan component [64]. This is evident from studying the normalised distributions of the angle θ * l , defined as the angle of the lepton between the lepton and Z direction in the Z rest frame is shown in Fig. 6. In the plot, we use the NLO Drell-Yan result, and plot the distributions with and without a cut of 100 GeV on the p T of the Z. The shape of the distribution without any cut, is significantly different, with the tree-level ZH giving a flat distribution while the gluon-fusion one peaks at 90 0 . The shape becomes similar for p Z T > 100 GeV, while a 200 GeV cut (not shown here) completely eliminates the difference. This behaviour is related to the polarisation of the Z which differs between the two production modes. This can be quantified by examin- Process  ing the relevant polarisation fractions in Table 2, as these are defined in [65]. The fact that the Z in gg → ZH is predominantly longitudinal leads to the central peak, while the small difference between f L and f R leads to a very mild asymmetry for qq. Setting a 100 GeV cut on the Z p T changes these values in agreement with the equivalence theorem, i.e., by increasing the longitudinal polarisation fraction. For completeness, we also mention here that the main background for this process, Z + b-jets leads to predominantly left-handed Z bosons [65,66] and therefore to different angular distributions, that could be potentially used as an additional discriminating handle to distinguish signal and background.
Further to the gg → ZH results that we have discussed above, interesting conclusions can be drawn by studying the loop-induced ZHj distributions. We have seen in Table  1 that these contributions are not negligible and their relative importance increases with the centre of mass energy. A complete NLO computation for gg → ZH would be fully inclusive in these contributions but as such a computation is not available yet, we aim to draw some conclusions by studying them independently. For such a study a minimum cut has to be set on the transverse momentum of the additional jet to avoid the divergent soft and collinear limit. We compare the distributions of the invariant mass of ZH and the p H T to those from 2 → 2 amplitudes by varying the p T cut set on the additional jet in Fig. 7. The ZH invariant mass distribution shows that for all values of the jet p T cut the bulk of the cross-section remains close to the 2m t threshold. The characteristic threshold behaviour due to the absorptive part of the amplitude remains visible for all cuts. At high invariant masses, we find that the amplitudes with an extra parton fall more slowly and overtake the 2 → 2 process.
In contrast to the invariant mass distribution where no extreme modification of the shape takes place, the p H T distribution is very much affected. First, we note that the threshold corresponding to m ZH ∼ 2m t is now not visible in the p T distributions for ZHj. The second and more striking observation is that above 300 GeV the 2 → 3 process leads to a much harder p T spectrum compared to the 2 → 2 one. Moreover, for p H T values above 400 GeV all three distributions for ZHj coincide. The fact that in this region the result is insensitive to the jet p T cut implies that hard jet emissions are dominating. This occurs because by allowing the emission of a parton new kinematic configurations open up. In this high p T region, the kinematic configuration in which a soft jet is emitted and the Z and H basically recoil against each other is not the most favourable one. Instead, the configuration in which a hard jet recoils against the H, with the Z remaining rather soft becomes the preferred one. We have explicitly confirmed this behaviour by setting a high cut on the p T of the Higgs, and studying the corresponding jet and Z transverse momentum distributions. A clear preference for the configurations where the jet is hard and the Z is rather soft is found when sufficiently far from the IR divergence. The behaviour of the 2 → 3 amplitudes at high p T can be traced back to the presence of t−channel gluon diagram such as the gg → ZHg one shown in the top right of Fig. 2, which becomes dominant in this region. The same behaviour is displayed by the qg → ZHq contributions, when these are considered separately, as they include diagrams of the same type as shown in the second row of Fig. 2.
In conclusion, we have found that, especially for the transverse momentum distribution of the Higgs, the emission of an additional jet can dramatically modify the shape, due to new allowed kinematic configurations. This effect might prove important in studies involving highly boosted Higgs as discussed for example in [32]. The 2 → 3 matrix elements are important and therefore need to be taken into account for accurate simulations. To combine the two in a consistent way and therefore provide a more realistic picture of the differential distributions, we will resort to merging and matching to a parton shower. In the following section we will discuss how this method allows us to provide more accurate predictions for the distributions.

Merging different jet multiplicities: setup
Given the lack of a complete NLO computation and the relevance of the 2 → 3 matrix elements, the best available procedure to accurately predict the distribution shapes is to employ the Matrix-Element-Parton Shower (ME+PS) procedure. ME+PS schemes allow the consistent combination of matrix elements with different jet multiplicities via their matching to a parton shower. In our study we employ the k T -MLM scheme as implemented in MadGraph5 aMC@NLO. Merged samples are then passed to Pythia 8 [22,67] for PS.
The implementation of MLM merging in MadGraph5 aMC@NLO/Pythia 8 comes in two variants: the traditional k T -MLM and the shower-k T schemes. The two give comparable results as discussed in [68]. In this study we will employ the shower-k T scheme. While this scheme has been used for phenomenological studies with Pythia 6 in the past, see for example [69], we employ the most recent implementation of the scheme in MadGraph5 aMC@NLO combined with Pythia 8. For a detailed description of this approach one can refer to [68], while here we only mention its main points.
In the shower-k T scheme, matrix element events are generated with a minimum separation p T min , between parton and the initial state (iB), and Q cut between final-state partons (ij), defined by the measure: where ∆R 2 ij = 2[cosh(η i −η j )−cos(φ i −φ j )] and p T i , η i and φ i are the transverse momentum, pseudorapidity and azimuthal angle of particle i. Short distance (parton level) events are then passed to Pythia 8 which evolves them down using the p T -ordered shower. In practice, for each event Pythia 8 records the scale of the hardest shower emission: Q P S hardest . This scale is then used to accept or reject the event as follows: for the low multiplicity events, the event is rejected if Q P S hardest > Q cut , while for the highest multiplicity the event is rejected if Q P S hardest > Q M E softest , with Q M E softest being the scale of the softest matrix-element parton in the event. The value of Q cut is selected on a process-by-process basis to ensure that there is a smooth transition between the ME and PS regimes. In practice, this is assessed by examining the differential jet rate distributions which show if the transition is indeed smooth.

Merged sample results for ZH in gluon fusion
Using the setup described in the previous subsection for the merging and matching, we present in this section our merged results for various observables. In our simulations we keep the H and Z stable. For the merging performed here, the shower-k T scheme is used with Q cut = p T min = 30 GeV. We have checked that this choice leads to smooth differential jet rate distributions, and therefore a smooth transition between the ME and PS regimes.
We start by presenting the results for the invariant mass of the ZH system and the p T of the Higgs in Fig. 8, while p ZH T and p j T distributions are shown in Fig. 9. A comparison is made between the gg → ZH sample showered with Pythia 8 and the merged 0 and 1jet matched sample, presented in combination with the uncertainties associated with scale choices for both the factorisation/renormalisation scale of the hard process and the shower starting scale. We set the central value for the renormalisation and factorisation scales to m ZH , as for the parton-level results. The shower starting scale in Pythia 8 can be set to either the kinematical limit (p T = √ŝ 2 ), corresponding to what we refer to as "power"shower or the factorisation scale of each event (m ZH in our case), i.e., "wimpy"-shower. Pythia 8 allows us, for the "wimpy"-shower case, to modify the shower starting scale in the range of 0.5µ F < Q P S < 2µ F . This gives us the possibility to systematically study the dependence of the results on the choice of the shower scale for both the merged and gg → ZH-only samples, as shown by the blue bands in the plots. To study the systematic uncertainties due to renormalisation and factorisation scale variations, we vary the scales between 0.5µ 0 < µ R,F < 2µ 0 , with µ 0 = m ZH . This variation is shown by the yellow bands in the plots (with the central prediction being the "power"-shower result). In the results shown for m ZH and p H T in Fig. 8, we also include the parton-level results for comparison purposes.
First, we notice that not all distributions are sensitive to the procedure of mergingmatching. In particular, the invariant mass of the ZH system shows no shape variation. In this process, we only have initial state radiation and therefore significant changes in the shape are not expected for an observable like m ZH . Other observables, on the other hand, are highly sensitive to the choice of shower parameters. The distributions for the transverse momentum of the Higgs, p H T , but more importantly the transverse momentum of the ZH system, p ZH T , and that of the hardest jet, p j T , which are trivially zero at parton-  TeV. The left column shows the results obtained for the gg → ZH case, with different starting scale for the shower: "wimpy" and "power" shower. The blue band shows the variation of the shower scale for "wimpy" shower in the range 0.5µ F < Q P S < 2µ F , while the yellow bands show the uncertainty associated with a factor of two variation of the renormalisation and factorisation scales with respect to their central value. The right column shows the same results for the merged sample. The green curves in the left column correspond to the parton level results before passing them through Pythia 8.
level, depend strongly on the shower parameters. We first notice that the shower produces a p H T distribution harder than the parton-level one for all shower scale choices. This is related to the harder behaviour of the 2 → 3 distributions discussed earlier.
Another interesting observation to be made is related to the shape changes associated with the shower scale choice. The "power"-shower leads to consistently harder distributions, while the "wimpy"-shower gives softer distributions. The different shower predictions start to diverge in a region correlated with the invariant mass of the ZH system, as this is the factorisation scale which is taken to be the starting scale of the "wimpy"-shower. The shower scale uncertainty bands become wider at larger p T values. This is more evident in the second set of observables, p ZH T and p j T , for which the non-merged predictions can vary by more than one order of magnitude between different shower scale options. At high transverse momentum, the shower uncertainty becomes more important than the intrinsic QCD one associated to the factorisation and renormalisation scale choice for the  hard process. We note here that despite the fact that the factorisation and renormalisation scale uncertainty is large, as evident from the yellow bands, it seems to mainly affect the normalisation of the curves. The advantage of the ME+PS procedure is then made obvious by noticing that the shower scale uncertainty is almost completely eliminated in the merged predictions. For all observables, the shower scale uncertainty bands remain well within the corresponding renormalisation and factorisation scale uncertainty ones, even at high transverse momentum. ME+PS predictions are therefore more accurate/precise and predictive than the parton shower alone as they include the exact 2 → 3 matrix elements. These play an important role in the phase space regions populated by highly boosted objects which is often the case for LHC searches.

Zφ production in the 2HDM
In the previous section we discussed gluon induced ZH production in the SM, employing the ME+PS merging method to improve the accuracy of the predictions for the differential distributions at the LHC. In this section, we will follow a similar approach for a beyond the SM scenario. The case we consider in this work is the 2HDM, as a minimal extension of the SM [70]. The 2HDM extends the scalar sector of the SM by introducing a second SU (2) L doublet Φ 2 , which leads to five physical Higgs bosons, i.e. in the case of CP conservation, the light CP -even one, h 0 , a heavier CP -even one, H 0 , a CP -odd one, A 0 , and two charged Higgs bosons H ± . Assuming no extra sources of CP -violation, seven input parameters fully characterise the model: with the convention 0 ≤ β − α < π (with 0 < β < π/2) fixing the sign of the Higgs coupling to the gauge bosons to be the same as in the SM. Depending on the structure of the Yukawa couplings, two main types of 2HDM setups can be considered: type I where all fermions couple to just one of the Higgs doublets and type II, where up-quarks (down) couple only to Φ 2 (Φ 1 ).
Various theoretical requirements such as stability, perturbativity and unitarity, impose constraints on the 2HDM parameter space. At the same time, electroweak precision measurements and recent LHC Higgs physics results further constrain the parameter space, as discussed in more detail in [51]. Nevertheless, 2HDM scenarios that satisfy these constraints and yet have a significantly different phenomenology than the SM exist and have been studied extensively in the literature. These scenarios arise in the "decoupling" limit [71], i.e., in the limit of cos(β − α) ≪ 1, which ensures that the masses of the additional Higgs bosons lie well above the light-Higgs one. Even in this case, significant shifts from the SM couplings are allowed, in particular in the tanβ ≫ 1 limit, known as "delayed decoupling" [72]. Interestingly, the "decoupling" limit is not the only 2HDM realisation which is consistent with all parameter space constraints. Scenarios with light additional Higgs bosons are also allowed in the so-called "alignment limit" [71].
Some examples of viable 2HDM scenarios will be considered in this section to explore possible 2HDM signatures in Higgs production in association with a Z boson. Interesting features can arise in this process, not only because of possible deviations of the light Higgs couplings from their SM values, but most importantly because of the presence of the heavier states, H 0 and A 0 , which can lead to resonant production of Zφ final states. Three neutral combinations of final states are possible: Zh 0 , ZH 0 and ZA 0 . These 2HDM processes have already been discussed in [29]. Similarly to ZH production in the SM, the production of the Zh 0 and ZH 0 final states can occur through Drell-Yan type diagrams, and in gluongluon fusion. The Drell-Yan like cross sections can be obtained straightforwardly by the appropriate rescaling of the SM cross-sections by the ratio of the g φ ZZ coupling over its SM counterpart, but the situation for the gluon fusion case is more involved. This can be inferred by considering the corresponding Feynman diagrams for the gluon fusion processes, shown in Figs. 10 and 11. The possibility of resonant production depends on the masses of A 0 and H 0 , while interesting interference patterns can arise due to relative sign of the A 0 φZ couplings. For completeness and to facilitate the discussion that follows, the dependance of the relevant Yukawa couplings on the 2HDM parameters is shown in Table 3, for type-I and type-II setups, as rescalings of their SM counterparts. We note that the following Coupling type-I type-IÎ couplings are also relevant for these process (valid for both type-I and type-II setups): withĝ φ V V being the rescaling of the φV V coupling compared to the HV V one in the SM, while the A 0 φZ couplings are proportional to: Figure 10: Representative Feynman diagrams for ZH 0 /Zh 0 production in gluon fusion in the 2HDM. Figure 11: Representative Feynman diagrams for ZA 0 production in gluon fusion in the 2HDM.
We stress here that several studies have been presented in the literature in particular for the A 0 Z process, mostly in the context of the MSSM [73][74][75][76][77]. In the case of the MSSM, there are more constraints on the values of the Higgs couplings, while the 2HDM allows more freedom that can lead to more striking signals. A particularly interesting cosmologically motivated 2HDM scenario leading to a A 0 → ZH 0 signature at the LHC is presented in [78], that finds very good prospects for discovery or exclusion even for the lowluminosity LHC. We also mention that various 2HDM scenarios allow significantly enhanced bottom Yukawas, and the Zφ states can be produced mainly through bb annihilation. This has been extensively discussed in the literature [76,79,80], and in relation with the subtleties of the treatment of the bottom quarks [81]. In this work, we will be focussing on the gluon fusion channel, presenting results for a series of 2HDM benchmarks.

Calculation setup
The calculation setup, regarding the reweighting and the ME+PS merging procedure, follows closely that described in the previous section for the SM. In this section we discuss the details specific to the 2HDM implementation. The computation is performed within the MadGraph5 aMC@NLO framework [20], where MadLoop [47] is used for the computation of the one-loop amplitudes. The 2HDM@NLO model obtained from the package NLOCT [82] is imported into MadGraph5 aMC@NLO for the computation of the 2HDM amplitudes. The model is based on the FeynRules [83] and UFO [53] frameworks. More importantly for our computation, it includes all the necessary UV counterterms and R2 vertices for the MadLoop calculation. The model allows the computation of tree-level and one-loop amplitudes within a completely general 2HDM setup. The 2HDM parameters for the different benchmarks are imported into MadGraph5 aMC@NLO [20] using a parameter card, constructed with the help of an in-house modification of the public calculator 2HDMC [84], in the same way as in [51].
We stress again here that as described in [51] the 2HDM benchmarks to be used here have been constructed in agreement with all up-to-date parameter space constraints, which we have included through an interface of the public tools 2HDMC [84], HiggsBounds [85,86], SuperIso [87,88] and HiggsSignals [89,90] along with additional routines of our own. Three benchmarks will be employed to present the 2HDM results, with the corresponding parameters shown in Table 4. Benchmarks B1 and B2 have been constructed and used already for our study of Higgs pair production in the 2HDM [51], they correspond to B1 and B4 in [51]. Benchmark B3 is a new one designed for this study. Here we briefly mention the main features of each benchmark. For completeness we also show the couplings relevant for gg → Zφ production in Table 5, as rescalings of the SM couplings, similarly to Table 3 Table 3. Yukawa couplings are normalised to their SM counterparts as discussed in the text, while for the A 0 ZH 0 and A 0 ZH 0 couplings we show the proportionality constants of Eq. 3.3.
• Benchmark B3: Another type-II 2HDM scenario with a reversed mass hierarchy between the heavy scalar H 0 and the pseudoscalar A 0 . The small tan β value allows us not to over-suppress theĝ A 0 tt coupling, while theĝ A 0 bb is enhanced. Thanks to the inverted mass hierarchy m h 0 < m A 0 < m H 0 the resonant production of A 0 with a Z boson due to the heavy neutral Higgs decay becomes kinematically allowed.

2HDM results
In this section we present our results for the three 2HDM benchmarks introduced in the previous paragraph. We start by considering the total cross section for each process, which is shown in Table 6. The heavy quark masses are again set to 173 and 4.75 GeV for top and bottom quarks, and the light Higgs mass to 125 GeV. The rest of the calculation details, such as the scale and PDF choices follow closely those of the SM calculation. We note here that where possible, we compared our results with the vh@nnlo version described in [29] and found very good agreement between the two implementations. Before moving to the discussion of some differential results, we first comment on the results in Table 6. First we notice that the cross-section for the Zh 0 process can be significantly enhanced. To be more precise, benchmark B3 leads to a cross section nearly twice the SM prediction, benchmark B1 to a 60% enhancement, while B2 is gives a smaller ∼20% increase. The main source of the increase in the cross-section is the presence of the resonant decay A 0 → Zh 0 , which is kinematically allowed in all three scenarios. The relative change in the Zh 0 cross section is strongly correlated with the mass of the pseudoscalar and the value of the A 0 Zh 0 coupling. We remind ourselves that this coupling is proportional to cos(β − α), i.e., it tends to zero in the alignment limit. For all scenarios considered here, its value remains small as seen in Table 5  to receive extremely large contributions from the resonance. This is in contrast with what we have seen in light Higgs pair production where the resonant decay of the heavy Higgs can lead to an enhancement of up to a factor of 60 for the gg → h 0 h 0 cross section [51].
The most interesting feature of Table 6, is the potential size of the cross section for the ZH 0 process. We find that this can exceed 1 pb when the pseudoscalar A 0 is sufficiently heavy to allow the resonant decay into the heavy Higgs and a Z. This has been noticed and discussed recently in [78], as a signature for a cosmologically motivated 2HDM scenario. It is remarkable that even if the production threshold lies significantly higher, this process can lead to larger cross sections compared to the Zh 0 . This is possible as the relevant coupling, ZH 0 A 0 , as shown in Table 5, is not suppressed by the "SM-like" light Higgs constraints. Despite the fact that the prospects for discovery depend strongly on the resulting decay products of the heavy Higgs, it is worth noting that even in the scenarios where H 0 decays predominantly into bb, the current experimental searches for ZH set a cut on the invariant mass of the bb pair close to the light Higgs mass and would therefore miss this signal. Finally, we note that the ZA 0 production cross section remains very small in the scenarios where the A 0 is heavier than H 0 , but can reach the picobarn level in a scenario such as benchmark B3, as a result of the inverted mass hierarchy.
Further interesting information on these processes can be extracted from the differential distributions. For brevity we present only those for the invariant mass of the system and the transverse momentum of the Higgs, but our setup is fully differential and any distribution can be plotted. We show these in Fig. 12, for the cases in which the cross section is not negligible. The results shown here are obtained with merged samples of 0 and 1-jet matched to Pythia 8 for parton shower, in the same setup as that described in Section 2 for the SM.
For the Zh 0 final state we also show the SM prediction for comparison. Resonance peaks arise in all scenarios for Zh 0 , each time located at the mass of the pseudoscalar A 0 . The sharpness of the peak varies with the mass of A 0 , as heavier A 0 have larger widths going from 0.01 GeV for B3, to 7 GeV for B1 and 35 GeV in B2. We also notice various interesting interference patterns, clearly visible for benchmarks B1 and B2. The A 0 -mediated diagram interferes with the SM-like amplitude, with the interference switching sign at √ŝ = m A 0 . Comparing scenarios B1 and B2, we see that the Zh 0 A 0 couplings have opposite signs and therefore in one case the dip appears right before the resonance peak,  and in the other right after. More subtle features are also visible in the plots away from the resonance peaks. These features can always be traced back to the 2HDM parameters and the value of the relevant couplings as shown in Table 5. One such example is the fact that the B2 m Zh 0 curve lies a bit lower than the SM one in the region below 350 GeV, which can be linked to the enhanced top Yukawa leading to a bigger box contribution. The box is in turn interfering destructively with the triangle leading to a smaller total amplitude for the gg → Zh 0 process.
For the ZH 0 process, only the two benchmarks that give measurable cross sections are shown. The plots shown for this process are dominated by the resonant decay of A 0 . This is more obvious in the B1 curve as the resonant peak is closer to the threshold. Scenario B2 receives some non-negligible off-peak contributions from the Z triangle and box diagrams, which in this case interfere constructively, as the H 0 top Yukawa sign is flipped. For B1 both the top Yukawa and ZZH 0 couplings signs are flipped, therefore the interference between triangle and box is destructive, and the result in the tails away from the resonant peaks, is suppressed compared to B2.
The situation is less complicated for ZA 0 for which in B3 a resonance very close to the m Z + m A 0 threshold dominates the plots, while the cross sections for B1 and B2 are extremely suppressed as no resonant decay is kinematically allowed. Moreover the production of a rather heavy ZA 0 pair probes the gluon luminosity at large partonic x values and is therefore suppressed.

Conclusions
Investigating the nature of the Higgs boson discovered at the LHC is a challenging task. While the results of the measurements undertaken so far show that the 125 GeV scalar agrees well with the SM prediction, there is still room for deviations from the SM and possibly an extended Higgs sector to be discovered at the LHC. The exploration of various Higgs production processes is of vital importance to exclude or confirm a non-minimal Higgs sector.
An important process yet to be precisely measured at the LHC is the associated production of a Higgs with a Z boson. In addition to the Drell-Yan type contributions, this process acquires a gluon fusion component at NNLO, which proves to be of particular importance in the boosted regime. In this work, we have reviewed the main features of the gg → ZH process, both at the matrix-element and cross-section level. We have examined the behaviour and the relative importance of the 2 → 2 and 2 → 3 matrix elements for the gluon induced component. We have found that in the high p T regions the 2 → 3 matrix elements behave in a different way from the 2 → 2 ones and therefore have to be taken into account to provide accurate predictions for the differential distributions. To achieve this, we have combined the two in a consistent way, by merging different jet multiplicity samples and matching them to a parton shower.
Our results have been obtained within the MadGraph5 aMC@NLO framework with the help of Pythia 8 for the parton shower. The ME+PS approach provides a more accurate description of the process compared to the parton shower alone. In particular, it significantly reduces the uncertainty associated with the shower scale choice. For observables such as the transverse momentum of radiated jets in the hard region, the prediction of the parton shower alone can be misleading as here the results are extremely sensitive to the shower parameters. We find that in the merged predictions this sensitivity is almost completely eliminated, with the shower uncertainty remaining well within the intrinsic QCD uncertainty due to the renormalisation and factorisation scale variations.
The reduction of the uncertainties associated with the SM prediction and especially the accurate description of differential distributions is crucial for searches for beyond the SM scenarios. One scenario that the LHC aims to explore is the 2HDM. In this paper, we have also provided predictions for the gluon fusion component of the Zφ associated production in the 2HDM. Following the same setup as in the SM, we have presented our predictions for three representative 2HDM benchmarks. We have considered all three neutral Higgs bosons, presenting results for the cross sections and the differential distributions.
In the production of the light Higgs in association with a Z, large enhancements can be achieved compared to the SM prediction if the resonant decay of the pseudoscalar A 0 is kinematically allowed. Moreover, interference patterns arise between the additional diagrams and the SM-like ones, leading to interesting features in the differential distributions. The resonant production of a H 0 Z pair also becomes important as the H 0 ZA 0 coupling is not suppressed, leading to large cross sections for gg → ZH 0 if the pseudoscalar A 0 is heavier then H 0 . Finally in scenarios where the pseudoscalar A 0 is lighter than the heavy Higgs, gg → H 0 → ZA 0 production is allowed and leads to large cross sections in still-to-be-excluded scenarios.