vh@nnlo-v2: New physics in Higgs Strahlung

Introducing version 2 of the code vh@nnlo, we study the effects of a number of new-physics scenarios on the Higgs-Strahlung process. In particular, the cross section is evaluated within a general 2HDM and the MSSM. While the Drell-Yan-like contributions are consistently taken into account by a simple rescaling of the SM result, the gluon-initiated contribution is supplemented by squark-loop mediated amplitudes, and by the $s$-channel exchange of additional scalars which may lead to conspicuous interference effects. The latter holds as well for bottom-quark initiated Higgs Strahlung, which is also included in the new version of vh@nnlo. Using an orthogonal rotation of the three Higgs CP eigenstates in the 2HDM and the MSSM, vh@nnlo incorporates a simple means of CP mixing in these models. Moreover, the effect of vector-like quarks in the SM on the gluon-initiated contribution can be studied. Beyond concrete models, vh@nnlo allows to include the effect of higher-dimensional operators on the production of CP-even Higgs bosons. Transverse momentum distributions of the final state Higgs boson and invariant mass distributions of the $V\phi$ final state for the gluon- and bottom-quark initiated contributions can be studied. Distributions for the Drell-Yan-like component of Higgs-Strahlung can be included through a link to MCFM. vh@nnlo can also be linked to FeynHiggs and 2HDMC for the calculation of Higgs masses and mixing angles. It can also read these parameters from an SLHA-file as produced by standard spectrum generators. Throughout the manuscript, we highlight new-physics effects in various numerical examples, both at the inclusive level and for distributions.


Introduction
With the ongoing Run II of the Large Hadron Collider (LHC), more and more properties of the Higgs boson H with mass 125 GeV, discovered in 2012 [2,3], are determined with increasing precision [4,5]. An important process in this respect is Higgs Strahlung, i.e. the production of the Higgs boson in association with a gauge boson V , which was used recently to access for the first time the decay rate H → bb [6,7], for example.
At tree-level, Higgs Strahlung is initiated through light quarks (Drell-Yan (DY)-like contribution). At O(α 2 s ), a loop-induced, gluon-initiated contribution enters ZH production which is particularly sensitive to new-physics effects such as modified Yukawa couplings, new particles in the loop, or resonant additional Higgs bosons. In extended Higgs sectors, also the bottom-quark initiated tree-level sub-process may become important for ZH production.
In order to study such new-physics effects at the level of total as well as differential cross sections, we extend the code vh@nnlo [8], which in its previous public release only included Higgs Strahlung within the Standard Model (SM). The main purpose of vh@nnlo so far was the efficient calculation of the total cross section for Higgs Strahlung including next-tonext-to-leading order (NNLO) quantum chromodynamics (QCD) [9,10] and next-to-leading order (NLO) electro-weak effects [11]. The NNLO QCD corrections to the DY-like terms were calculated by employing ZWPROD's implementation of the total NNLO Drell-Yan cross section [12]. Soft-gluon resummation effects are small compared to the NNLO fixed-order result [13], which indicates a good convergence of the perturbative series. The SM electroweak effects to the DY-like terms are implemented in terms of numerical tables, obtained from Ref. [11]. Also included are contributions to DY-like production involving a closed fermion loop, first calculated in Ref. [14] (see Sect. 2.2 for details). The gluon-initiated contribution gg → ZH was originally implemented with the help of FeynArts [15] and FormCalc [16] at the level of squared amplitudes. The current implementation, on the other hand, is at the level of amplitudes: For the SM, the 2-Higgs-Doublet Model (2HDM), and Minimal Supersymmetric Standard Model (MSSM), we use the result of Kniehl and Palisoc [17]; other new-physics effects have been taken into account with the help of FormCalc. For the SM, the calculation of the gluon-initiated contribution is supplemented by the NLO QCD K-factor as described in Ref. [18]. Ref. [19] discussed the resummation of threshold effects for gg → ZH in the SM, which reduce the scale uncertainty substantially; they are not included in the current version of vh@nnlo though.
As already reported in Ref. [20], vh@nnlo has been extended to study Higgs Strahlung pp → Wφ/Zφ in a general 2HDM, one of the simplest extensions of the SM Higgs sector, see Refs. [21][22][23][24][25][26]. Assuming the absence of tree-level flavor-changing neutral currents and CP violation, the two Higgs doublets form three physical, neutral Higgs bosons φ: two CP-even Higgs bosons h and H 0 with m h < m H 0 and one CP-odd Higgs boson A. The different ways to couple the two Higgs doublets to fermions imply four types of 2HDMs (see Ref. [20], for example). A brief report on Higgs Strahlung within the 2HDM is contained in Ref. [27]; more detailed studies can be found in Refs. [20,28]. The DY-like contribution to the cross section in the 2HDM can simply be obtained from the SM result by re-weighting with the corresponding V V φ couplings. The gluon-initiated contribution, on the other hand, also depends on the Yukawa couplings. Moreover, with an increased bottom-quark Yukawa coupling, the bottom-quark initiated contribution, where the final state Higgs φ couples directly to the initial state bottom-quarks, is of relevance. In addition, both latter processes involve contributions from the s-channel exchange of a (resonant or virtual) scalar φ = φ, which may alter the cross section decisively. The program vh@nnlo takes them into account, including all interferences. While these statements also hold for the MSSM, the latter involves additional effects in the gluon-initiated component due to squark-loops [17,[29][30][31][32][33][34], albeit only for ZA production. The current release of vh@nnlo includes those as well. Furthermore, both in the 2HDM and the MSSM, vh@nnlo allows for a mixing of the Higgs CP eigenstates to three neutral mass eigenstates.
The implementation of the gg → Zφ component at the amplitude level clearly facilitates the inclusion of new-physics effects in vh@nnlo. Therefore, aside from the 2HDM and the MSSM, the new version of vh@nnlo currently also includes vector-like quarks (VLQs) and effective Higgs couplings through dimension-6 operators. The latter were checked successfully against earlier work in the framework of MadGraph [35,36]. For the former, we follow the general parametrization of Ref. [37] of the seven relevant representations of the VLQs.
Whereas most of our discussion aims at the prediction of the inclusive cross section for Higgs Strahlung, our implementation of the gg → Zφ and the bb → Zφ components also allows for differential quantities such as the Higgs boson's transverse momentum (p T ), and the ZH invariant mass distribution. Within the SM, fully differential predictions for Higgs Strahlung at NNLO QCD were provided in Ref. [38][39][40]. NLO electro-weak effects were calculated in Ref. [41]. New-physics effects on the boosted regime (p H T 150 GeV) for the gluon-initiated contribution were first discussed in Ref. [20,42,43]. In order to allow for the prediction of kinematical distributions for all sub-contributions to Higgs Strahlung, the new release of vh@nnlo provides a link to MCFM [44][45][46][47].
The paper is organized as follows: We start with an outline of the Higgs-Strahlung process and its partonic sub-processes as they are implemented in vh@nnlo in Sect. 2. Section 3 describes the generic settings of vh@nnlo which are relevant for every run of vh@nnlo. Section 4 focuses on how to obtain kinematic distributions in vh@nnlo. In Sect. 5, we address the newly implemented new-physics scenarios and consider their numerical effects through various examples. Sect. 6 contains our conclusions. Details about the installation and compilation of the code can be found in Appendix A, and options for links to external codes are described in Appendix B.
2 General structure and Standard Model mode 2

.1 Contributions to the cross section
The implementation of the total inclusive cross section for Higgs Strahlung in vh@nnlo has the form where σ Vφ DY denotes the DY-like terms, whose leading order diagram is shown in Fig. 1 (a). By definition, we only count diagrams as DY-like if they can be obtained from Fig. 1 (a) by dressing it with real or virtual gluons and quarks via strong interactions, and possibly crossing them between the initial and final state. Consequently, the structure of the QCD corrections to σ Vφ DY is the same for both Wφ and Zφ production; it is determined by the QCD corrections to the DY process q q → V * , provided through NNLO in Ref. [12] (see Ref. [10] for details). The other terms in Eq. (1) will be discussed in more detail in subsequent parts of this paper.

Standard Model cross section
In the SM mode of vh@nnlo, all terms shown in Eq. (1) can be included. All but the first and the last one only contribute at NNLO QCD and beyond, i.e. at O(α 2 s ), and involve loops of top-or bottom-quarks. Eq. (1) assumes that electro-weak effects δ Vφ EW [11,41] to the DY-like contribution fully factorize from the QCD effects. The term σ VH I is similar for WH and ZH production and involves diagrams where a top-quark loop is inserted into a gluon line of the real and virtual NLO QCD diagrams to pp → V * , from which the Higgs boson is radiated. The effect of these contributions on the inclusive cross section is at the 1-2% level [14]. The contribution σ ZH II , which only exists for ZH production, represents qq-induced diagrams where the Z boson couples to a closed top-or bottom-quark loop. Its contribution on the total rate is even smaller than σ VH I [14]. We will refer to the latter two contributions collectively as σ VH I+II in what follows. Much more relevant for ZH production is the gluon-initiated contribution σ gg→ZH , whose generic LO diagrams are displayed in Figs. 1 (b) and (c). Here, a closed fermion loop is connected to two initial state gluons and the Z boson is attached to this loop, whereas the Higgs boson H is either radiated off the Z boson or off the closed fermion loop, leading to triangle and box diagrams, respectively. It is well known that these two contributions interfere destructively in the SM [48,49]. Ref. [18] provided the NLO QCD K-factor for the gg → ZH sub-process in the heavy-top limit, denoted K ∞ NLO in Eq. (1), which enhances this contribution by roughly a factor of two. Corrections beyond the strict heavy-top limit were calculated in Ref. [50], and Ref. [19] supplied the soft-gluon resummation which leads to a decrease of the residual renormalization scale dependence. As pointed out in Sect. 1, the current implementation of σ gg→Zφ in vh@nnlo uses the results of Ref. [17] and the NLO K-factor of Ref. [18], as well as new-physics amplitudes obtained with FormCalc [16]. Soft-gluon resummation [19] is currently not included.
Eq. (1) involves two bottom-quark initiated contributions. The first one is contained in σ Vφ DY and depends only on the gauge couplings of the bottom quarks. The second one, σ bb→Zφ , is proportional to the bottom Yukawa coupling; sample diagrams which contribute to the latter in the SM are shown in Fig. 2 (a) and (b). Currently, σ bb→Zφ is implemented only at LO QCD. Since we assume the initial-state bottom-quarks to be massless (while keeping the Yukawa coupling non-zero), there is no interference between σ bb→Zφ and σ Vφ DY . The numerical effect of σ bb→Zφ becomes important in certain models beyond the SM (BSM), see below.
In vh@nnlo, each of the terms in Eq. (1) is provided separately for the total inclusive cross section, as will be explained in more detail in Sect. 3.

Generic options of vh@nnlo
In this section we focus on the basic input for a generic vh@nnlo run. The installation of the code is described in Appendix A. The code itself can be downloaded from Ref. [1].
Compared to version 1.0, vh@nnlo has significantly extended its functionality and capability. This implies an extended set of input parameters. The input of vh@nnlo is controlled through SLHA-inspired input blocks (see Refs. [51,52]). We begin with the description of the blocks VERSION, VHATNNLO, 1 ORDER, SCALES, PDFSPEC, and SMPARAMS, which define the central parameters for each run of vh@nnlo and need to be part of all input files. The settings are summarized in Tables 1-7. In these tables and throughout the draft, the notation NAME(i)=k means that the entry i of Block NAME is set to k.
For all input blocks which are required in the SM mode (i.e. VHATNNLO(1)=0), vh@nnlo provides default entries which are automatically set if the corresponding entry in the input file is missing. Their values are specified in Tables 1-7 and 10. In contrast, blocks which specify new physics do not provide default values, except for a few documented cases. The units of the input parameters are given in square brackets in the column "meaning" of these tables. All dimensionful quantities are of data type "real". For input values of data type "integer", the table lists the set of allowed values in column "range". Since such restrictions are much more difficult to define for most of the input parameters of type "real", we refrain from any explicit restrictions in this case and appeal to the user's common sense when setting these parameters.
In order to check the consistency of the input file with the version of vh@nnlo, the user needs to provide the version number in VERSION(1). Accordingly, input files of version 1.0 of vh@nnlo will not run with the current version 2.    (3)=-1/0, the bottom-quark induced contribution bb → Zφ by ORDER(4)=-1/0, and the electroweak correction factor δ VH EW by ORDER(5)=-1/0. We note that the settings ORDER(2)=1 and ORDER (3,5)=0 are only allowed in the SM. The entry ORDER(10)=-1/0 toggles ∆ b resummation in the bottom-Yukawa coupling for tan β-enhanced sbottom contributions, in case the link to FeynHiggs is active. We refer to Sect. 5.2 for details.
By default, the renormalization and factorization scale (µ R , µ F ) is set relative to the invariant mass of the Vφ system through the entries SCALES(1,2), respectively. This means that these scales vary with the integration variableŝ when convolving the partonic cross section with the parton densities. With the exception of the DY-like terms, µ R and µ F can also be fixed in absolute terms to the values given in SCALES (3,4), respectively, by setting SCALES(10)=1. The DY-terms must be switched off in this case, i.e. ORDER(1)=-1.
PDFSPEC(1) sets the name of the employed parton distribution function (PDF), see Table 5.
Since vh@nnlo needs to be linked to LHAPDF [53], all PDF sets of LHAPDF can be loaded. For each run, only one PDF set can be specified, which will be used throughout the calculation of all sub-processes. PDFSPEC(10) fixes the PDF set number. vh@nnlo will automatically use the value of α s (M Z ) which corresponds to the specified PDF set for the calculation of the partonic cross section.
Block MASS defines the masses of the Higgs bosons of the theory by employing the corresponding PDG codes; in the SM mode, only MASS(25) is relevant. If run with FeynHiggs or 2HDMC, this block is overwritten though, see Sect. 5.1 and Sect. 5.2. In this case, the Block ORDER entry default range meaning  Table 4: Defining the unphysical scales. M Vφ is the invariant mass of the Vφ system. output file will also list the total decay widths of the h/H/A boson in MASS(250/350/360) as calculated by these external programs.
The block SMPARAMS defines numerical values for the central SM parameters. SMPARAMS(2) sets the numerical value for Fermi's constant. SMPARAMS(3) allows to provide a value for α s (M Z ) which is independent of the one associated with the PDF set specified in Block PDFSPEC, see below. It will only be used by FeynHiggs or 2HDMC, however; the partonic cross section is always calculated using the value of α s (M Z ) as taken from the employed PDF set. In entries 6,10,5,4,11 of Block SMPARAMS, the masses of the topquark, the bottom-quark (on-shell and MS m b (m b )), the Z-boson and the W -boson are defined, respectively. The on-shell bottom-quark mass M b is used to define the SM Yukawa Block PDFSPEC entry default meaning 1 'PDF4LHC15 nnlo mc.LHgrid' PDF set name (from LHAPDF) 10 0 PDF set number   = M b /v in the gg → Zφ process, and also enters the associated loop integrals. The MS bottom-quark mass, on the other hand, enters as input for FeynHiggs, 2HDMC, and the sbottom-Higgs couplings, see Sects. 5.1 and 5.2. For the Yukawa coupling four-loop renormalization group running. SMPARAMS (12,13) contain the total widths of the gauge bosons. vh@nnlo internally calculates the weak mixing angle from sin 2 θ W = 1 − M 2 W /M 2 Z , and the fine structure constant from α = G F √ 2M 2 W sin 2 θ W /π, which are listed in entries 15 and 1 of Block SMPARAMS in the output file, respectively. Finally, the Cabbibo angle is specified in SMPARAMS (14). It enters the calculation of the DY-like terms and top-quark induced contributions.
Another potential input block available in all models is Block FACTORS which allows to alter the couplings of all involved Higgs bosons to quarks and gauge bosons. The user can either specify common factors identical for all Higgs bosons through entries 1-3 for the couplings to the top quark, the bottom quark and the gauge bosons, respectively. Alternatively, the entries {111, 211, 311}, {112, 212, 312} and {121, 221, 321} separately change the couplings of the light Higgs boson, the heavy Higgs boson and the pseudoscalar, respectively (see Table 8). If not specified, these entries are replaced by their default value (= 1.0). The AhZ and AHZ coupling of the 2HDM and the MSSM (see Eq. (3) below) are Block FACTORS entry default meaning 1 1.0 factor for htt coupling 2 1.0 factor for hbb coupling 3 1.0 factor for hV V coupling 111 entry 1 factor for htt coupling 211 entry 2 factor for hbb coupling 311 entry 3 factor for hV V coupling 121 entry 1 factor for Att coupling 221 entry 2 factor for Abb coupling 321 entry 3 factor for AV V coupling 112 entry 1 factor for Htt coupling 212 entry 2 factor for Hbb coupling 312 entry 3 factor for HV V coupling re-scaled in the same way as the V V H and the V V h couplings, respectively. Note that, since currently vh@nnlo always assumes a vanishing V V A coupling in accordance with a 2HDM with real parameters, FACTORS(321)is ineffective at the moment. The couplings of the CP-even scalars to squarks are currently fixed to their values determined by SUSY and cannot be changed.
Finally, Block VEGAS allows to control the numerical integration parameters, separately for the DY-like, the gluon initiated, and bottom-quark initiated contributions. More details on this block can be found in Appendix B.1.
Before discussing the new-physics input, we shortly explain the output file vh.out of vh@nnlo. It contains all input blocks, partially extended by parameters calculated internally. In addition, Block MASS and Block HIGGSCOUP provide information about the Higgs mass(es), their widths (if calculated by 2HDMC or FeynHiggs) as well as its/their coupling(s) to third generation quarks and gauge bosons. Block SIGMA summarizes the calculated production cross section, split into the individual sub-processes, including the integration errors. Its entries are presented in Table 9. The total cross section presented in SIGMA(1) is calculated according to Eq. (1). Block SIGMA also includes the value of α s (M Z ) associated with the employed PDF set.

Kinematic distributions
Fully differential predictions for the Higgs Strahlung process within the SM through NNLO QCD have been available for some time [38][39][40]. Since recently, they are also publicly available in MCFM [44][45][46][47]. Since the new-physics effects discussed here affect the DY-like component of the Higgs-Strahlung process only through a global factor, we link MCFM to vh@nnlo and rescale this component accordingly. 2 The kinematic distributions of the gg → Zφ and bb → Zφ components are calculated within vh@nnlo at LO QCD, including all available new-physics effects. Here, one can choose between p T or M Vφ distributions, where p T denotes the transverse momentum of the Higgs boson, and M Vφ is the invariant mass of the V φ system.
If the input file contains a Block DISTRIB (see Table 10), vh@nnlo produces a kinematic distribution which is written to the file ptDist 'outputfile' or mhvDist 'outputfile', where outputfile is the name of the output file, see Appendix A for details. The prefix of the file name depends on the distribution type selected in DISTRIB(1). Block DISTRIB also allows to set a minimal and maximal value for the kinematic variable under consideration in DISTRIB(2/3) (histo start/histo end), together with the bin width in DISTRIB(4). Unless vh@nnlo is linked to MCFM, the kinematical distribution only includes  the contributions from gg → Zφ and bb → Zφ (listed separately). We emphasize that, rather than integrating the distributions within each bin, vh@nnlo directly evaluates the distributions dσ/dp T or dσ/dM Vφ for these contributions at the center of the specified bins. The DY-like contribution requires a link to MCFM and, if established, is shown in a third column, potentially rescaled by the couplings of the specified BSM theory. MCFM needs a number of additional input files to be put in the folder MCFM of the vh@nnlo directory. In order to obtain results which are compatible with the settings of vh@nnlo, the MCFM flags removebr and zerowidth are set to true such that the Higgs and the W/Z bosons do not decay.

Beyond the Standard Model
In the following sub-sections we list the different new-physics models and contributions supported by vh@nnlo. A couple of sample diagrams relevant for the gg → Zφ contribution are shown in Fig. 4. We start with the 2HDM in Sect. 5.1 including a discussion of resonant (pseudo)scalars, and continue with the MSSM in Sect. 5.2. In both sub-sections we explain the potential links to external codes for the calculation of a consistent set of model parameters. In Sect. 5.3, we address CP mixing among the three neutral 2HDM scalars of these two models.  in Sect. 5.5 we go beyond a concrete model implementation and describe the incorporation of higher-dimensional operators. Since σ VH I+II as well as the NLO K-factor to gg → ZH are only known for the SM, they need to be switched off whenever BSM effects are requested. Where appropriate, we discuss numerical examples to highlight the new-physics effects in Higgs Strahlung.

2-Higgs-Doublet Model
The program vh@nnlo allows to calculate the production cross section for any of the three electrically neutral Higgs bosons in the 2HDM. The effects which distinguish this case from the SM calculation have been described in detail in Ref. [20]. Therefore, they shall be summarized only briefly here. For once, the couplings of the CP-even Higgs bosons to the weak gauge bosons and the quarks are different from the SM Higgs couplings. While the former are the same for all 2HDMs, the Higgs-quark couplings depend on the specific realization of the Z 2 symmetry which is imposed in order to avoid tree-level flavor-changing neutral currents. Specifically, for the gauge-Higgs couplings we have where α parametrizes the mixing of the CP-even Higgs bosons from the isospin to the mass eigenstates, and tan β is the ratio of the vacuum-expectation values of the two Higgs doublets. The quark-Higgs couplings are summarized in Table 11. In types III and IV, the quark-Higgs couplings equal the ones of types I and II, respectively. Only the couplings to leptons differ; this affects the total widths of the Higgs bosons, and thus indirectly the Higgs-Strahlung cross section through diagrams with Higgs bosons in the s-channel.
Type I cos α/ sin β sin α/ sin β cot β cos α/ sin β sin α/ sin β − cot β Type II cos α/ sin β sin α/ sin β cot β sin α/ cos β cos α/ cos β tan β    Table 11 also show the couplings of the CP-odd Higgs boson A to quarks. Its coupling to the weak gauge bosons vanishes, but there are tri-linear couplings to Zh and ZH, also included in Fig. 5, where These latter couplings imply new contributions to the Higgs-Strahlung process, with virtual Higgs bosons in the s-channel. When Higgs-quark couplings are neglected for the first two quark generations, such contributions can appear only in the partonic processes bb → Zφ and gg → Zφ. Their numerical effects have been studied in detail in Ref. [20]. They can be huge, for example for gg → A → Zh in the wrong-sign limit (g h d = −1) of the 2HDM [54]. We note that using Block FACTORS, see Sect. 3, implies that the relative factors of g Ah Z and g AH Z are set to the ones of g H V V and g h V V , respectively. If the mass M φ of the s-channel Higgs boson φ is larger than the threshold for the Zφ production process under consideration, the integration over the parton density functions would encounter a pole atŝ = M 2 φ which, however, is regulated by the width Γ φ of φ through the replacement 1 The numerical value for Γ φ is required by vh@nnlo as input. If FeynHiggs or 2HDMC are linked to vh@nnlo, Γ φ will be taken from their output though. The same effect occurs in the bottom-quark initiated contribution σ bb→Zφ due to Fig. 2 (c). We emphasize that our implementation allows to take into account the interferences of the s-channel diagrams gg/bb → A → Zh or gg/bb → {h, H} → ZA with all other Feynman diagrams. Until now, searches for heavy pseudoscalars decaying to Zh have applied the narrow-width approximation, which means that such interference effects have been neglected. This allowed to include higher-order corrections to the production mechanisms gg/bb → A, as implemented in SusHi [55,56], which can also be linked to 2HDMC to obtain the branching ratio A → Zh. However, it is well-known from decays of heavy scalars into a pair of gauge bosons VV [57,58] or a pair of Higgs bosons hh [59] that interference effects in a 2HDM can be large. Thanks to the option to obtain the invariant mass distribution M Zφ , vh@nnlo allows to study the shape of the interferences around the mass of the intermediate resonance. We will show an example at the end of this section.
One option to evaluate the cross section in the 2HDM with vh@nnlo is to define the three neutral Higgs boson masses and their total decay widths in Block MASS, see Table 12.
The other relevant 2HDM parameters need to be provided in Block 2HDM, also described in Table 12. Recall that various contributions to the cross section can be switched on and off in Block ORDER, see Table 3. Table 11 may serve as definition of tan β and the CP-even Higgs mixing angle α (up to shifts of α by ±π, which leave the cross section invariant). As explained before, 2HDM types I and III as well as II and IV have identical quark-Higgs couplings. Thus, using Block 2HDM will reveal identical results for 2HDM(2)=1 and 3 as well as =2 and 4, provided the total Higgs widths are left unchanged.
Another option to define the 2HDM input parameters is through the link of vh@nnlo to 2HDMC [60,61] (see Appendix B.3 how to establish the link). If vh@nnlo finds the Block 2HDMC in its input file, it will read 2HDMC input from this block, and ignore the block 2HDM. Linking to 2HDMC gives you access to the features of this program within vh@nnlo, meaning that you can use different parameterizations of the 2HDM, check for theoretical consistency of the 2HDM, and get consistent numerical values for the decay widths of the Higgs particles.
The program 2HDMC produces its own additional output file, named 2HDMC.out, which contains all relevant information about the Higgs spectrum and decay widths. The definition of the entries in Block 2HDMC is described in Table 13. The specific 2HDM parameterization is set by 2HDMC(1), which currently can assume the values 1, 2, and 3, entry  default  meaning  25 125.

Block MASS
see text CP even Higgs mixing angle α corresponding to the "λ basis", the "physical basis", and the "H 2 basis", respectively (see Refs. [60,62] for the definition of these bases If the user works with the convention 0 ≤ β − α ≤ π, such that cos(β − α) covers the whole range [−1, 1] and sin(β − α) ∈ [0, 1], then for values between π 2 ≤ β − α ≤ π, the 2HDMC convention is obtained by shifting β − α by −π. The other parameters are in principle not constrained, but 2HDMC can check whether a parameter point obeys stability of the vacuum, unitarity and perturbativity requirements. If 2HDMC(10)=1, vh@nnlo will terminate if a 2HDM parameter point does not obey all of these criteria. This is quite useful for large parameter scans. In contrast to the run with Block 2HDM, the link to 2HDMC will reveal differences among the 2HDM types I and III as well as II and IV. This is due to the fact that the decay widths obtained from 2HDMC depend on all couplings, including the Higgs-lepton couplings.
A thorough discussion of Higgs Strahlung in the 2HDM can be found in Ref. [20]. We therefore focus only on a single numerical example which shows the interference effects of intermediate (pseudo)scalars with non-resonant diagrams as discussed above. In Fig. 6, we exemplify the two cross sections gg → Zh and bb → Zh as a function of the invariant mass of the Zh system around the pseudoscalar mass m A = 500 GeV. The light Higgs mass is {0,1} ignore theory inconsistencies: {yes, no} 2HDMC input key entry meaning chosen to be m h = 125 GeV, and the other relevant parameters are fixed to tan β = 10 and sin(β − α) = 0.999. The total pseudoscalar width obtained by 2HDMC is Γ A = 1.13 GeV. In this region of the parameter space, the "signal" process gg → A → Zh yields a small cross section and induces large interference effects with all other diagrams as shown in Fig. 6. For bb → A → Zh, the effect is not as pronounced, since also the contribution of the non-resonant diagrams is small. We leave a thorough analysis of such interference effects to future work. A generic study using vh@nnlo was already carried out in Ref. [63]. An effective coupling of the pseudoscalar to two gluons through the operator L 6 , see Sect. 5.5, was employed therein.

Minimal Supersymmetric Standard Model
Concerning the quark-loop contributions to the gg → Zφ process, the MSSM amplitudes are identical to those of a type-II 2HDM. As explained in Ref. [17], the squark-loop contributions to CP-even Higgs production vanish, and of the squark-loop contributions to CP-odd Higgs production, only triangle-type diagrams are non-vanishing. We compared the results for the quark-and the squark-loop diagrams for both CP-even and -odd Higgs production from Ref. [17] to our own calculation of the gg → Zφ cross section based on FormCalc [16] and found agreement after fixing some typos in Ref. [17].
Similar to the 2HDM, the parameters for the MSSM mode VHATNNLO(1)=1 can be defined in several ways. The first one is through the blocks 2HDM, MASS, and SQUARK. The block SQUARK is only required for CP-odd Higgs production, i.e. for VHATNNLO(2)=21. As described in detail in Table 14, its entries define the masses and mixing of stops and sbottoms. For the squark mixing angle θ q (q ∈ {t, b}) both conventions cos θ q ∈ [−1, 1] and sin θ q ∈ [0, 1] or vice versa are possible. If sin θ q is not specified, the convention with cos θ q ∈ [−1, 1] is assumed, and sin θ q ≥ 0 is obtained from 1 − cos 2 θ q . Since we work at leading order for what concerns squark effects, the renormalization scheme of these parameters is arbitrary. Furthermore, vh@nnlo does not restrict the range of values for A t , A b and µ. However, we urge the user to choose these parameters in accordance to the squark masses and mixings, which limits their ranges. The Feynman rules for the squark couplings can be found in Ref. [55]. They involve also the quark masses m b and m t , for which we insert the MS value m b (m b ) provided in SMPARAMS (5), and the on-shell value m t of SMPARAMS (6), respectively.
Another way to provide the MSSM parameters is through the link of vh@nnlo to the code FeynHiggs [64][65][66][67][68][69][70] (see Appendix B.3 on how to establish this link). If the input file contains the Block FEYNHIGGS, vh@nnlo will ignore the blocks 2HDM, MASS, and SQUARK, but rather call the program FeynHiggs to determine the spectrum and the couplings of the specific MSSM parameter point. A number of flags for FeynHiggs can be accessed through the block FEYNHIGGSFLAGS. A detailed description of the two blocks is given in Table 15 and 16. 3 The blocks associated with the link to FeynHiggs generally understand all input parameters as on-shell. Accordingly, stop and sbottom masses returned from FeynHiggs are renormalized on-shell. In case the link to FeynHiggs is active, vh@nnlo also allows to incorporate the resummation of tan β-enhanced sbottom contributions in the bottom-quark Yukawa coupling, known as ∆ b resummation [71][72][73][74][75][76]. Setting ORDER(10)=1 will activate the so-called "full resummation" as presented in Eqs. (16c)-(16e) of Ref. [55]. Note that one may also define an independent value of ∆ b by adjusting the bottom-quark Yukawa coupling through Block FACTORS, see Table 8.
A third option to provide the relevant MSSM input is through an additional SLHA-style input file, referred to as "spectrum file" in what follows, which contains the MSSM particle spectrum. The relevant content and a list of potential spectrum generators providing a spectrum file are given in Appendix B.3. This mode of operating vh@nnlo is activated by including a Block SPECTRUMFILE in the actual input file of vh@nnlo, and specifying the path to the spectrum file in SPECTRUMFILE(1). Apart from that, the vh@nnlo input file only requires the generic blocks described in Sect. 3    We present a numerical example for light Higgs production in association with a Z boson in the MSSM in Fig. 7 for the 13 TeV LHC. The Higgs mass and mixing is obtained from to Zh production at order O(α 2 s ), the features correspond to the well-known decoupling behavior for the light Higgs boson. Since all SUSY states are above 1 TeV, they do not enter the pseudoscalar total decay width, which is relevant for the sub-processes gg → Zh and bb → Zh. Only for large enough values of m A and tan β, a light Higgs mass m h in the range of (125 ± 3) GeV is reached through higher-order corrections to the light Higgs boson mass, which at tree-level is bound to be below the Z-boson mass m Z . This is the decoupling region, i.e. all light Higgs boson couplings to SM particles are close to their SM value. This explains the behavior of the DY-like component to Higgs Strahlung, shown in Fig. 7 (b), since it directly follows the squared coupling (g h V V ) 2 of the light Higgs boson h to gauge bosons V ∈ {W, Z}. All components of the cross sections in Fig. 7  1, so that the contribution from triangle diagrams is small, whereas the box diagrams proportional to the Higgs Yukawa couplings depend on tan β. In fact, the minimal cross section is obtained for intermediate   Finally, we present the component bb → Zh in Fig. 7 (d). The MSSM cross section exceeds the SM cross section (which is tiny) by one to two orders of magnitude almost throughout the entire m A -tan β-plane. For m A 150 GeV, i.e. outside the decoupling region, the light-Higgs-boson bottom-quark Yukawa coupling g h b is strongly dependent on tan β, which explains the steep increase of the cross section with tan β. Relevant Feynman diagrams are the t-and u-channel contributions shown in Fig. 2 (a) and (b). For m A 200 GeV, those Feynman diagrams are less relevant, since the light-Higgs-boson bottom-quark Yukawa coupling g h b approaches its SM value-though slightly slower than the top-quark Yukawa coupling g h t . On the other hand, the pseudoscalar contribution bb → A → Zh is manifest (see Fig. 2 (c)), since the pseudoscalar can be resonant. Here the strong enhancement of the pseudoscalar bottom-quark Yukawa coupling g A b with increasing tan β leads to observable differences also at high values of tan β despite the fact that g Ah Z almost vanishes, i.e. decoupling is reached. We note that for this example we did not include ∆ b resummation.
Before we close our discussion on the MSSM, we show p T distributions and M Zφ distributions for the gluon-induced production of a pseudoscalar with a mass of 500 GeV to  Table 17.
demonstrate the influence of squark contributions. The cross sections are produced for the 13 TeV LHC. We choose tan β = 10 and µ = 1 TeV and pick six choices of squark masses shown in Table 17. Again we link to FeynHiggs version 2.13.0 to obtain the light Higgs mass as well as the squark masses. Only scenarios 1 and 3 yield a light Higgs mass in accordance with experimental data, i.e. within 125 ± 3 GeV. Scenario 6 yields the lowest mass with 110 GeV. Such scenarios can be considered pedagogical examples to demonstrate squark mass thresholds and the influence of X t . The process gg → ZA also has a light and heavy Higgs contribution gg → {h, H} → ZA, where the former vanishes for exact decoupling g h V V = 1, i.e. g Ah Z = 0. The squark contributions in exact decoupling are thus relevant only for gg → H → ZA. In Fig. 8 (a) and (b) we first show the influence of the choice of tan β in a 2HDM of type II and compare it with the distribution of a "SM Higgs" with the same mass of 500 GeV. Also the influence of gg → {h, H} → ZA can be inferred from the difference between sin(β − α) = 1 and sin(β − α) = 0.999. The choice of tan β = 10 minimizes the cross section, since again g A t is reduced, while g A b is only moderate. Therefore one expects large relative squark effects in this region of the parameter space, which we show in Fig. 8 (c) and (d). Indeed, in the invariant mass distribution, the squark mass thresholds are clearly visible in scenarios 5+6. Given that the partonic center-of-mass energy √ŝ has to fulfillŝ ≥ m 2 , also the p T distribution in scenario 5 probes the squark mass thresholds 2mt 1 and 2mt 2 at around p T ≈ 320 − 330 GeV and p T ≈ 660 − 670 GeV, respectively. Also the inclusive cross section in scenario 5 is enhanced by a factor of more than 20, due to resonant squark contributions. It is remarkable that even larger squark masses still lead to a significant deformation at low p T and M Zφ of the distributions, in particular in case of large squark mixing X t = 0. However, the cross sections are very small and hard to access experimentally.

CP mixing among Higgs bosons
By default, vh@nnlo assumes the Higgs CP eigenstates (h, H, A) of the 2HDM and the MSSM to be identical to the mass eigenstates (h 1 , h 2 , h 3 ). However, it does provide the option for a naive mixing according to where the orthogonal rotation matrix R is parametrized in the form with c ij = cos θ ij and s ij = sin θ ij . The three independent angles θ 12 , θ 13 , and θ 23 can be specified in the input file. For all but the gg → Zφ and the bb → Zφ contributions, this Block HCPMIX entry default meaning 1 0. simply results in a modification of the V V φ couplings according to where R ih , R iH , R iA refer to the first, second, and third column of the matrix R, respectively, the g φ V V (V ∈ {W, Z}) are the couplings of the specified model without CP mixing (see Eq. (2) and Fig. 5), and g i V V is the coupling constant for the h i V V vertex. Note that, up to the sign, θ 12 has the same meaning as α, defined in Eq. (2) for the 2HDM, since both angles determine the rotation among the two CP-even Higgs eigenstates. We keep this redundancy in vh@nnlo for the sake of a transparent implementation of the mixing.
Finally, we note that the rotation defined through Block HCPMIX is applied after a possible rescaling of the couplings by the entries in Block FACTORS (see Table 8 in Sect. 3).
The gg → Zφ component of the partonic cross section for a specific φ involves also the other two Higgs bosons through s-channel exchange. vh@nnlo combines the corresponding helicity amplitudes by assuming a CP-even and a CP-odd coupling of each Higgs mass eigenstate to the top and the bottom quarks, according to and using the appropriate h i h j Z couplings where g φ q and g Aφ Z are again the qqφ and ZAφ coupling constants of the specified model without CP mixing (see Table 11, Eq. (3), and Fig. 5). Currently, the bb → Zφ component cannot be evaluated with CP mixing within vh@nnlo, so ORDER(4) needs to be set to -1 when CP mixing is switched on. Mixing is activated by including the Block HCPMIX in the input file, where the mixing angles are specified, see Table 18. In order to obtain the cross section for h 1 , h 2 , or h 3 production, one needs to set VHATNNLO(2) to 11, 12, or 21, respectively. Similarly, the masses m h 1 , m h 2 , and m h 3 are defined through MASS (25,35,36), respectively.
Instead of a consistent implementation of CP violation in the MSSM through complex parameters, vh@nnlo simply assumes that the quark and vector-boson couplings of the Higgs are modified as described above, and in addition, the Higgs-squark couplings are affected by an analogous rotation: where g q,φ kl is the coupling constant of φ to the squarksq k andq l in the CP-conserving MSSM, see Ref. [55]. This approach misses contributions due to CP violation in the squark sector, which induces new couplings of the three neutral eigenstates {h, H, A} to squarks.
As an example, Fig. 9 shows the dependence of the gg → Zφ contribution to the cross section on one of the relevant mixing parameters successively for φ = φ 1 , φ 2 , and φ 3 . No mixing corresponds to a type-II 2HDM with tan β = 20 and sin(β − α) = 1 in this figure. The Higgs boson masses are m h 1 = 125 GeV and m h 2 = m h 3 = 300 GeV. As the mixing angle θ ij increases, the produced Higgs boson h i gradually assumes the CP properties and the couplings of what is h j in the unmixed case, while the mass remains unchanged.

Vector-like quarks
Vector-like quarks (VLQs) impact the gluon-initiated contribution to Higgs Strahlung. We adopt the general parametrization of Ref. [37] for adding a single multiplet of vector-like quarks (VLQs) to the SM Lagrangian. There are seven possible representations which can mix with the SM quarks through Yukawa interactions [78,79]. They all involve vector-like quarks T and/or B with electric charges +2/3 and −1/3, respectively, and possibly one of X and Y , with electric charges +5/3 and −4/3. In this way, one arrives at two singlets, three doublets, and two triplets: phase φ q : where the superscript "0" denotes the weak eigenstates in order to distinguish them from the mass eigenstates. This parametrization of the mixing applies to (q, Q) = (b, B) and (q, Q) = (t, T ), and separately to their left-and right-handed components. The left-and right-handed mixing angles are related by for singlets and triplets, and The phases φ q are the same for the left-and right-handed components. For the triplet representations, one finds an additional constraint of the form where a = √ 2 for the triplet (X, T, B), and a = 1/ √ 2 for the triplet (T, B, Y ). Since the T and B quark share the same mass term in the Lagrangian, there is another constraint on the mass splitting of the form Without mixing, the VLQs as introduced here would not couple to the Higgs boson, and their couplings to the Z boson would be purely vector like. This means that they would not contribute to the gg → ZH process at LO due to charge conjugation invariance 4 (Furry's theorem). In fact, the LO amplitude for gg → ZH does not involve X and Y for this reason.
However, aside from modified ttZ, bbZ, ttH, and bbH coupling constants with respect to their SM values, the mixing in the top and bottom sector leads to tT Z, bBZ, tT H, and bBH vertices, an axial-vector component of the TT Z and BBZ coupling, and nonvanishing BBH and TT H couplings. The LO effects of the VLQs on the gg → ZH process are therefore as follows: (i) a modification of the SM box contributions; (ii) additional box contributions involving one to four T or B propagators, see Fig. 10 Figure 11: Feynman rules for the Higgs and Z couplings due to a vector-like T quark calculated with FeynRules [80,81].
to the vertices of two T quarks with a Higgs or Z, and a single T quark with a top quark, a Higgs or a Z boson are listed in Fig. 11 with generic couplings. Here, Q t is the electric charge of the t (and the T ) quark, X L/R T T and X L/R tt are the left-and right-handed couplings of the T and t quark to the Z boson, respectively; Y T T and Y tt are the Yukawa couplings of the T and t quark; and X L/R tT and Y L/R tT are the left-and right-handed couplings of a T and t quark to a Z boson and a Higgs. Note that vertices for a vector-like B quark can be obtained by introducing an additional minus sign to the couplings and changing T → B, t → b. The generic couplings assume values which are specific to the various representations. A full list can be found in Ref. [37]. The couplings of VLQs to gluons are the same as for SM quarks. The amplitudes including VLQs have been evaluated with the help of FormCalc.
The input block for VLQs in vh@nnlo is described in Table 19. VLQ(1) fixes the representation for the VLQs according to Eq. (11); vh@nnlo sets the generic couplings discussed above according to this representation; VLQ(2) and VLQ(3) set the mass of the T and B quark, respectively; the left-and right-handed mixing angles for the bottom sector are set in VLQ(4) and VLQ (5), respectively, and in VLQ (7) and VLQ(8) for the top sector; finally, the phases for the bottom and top mixing can be set in VLQ(6) and VLQ (9). Note that there are conditions like Eq. (15) which constrain these parameters; the input will be checked for consistency by vh@nnlo. In fact, most conveniently the user may provide only an independent subset of the parameters, and let vh@nnlo determine the remaining parameters from the theoretical constraints described above. Specifically, for the (T B) doublet, it is sufficient to provide one of the following sets: 1. m B , m T , and one mixing angle; 2. m B , (θ b L or θ b R ) and θ t R ; 3. m T , (θ t L or θ t R ) and θ b R .
For the triplets (T BY ) and (XT B), possible input sets are: Note, however, that the left-or right-handed mixing angle will always be calculated for any representation, if only one of these mixings and the corresponding mass of the vectorlike quark is set as input. The complete set of parameters, including the coupling factors appearing in Fig. 11, are listed in Block VLQ in the output file.
Block VLQ entry range To study the impact of VLQs on the total cross section, we define the ratios R θ and R φ , In R θ we normalize the total cross section of gg → ZH with VLQs at a given value of θ t/b L/R to the total cross section of the gluon induced contribution of the SM. Moreover, we restrict ourselves to cases with φ t/b = 0. For R φ we choose θ t/b L/R = 0 since otherwise the VLQ will have no effect on the cross section. Thus, instead of normalizing this result to the SM alone, we add the contributions of the VLQs for φ t/b = 0 in the denominator. Additionally, for R φ we set φ t = φ b . We choose the masses and the dominant mixing angle of the vector-like T partner as our input scheme for singlets and doublets, and m T with θ t L as the input for the triplets. Fig. 12 shows these ratios for four different representations of VLQs by varying θ t L (for singlets and triplets) or θ t R for doublets in (a) and φ t/b in (b). The other settings are described in the caption of the figure. We first discuss Fig. 12 (a) in detail: All curves start at the SM cross section for zero mixing angle, where the VLQs have no couplings to the SM particles. For the vector-like top singlet (T), the total cross section for gg → ZH steadily increases with the mixing angle, until it reaches a maximum of about 1.5 times the SM cross section at maximal mixing, i.e. θ t L = π 2 , where the couplings of the t and the T are exactly interchanged w.r.t. the SM. Therefore, the cross section at θ t L = π 2 is the SM value for a hypothetical top quark mass of 600 GeV. This is also true for the (XT ) doublet, where the contributing diagrams are the same as for the (T ) singlet. Since the couplings are different, however, and the dominant mixing is the right-handed one, the intermediate evolution of the cross section differs from the (T ) singlet. It reaches almost ∼ 200% of the SM cross section at about θ t R = 1.1, which is mainly due to the enhancement of the coupling to t in this region (as can be seen in Fig. 13).
For the (T B) doublet, the behavior is strikingly different. While it resembles the behavior of a (T ) singlet at low mixing angles, which is due to the dominance of the top quark in this regime, it quickly drops down to very small values as the mixing increases towards π 4 , before it nearly vanishes at maximal mixing, where T and B assume the couplings of their SM partners, while those of t and b vanish. Although both VLQs carry the same mass, at maximal mixing their couplings to the SM bosons slightly differ due to Eq. (14). If they coupled equally, their amplitudes would be of equal magnitude but of opposite sign for maximal mixing, such that they would completely cancel in the cross section. The same behavior would be observed in the SM, if the b and the t quark had identical masses.
Finally, the contributing diagrams for the (T BY ) triplet are the same as for the (T B) doublet, but with different couplings. We therefore do not see the same cancellation as for the (T B) doublet, which is even more suppressed since only at minimal and maximal mixing the T and B masses are degenerate. Instead, one observes similar features as in (XT ), with the maximum shifted towards θ t R = 1.1. At maximal mixing, the cross section again assumes the value found for (T ) and (XT ). In this case, however, all SM and T quark contributions vanish, and the B quark is the only contribution left. Since for maximal mixing one finds m B = m T = 600 GeV and the couplings are the same as for the (T ) singlet, their results coincide.
In Fig. 12 (b) we can observe the impact of φ t/b for the four different representations. The mixing angle of the top and bottom sector is set to θ t/b L = π 16 , where in the case of the triplet θ b L and m B are calculated according to Eq. (15) and Eq. (16b). Since the phase φ t/b only appears in mixed couplings of quarks, vector-like-quarks and SM bosons, one can probe the magnitude of the contributions in which these mixed vertices appear. For the case of gg → ZH this only affects box diagrams. Here, we set φ t = φ b for doublets and triplets which contain both T and B quarks. The shape is determined by the dependence on exp(iφ t/b ) in the couplings, and varying the phase can lead to an enhancement of up to 50% for the (XT ) doublet and a reduction of −30% for the (T B) doublet. For the studied (T ) singlet and the (T BY ) triplet the effect of the phases is rather mild.
Kinematic distributions are typically much more sensitive to effects of new physics than the total cross section. Thus we study the p T and M ZH distributions with an additional multiplet. Again we emphasize that p T refers to the transverse momentum of the Higgs boson. For the (T ), (XT ) and (T B) multiplets, we fix m T = m B = 600 GeV, set the mixing angle such that we are in a maximum of Fig. 12  input parameters with vh@nnlo. For the (T BY ) triplet, we only fix m T = 600 GeV and θ t L = 1.06. Moreover, for all multiplets we set φ t/b = 0. In Fig. 13 we show (a) the p T and (b) the M ZH distributions of the gluon initiated process in comparison to the SM expectation. In the p T spectra we see an overall enhancement of the SM distribution which is more distinct in the boosted regime, i.e. p T 200 GeV. The threshold of the top quark at p T ∼ 150 GeV and the threshold of the VLQs at about p T ∼ 590 GeV is visible for all representations and more pronounced for (T ) and (XT ) than the other studied cases.
The invariant mass is more sensitive to threshold effects than the transverse momentum, thus showing distinct peaks at M ZH ≈ 2m t and at twice the VLQ masses. It shares the same features as the p T distributions, i.e. the VLQ contributions lead to an overall enhancement in comparison to the SM, which is more distinct in the boosted regime. Since the VLQ masses are not degenerate for the scenario of the (T BY ) triplet studied here, we observe in total three thresholds, one for the top quark, one for the B quark at about 900 GeV, and one for the T quark at about 1.2 TeV.

Higher-dimensional operators
Finally, we describe the implementation of higher-dimensional operators, which is independent of a concrete model implementation. In the SM, our setup allows to set bounds on the implemented higher-dimensional operators. On the other hand, our implementation can also be used in the MSSM or the 2HDM for the production of a CP-even Higgs boson and a Z boson, where additionally one dimension-5 operator coupling two gluons to a pseudoscalar can be taken into account. In our choice of the dimension-6 operators, we adopt the so-called gauge basis of Refs. [82,83] with their corresponding FeynRules implementation [84], but restrict ourselves to the operators which involve third-generation quarks. We allow these operators to be added to any implemented model, i.e. the SM, the 2HDM and the MSSM, but note that we only implemented operators relevant for the production of the light SM-like CP-even Higgs boson gg → Zh. The effective Lagrangian where the sum is over the following higher-dimensional operators: Herein, we define the dual gluon field-strength tensorG a µν = 1 2 µναβ G a,αβ ( 0123 = +1), Eq. (19) is set up as an extension of the SM Lagrangian, which means that, in the SM mode of vh@nnlo, Φ is the SM Higgs doublet. However, vh@nnlo also allows to use Eq. (19) for CP-even Higgs production in BSM models, in which case the same Feynman rules for Eq. (19) as those derived in the SM are assumed, with the Higgs particle being interpreted as the one specified in VHATNNLO (2). Note in particular that Y SM t in Eq. (19) is always replaced by √ 2m t /v, while Y t is set to the Yukawa coupling of the model specified in VHATNNLO (2). All other couplings of the Higgs are taken into account according to the model under consideration.
Except for L 4 , which differs by the sign, our notation follows Ref. [84]. Fig. 14 depicts selected Feynman diagrams arising from the inclusion of higher dimensional operators. The operator L 6 is a dimension-5 operator which can have a significant impact in extended Higgs sectors and is therefore added to the Lagrangian. In the 2HDM or the MSSM, it corresponds to an effective coupling among two gluons and the pseudoscalar A at LO, relevant for the contribution gg → A → Zh. This involves the g Ah Z coupling, see Eq.
(3), as well as finite width effects for the pseudoscalar A, see the discussion in Sect. 5.1.
The coefficients of the operators in L i are specified through entries i in Block DIM in the vh@nnlo input file. We list them in Table 20. We note that by default all coefficients are understood to be constants. However, by specifying an input scale through DIM (11), the running of c tG to the renormalization scale following Ref. [85] is included. The strong coupling constant g s = √ 4πα s is evaluated at the renormalization scale µ R . In contrast to other input blocks which parametrize new physics, vh@nnlo inserts the default values given in Table 20; i.e. any unspecified coefficient of Eq. (19) is set to zero. The range of values for the coefficients is unrestricted in principle, but it is clear that very large values are potentially non-perturbative, and thus the results of vh@nnlo will be unreliable.
The effect of the dimension-6 operators is two-fold in general: on the one hand, they lead to additional vertices which do not occur in the model under consideration; on the other hand, they may modify the coupling constants of the existing vertices. vh@nnlo takes both of these effects into account. Some of these operators have been already implemented in Ref. [36] with a different normalization. To compare our results to Tab. 7 of Ref. [36] one has to transform the effective couplings as follows: where Λ is the energy scale of new physics. On the right-hand side of Eq. (20), the couplings are normalized and named as in Ref. [36].
The total cross section can now be written as where σ SM is the total cross section of the Standard Model, σ (1) contains the interference of the SM amplitudes which involves one modified vertex, and σ (2) contains the interference of amplitudes of different operators, where at most one modified vertex is inserted, or the square of an amplitude with one modified vertex. Again, the amplitudes have been calculated with the help of FormCalc.
In Table 21 we show the results for the total cross section of gluon initiated Zh production at √ s = 8 and 13 TeV by taking into account L 1,2,4,5 , with the couplings of Eq. (21), and c t = 0.1. The effect of L 3 is not shown since it leads to the same results as L 2 . We only enable one dimension-6 operator at a time. Furthermore, we fix the scales to µ R = µ F = 125 GeV, use the LO MSTW2008 [86] parton distribution functions, and set m b = 0. We observe full agreement with Ref. [36] within the numerical accuracy. As is wellknown, the largest effects of such operators typically occur at high transverse momenta or high invariant masses, where the validity of the effective field theory description becomes questionable. In the case of Zh production, however, the operators under consideration only affect the gg → ZH component of the cross section, which is known to affect mostly the kinematic region at and slightly above the top quark threshold. In order to investigate the impact of the dimension-6 operators in more detail, we employ the ratio of the full This quantity was shown in Ref. [87] to be particularly suited for a data-driven extraction of the gluon-initiated component σ gg→ZH from experiment. Using the example of L 1 and L 4 , Fig. 15 shows that large effects of the dimension-6 operators can already be observed in the few-hundred-GeV region of both p T and M ZH , i.e., well in the validity range of the effective theory. This observation adds to the virtues of the ZH process in the search of New Physics.

Conclusions
We presented version 2 of the code vh@nnlo, which allows to study various new-physics aspects in Higgs Strahlung. In detail we described the general structure of the code and its control through an input file containing SLHA-inspired blocks. We implemented the Higgs sectors of a 2HDM, which-in type II-also allows for the calculation of MSSM cross sections for Higgs Strahlung. For this purpose, we added the relevant squark amplitudes to the gluon-initiated contribution of Higgs Strahlung. In addition, the bottom-quark initiated contribution is given at leading-order in perturbation theory and resonances of (pseudo)scalars are included in both gluon-and bottom-quark initiated contributions. We demonstrated their relevance for light Higgs production in the MSSM. The shape of their interferences with the non-resonant Feynman diagrams can be studied, but a thorough analysis is left for future work. CP mixing among the three 2HDM scalars can be taken into account. Vector-like quarks can be studied in the SM, which-depending on the representation and the mixing angles-can have a significant effect on the inclusive cross section and the kinematical distributions. vh@nnlo can also provide the transverse momentum distribution of the final state Higgs boson, and invariant mass distributions of the two-particle final state. For the DY-like component, this requires a link to MCFM. We showed threshold effects of vector-like quarks and squarks in both distributions. Lastly, beyond such concrete model implementations, vh@nnlo now includes higher-dimensional operators. We compared our implementation to the literature and found agreement. It turns out that their effect is particularly pronounced directly above the top-quark threshold, which is well below the cut-off which is typically assumed for the underlying effective field theory. Apart from the already mentioned internal link to MCFM, the new version of vh@nnlo can also be linked to FeynHiggs and 2HDMC for the calculation of Higgs masses and mixings. Other dependencies on external codes are described throughout the manuscript.
A Installation and execution of vh@nnlo vh@nnlo can be downloaded from Ref.
[1]. After unpacking it, the user has to run ./configure, which determines local Fortran and C++ compilers and library dependencies. Thereafter, the user should open the Makefile of vh@nnlo and adjust the paths to the various needed external codes, see Appendix B.1. Then make will result in an executable, which is placed in the newly generated /bin folder. This executable is run by typing ./x.vhnnlo input.in output.out, where input.in has to be an input file, which contains the various blocks described throughout the manuscript. Example input files for the various new-physics models and for distributions can be found in the folder /example. If the name of the output file output.out is not specified, it is named out.vh. In case vh@nnlo is recompiled we recommend to type make recompile to make sure that a few dependencies are newly set up. If this does not result in the desired behavior, a previous make (some)clean can be helpful. If vh@nnlo should be linked to 2HDMC or FeynHiggs, the Makefile has to be processed with an additional predef=2HDMC or FH. If vh@nnlo was compiled beforehand, it has to be accompanied by recompile, i.e. make recompile predef=2HDMC or FH. Note that adding GGZH=NO allows to compile without the gluon-induced component. In that case, also the link to LoopTools is not required.

B Links to external codes
First, vh@nnlo has to be linked to LHAPDF for the usage of modern PDF sets, the CUBA package for numerical integration, LoopTools for the calculation of loop integrals and potentially MCFM for the calculation of kinematic distributions. As we described in Sect. 5.1 and Sect. 5.2 vh@nnlo can be linked to external codes for the calculation of Higgs boson masses and decay widths. Those are 2HDMC for the 2HDM and FeynHiggs for the MSSM. We subsequently describe the various needed and potential links to external codes.
B.1 LHAPDF, CUBA and LoopTools vh@nnlo needs to be linked to LHAPDF [53], which allows to make use of up-to-date PDF sets. We recommend to use LHAPDF version 6.2 or higher. Also older versions ≥ 6 work, but might produce an unpleasant flow of LHAPDF screen messages. After downloading [88] and installing LHAPDF into a local directory, the user has to specify the paths to the LHAPDF installation folder in the vh@nnlo Makefile variables LHAPDFDIR, LHAPDFINCDIR and LHAPDFLIBDIR, where the latter two folders have to include relevant header files and contain the LHAPDF library, respectively. If LHAPDF is installed in a global directory like /usr/share, etc., the specification of the paths might be obsolete. Note that the PDF sets, which were specified in the input file of vh@nnlo, need to be installed, i.  For numerical integration vh@nnlo makes use of the CUBA package [89] 6 which needs to be installed separately. The paths to the CUBA installation directory is then specified in the variable CUBADIR in the Makefile of vh@nnlo. For the numerical integration, vh@nnlo uses CUBA's implementation of VEGAS algorithm. The central integration parameters can be controlled through Block VEGAS in the input file of vh@nnlo, see Table 22. These are the number of evaluations the numerical integration is starting with, the increase in the number of evaluations, and the minimal and maximal number of evaluations. They are specified in VEGAS(10+10·i), VEGAS(11+10·i), VEGAS(12+10·i) and VEGAS(13+10·i), respectively, where i ∈ {0, 1, 2} labels the DY-like component, the gluon-and the bottomquark initiated contribution, respectively.
Finally, for the evaluation of loop integrals, vh@nnlo requires a link to LoopTools [16,91], which can be obtained from Ref. [92]. We recommend to use LoopTools version 2.13 or higher. 7 The link to the LoopTools installation directory and the library is specified in the variables LTDIR and LTLIBDIR of the vh@nnlo Makefile, respectively. If vh@nnlo is compiled with GGHZ=NO, i.e. the gluon-induced component is not included, the link to LoopTools is not needed.

B.2 MCFM
MCFM can be downloaded from its webpage [93]. It has to be compiled and installed separately. Since the implementation of vh@nnlo is currently not compatible with OpenMP, the MCFM library should be built and linked without OpenMP, if a link to vh@nnlo should be established. vh@nnlo can be linked to MCFM 8 by adjusting the path specified in the variable MCFMDIR and setting the corresponding flag MCFM=YES in the Makefile of vh@nnlo.
Since MCFM needs some additional files in the working directory, one has to call make MCFM to copy them into the vh@nnlo directory. The interface to MCFM does not rely on input files, thus all parameters needed by MCFM are set internally, except the input given in the input block DISTRIB. This includes the start, end and bin width of the histogram, cuts on the invariant mass of M VH and the integration points for VEGAS, see Table 10. All SM parameters, the factorization and renormalization scale and the parton distribution function are set to the same values as given in the blocks SMPARAMS, SCALES and PDFSPEC, respectively. To produce results comparable to vh@nnlo we disable the decays of H, Z, and W . This disables automatically any jets so that all input regarding jets is filled with dummy variables. The remaining input values for MCFM, which are automatically set by vh@nnlo, are summarized in Table 23. Remove decay width removebr .true.
Disable decay Another option is to use vh@nnlo through an additional SLHA-style spectrum file. Its name is provided in the vh@nnlo input file through SPECTRUMFILE(1). The additional SLHA-style spectrum file needs to contain the following blocks: • Block MASS, for the on-shell masses of the gluino, h, H, A boson, stop and sbottom quarks, • Block AU and Block AD, for the soft-breaking trilinear stop-Higgs and sbottom-Higgs coupling, respectively, • Block STOPMIX and Block SBOTMIX to evaluate the stop and sbottom mixing angle, • Block ALPHA, for the scalar mixing angle α, • Block HMIX to get the values of the µ parameter and tan β.
Note that these blocks are part of the standard output of typical SUSY spectrum generators such as SoftSusy [95], SPheno [96,97], or FlexibleSUSY [98][99][100]. Combining the latter with Himalaya [101] allows to run vh@nnlo in the MSSM by taking into account threeloop corrections to the light Higgs boson mass [102,103]. For now, vh@nnlo does not read in the Higgs boson decay width from the spectrum file, since not all external codes provide such numbers. Thus, the user needs to fill the entries MASS(250/350/360) in the vh@nnlo input file with the Higgs boson decay widths of the light CP even, the heavy CP even Higgs boson and the pseudoscalar, respectively. We refer to the example vh@nnlo input file together with a FlexibleSUSY spectrum file in folder /example.