NLO electroweak corrections in extended Higgs sectors with RECOLA2

We present the computer code RECOLA2 along with the first NLO electroweak corrections to Higgs production in vector-boson fusion and updated results for Higgs strahlung in the Two-Higgs-Doublet Model and Higgs-Singlet extension of the Standard Model. A fully automated procedure for the generation of tree-level and one-loop matrix elements in general models, including renormalization, is presented. We discuss the application of the Background-Field Method to the extended models. Numerical results for NLO electroweak cross sections are presented for different renormalization schemes in the Two-Higgs-Doublet Model and the Higgs-Singlet extension of the Standard Model. Finally, we present distributions for the production of a heavy Higgs boson.


Introduction
Since the discovery of a Higgs boson at the Large Hadron Collider (LHC) [1,2] the community is moving forward focusing on precision. Precision is the key to probe the Standard Model (SM) and Beyond Standard Model (BSM) physics and potentially allows, together with automation, to disprove the SM or even to single out new models. State of the art predictions involve typically two-loop and occasionally three-loop QCD and one-loop electroweak (EW) corrections for many processes of interest at the LHC. As the aim is to cover all accessible processes at the LHC and future colliders, a lot of effort has gone into the full automation of one-loop amplitudes. With one-loop QCD amplitudes being available since a long time, more recently much effort has been spent on the automation of EW one-loop corrections, which are more important than ever in view of the recent progress in multi-loop QCD calculations. SM EW corrections are nowadays available in various approaches, e.g. OpenLoops [3], MadGraph5 aMC@NLO [4], GoSam [5,6], FeynArts/FormCalc [7,8], and in our fully recursive approach RECOLA [9,10]. For BSM physics precision is important, and especially EW corrections should not be underestimated as they can be comparable to QCD corrections in certain BSM scenarios.
The automation for one-loop BSM physics requires three ingredients: First, new models need to be defined, typically in form of a Lagrangian and followed by the computation of the Feynman rules. For this kind of task Feynrules [11] and SARAH [12] are established tools. Then, a systematic and yet flexible approach to the renormalization and computation of further ingredients is required to deal with generic models. Finally, the renormalized model file needs to be interfaced to a generic one-loop matrix-element generator. As for the automation of renormalization, there has been progress in the Feynrules/FeynArts approach [13]. In this paper we present an alternative and fully automated procedure to the renormalization and computation of amplitudes in general models, thus, combining the second and third step. Our approach makes use of tree-level Universal FeynRules Output (UFO) model files [14] and results in renormalized one-loop model files for RECOLA2, a generalized version of RECOLA, allowing for the computation of any process in the underlying theory at the one-loop level, with limitations only due to available memory or CPU workload.
As an application of the system, we focus on two BSM Higgs-production processes at the LHC, namely Higgs production in association with a vector boson, usually referred to as Higgs strahlung, and Higgs production in association with two jets, known as vector-boson fusion (VBF), in the Two-Higgs-Doublet Model (2HDM) and the Higgs-Singlet extension of the SM (HSESM). Those processes are particularly interesting for an extended Higgs sector, as they represent the next-to-most-dominant Higgs-production mechanisms at the LHC. There has been enormous progress in higher-order calculations to Higgs strahlung and VBF in the SM and BSM. For Higgs strahlung the QCD corrections are known up to NNLO for inclusive [15][16][17] and differential [18,19] cross sections. On-shell EW corrections were computed in Ref. [20] and followed by the off-shell calculation in Ref. [21]. Higgs strahlung has also been investigated in the 2HDM for QCD [22] and EW [23] corrections. NLO QCD corrections matched to parton shower have been presented in Ref. [24] in an effective field theory framework. For VBF, the first one-loop QCD corrections were obtained in a structure function approach [25] followed by the first two-loop prediction [26,27] in the same framework. As for differential results, the first one-loop QCD and EW corrections were calculated in Ref. [28] and Refs. [29,30], respectively. Since recently also the differential two-loop [31] and three-loop [32] QCD corrections are available. VBF has been interfaced to parton showers [33,34] and has been subject to studies for a 100 TeV collider [35]. In view of BSM, VBF has been studied in the MSSM [36]. Higgs strahlung and VBF are nowadays available in public codes, such as V2HV [37], MCFM [38], HAWK2.0 [39] and vh@nnlo [40]. This paper is organized as follows. In Section 2 the computer program RECOLA2 is presented as a systematic approach towards the automated generation of one-loop processes. RECOLA2 relies on one-loop renormalized model files which are automatically generated with the new tool REPT1L from nothing but Feynman rules. The computation steps are explained in different subsections, where we discuss the translation from UFO to RECOLA2 model files (Section 2.1), the counterterm expansion and renormalization, and the computation of rational terms of type R 2 (Section 2.2). In Section 3 we give details on the HAWK 2.0 interface with RECOLA2, which has been used for the phenomenology. In Section 4 we list our conventions for the 2HDM and the HSESM, focusing on the physical input parameters. In Section 5 we discuss the application of the Background-Field Method (BFM) in RECOLA2. We present the renormalization for extended Higgs sectors in the BFM and give details on the implementation in REPT1L. In Section 6 we fix the calculational setup and define the benchmark points, which were mainly taken from the Higgs cross section working group (HXSWG). For the numerical analysis we use different renormalization conditions for the mixing angles, which we introduce in Section 6.3. In Section 7 we present the numerical results, discussing total cross sections in view of different renormalization schemes and distributions for heavy Higgs-boson production. After the conclusions in Section 8, we illustrate in App. A how the colour flow is derived and provide additional information on the derivation of a minimal basis for off-shell currents in App. B. Finally, in App. C we discuss the application of on-shell renormalization schemes combined with different tadpole counterterm schemes focusing on the gauge dependence.
2 RECOLA2: RECOLA for general models RECOLA2 is a tree-level and one-loop matrix-element provider for general models involving scalars, fermions and vector particles. It is based on its predecessor RECOLA [9,10], which uses Dyson-Schwinger equations [41][42][43] to compute matrix elements in a fully numerical and recursive approach. The implementation at tree level follows the strategy developed in Ref. [44], supplemented by a special treatment of the colour algebra. The one-loop extension, inspired by Ref. [45], relies on the decomposition of one-loop amplitudes as linear combination of tensor integrals and tensor coefficients. The former are evaluated by means of the library COL-LIER [46], while the latter can be computed by making use of similar recursion relations as for tree amplitudes. The key point is the construction of the proper tensor structure of the coefficients at each step of the recursive procedure, which has been implemented in RECOLA relying on the fact that in the Standard Model in the 't Hooft-Feynman gauge the combination (vertex)×(propagator) is at most linear in the momenta. RECOLA2 circumvents these and other limitations of RECOLA. In the following we give an introduction to RECOLA2 and its capabilities, focusing on the generalization with respect to RECOLA and on the applications presented in Section 7.
The generalization of RECOLA has required to remove all SM-based pieces of code, replacing them with generic structures which are able to retrieve any necessary information from the model file. Furthermore, the process-generation algorithm makes use of recursive functions dealing with different cases on equal footing. This has produced a more compact code as no model-dependent information has been hard-coded. Finally, RECOLA2 just needs the Feynman rules to be provided by model files in a specific format to directly evaluate NLO amplitudes in the model under consideration by using similar recursion relations to those of the SM. As for RECOLA, the key ingredients are the so-called off-shell currents w i (P, {n}) = n P (2.1) defined as the sum of all Feynman graphs which generate the off-shell particle P combining n external particles. 1 The generic index i is related to the spin. For example, in the case of a vector field i is a Lorentz index or in the case of a fermionic field i is a spinor index. Other indices are suppressed and not relevant for the following discussions.
The off-shell currents (2.1) are build recursively according to the Berends-Giele recursion relations (BGR) [47] = which constitute a generalization of Eq. (2.2) of Ref. [9] for general models where elementary couplings with more than four fields are present. Note that in RECOLA [10] the terms with λ i , i > 4 are absent as only 2-, 3-, and 4-point interaction vertices are supported. Practically, each term on the right-hand side of the BGR equation (2.2) combines off-shell currents, referred to as incoming currents, and contributes to the construction of the current on the left-hand side, referred to as outgoing current. An outgoing off-shell current with n external particles is calculated using the vertices of the theory connecting incoming off-shell currents with less than n external particles, which, when combined, add up to n external particles. This can be realized for tri-linear, quadri-linear, quinti-linear, or even higher n-point vertices if present in the theory. The contribution to the outgoing current generated in each term of equation (2.2) can be formally seen as the result of the action of the BGR operator defined by 3) can be written in a model-independent way as a linear combination of Lorentz structures from which the couplings, colour structures and other relevant information that needs to be propagated from the left to the right is factorized. RECOLA2 is fully relying on the model file to provide those rules, in addition to recursive rules for the colour-flow and helicity-state propagation. One could argue, that not too many different operators are required, at least for the renormalizable theories, which could have been hard-coded. However, in view of different conventions, different gauges and non-renormalizable theories, we decided for a flexible system by moving this dependence to the model file. As now the model file provides the rules for computing off-shell currents, we can easily incorporate the BFM and R ξ -gauge for the SM and BSM models for NLO computations which is discussed in Section 5. In addition, RECOLA2 has been generalized to deal with arbitrary n-point vertices, 2 and, thus, can compute processes with elementary interactions between more than four fields. Dealing with higher n-point vertices required to improve, among other parts of the code, the generation of the tree graphs of the process. The generation of those graphs is a combinatorial problem which is practically solved in the binary representation as introduced in Ref. [48] (see also Ref. [44]). For elementary interactions involving an arbitrary number of fields the method requires to compute distinct ordered integer partitions of arbitrary size with no bitwise overlap between elements.
Further, RECOLA2 allows for arbitrary powers of momenta 3 in Feynman rules, which is crucial for EFTs and the R ξ -gauge at one-loop level. In order to implement this important generalization, we had to generalize the construction of the tensor structure of loop currents (i.e. of the coefficients of the tensor integrals), allowing the combination (vertex)×(propagator) to contain any power of momenta.
New theories may involve new fundamental couplings, and RECOLA2 can deal with an arbitrary number of them. 4 The computation of matrix elements is ordered according to powers of fundamental couplings, and RECOLA2 provides methods to automatically compute amplitudes and interferences for all possible orders of these couplings. For instance, this feature can be used to control the number of insertions of a higher dimensional operator in a given amplitude.
Finally, RECOLA2 comes with almost all features and optimizations as provided by RECOLA. It is designed to be backward compatible in the sense that a program which successfully runs with RECOLA can be linked to RECOLA2 and a SM (or SM BFM) model file and is guaranteed to run without any code adaptation. This is realized by a dedicated SM interface which has been developed on top of the general interface to model files. The most notable optimizations concern partial factorization in colour-flow representation, the use of helicity conservation and the identification of fermion loops for different fermions with equal masses.
2.1 RECOLA2 model-file generation RECOLA2 model files are generated with the tool REPT1L (Recola's rEnormalization Procedure Tool at 1 Loop) which is a multi-purpose tool for analytic computations at the oneloop order. REPT1L is written in Python 2.7 5 and depends on other tools, most notably RECOLA2 for the model-independent current generation, which is used in combination with FORM [49] to construct analytic vertex functions or S-matrix elements, and SymPy [50], which is a computer-algebra system (CAS) for Python.
REPT1L requires the Feynman rules in the UFO format [14] which can be derived via Feynrules [11] or SARAH [12]. As there has been progress for an automated renormalization in the Feynrules framework [13], we stress that we do not require any results for counterterms or rational terms. Those terms are automatically derived from the tree-level Feynman rules in a self-contained fashion as explained in Section 2.2.
The RECOLA2 model-file generation consists of two phases. In the first phase REPT1L loops over all vertices in the UFO model file, disassembling each into the vertex particles, Lorentz and colour structures, and couplings. The colour structure is transformed to the colour-flow basis possibly rearranging Lorentz structures and couplings. This is discussed in more detail in App. A. The resulting Lorentz structures are used to derive the BGR operators in a model-independent way. For every Feynman rule REPT1L tries to map the encountered Lorentz structure onto one of those operators. If a new structure cannot be mapped onto an existing operator a new operator is added. In an optional second pass, the existing base of operators is minimized (see App. B for more details).
In the second phase of the model-file generation the information is exported as Fortran95 code in form of a model-file library as depicted in Fig. 1. Particle configurations are linked to the individual contributions on the right-hand side of (2.2), which differ in the underlying BGR (2.3), colour flow, colour factors, couplings, coupling orders or other information, via a Fortran95 hash table, allowing for a flexible and efficient access. The actual BGR are computed and exported as Fortran95 subroutines in different forms. For the numerical evaluation tree and loop BGR are used to construct tree-level and one-loop amplitudes as it is done in RECOLA. The tree BGR are a special case of the loop BGR, with no loop-momentum dependence. As a new feature in RECOLA2, an analytic version of the BGR allows to generate amplitudes as FORM code. 6 In this way the analytic expressions for the amplitudes needed in the renormalization conditions are derived in the same framework as the loop amplitudes of the computed processes, ensuring that properly defined renormalization schemes automatically imply UV-finite results in numerical computations. In general, the UV finiteness of the theory can (and should) be verified numerically in RECOLA2 process by process by varying the scale µ UV related to the dimensional regularization of UV singularities [10]. This check also works in combination with MS subtraction schemes, even though in this case amplitudes have an intrinsic scale dependence. To this end, we separate the scale dependence originating from the MS subtraction from the one of regularization.
Finally, RECOLA2 requires particle information such as the mass, spin, and colour of particles. This information is directly obtained from UFO particle instances and is translated to Fortran95 code. These steps conclude the tree-level model-file generation. In the next section we discuss the counterterm generation and renormalization and the computation of rational terms of type R 2 .

Counterterm expansion, renormalization and computation of R 2 terms
REPT1L supports an automated renormalization of model files following the standard procedure (see e.g. Ref. [51]). Here we give a short summary of all the steps, followed by details on the counterterm expansion, the renormalization conditions, and the computation of rational terms of type R 2 .
The starting point is a tree-level UFO model file. In the first step an independent set of parameters is identified, followed by a counterterm expansion. The RECOLA2 model file is derived, enabling the formal counterterm expansion in REPT1L and leaving the values for counterterm parameters unspecified. Renormalization conditions are used to fix the counterterm parameters. REPT1L allows to renormalize counterterm parameters in various schemes, and

Counterterm expansion
In the default setup, REPT1L defines the counterterm expansion rules of the masses M V , M S , m f , associated to scalars (S), vector bosons (V ) and fermions (f ), of the not necessarily physical bosonic (φ) and fermionic fields (ψ), and of a set of external couplings g k , according to 7 with δZ jl being, in general, a non-diagonal matrix and L, R denoting the left-and right-handed components of fermionic fields, which, by default, are assumed to be diagonal. REPT1L automatically deals with counterterm dependencies if the parameters, being assigned a counterterm expansion, are declared as external parameters in the UFO format. Here, an external parameter is an independent parameter, whereas internal parameters depend on external ones and their counterterm expansion can be determined by the chain rule. Once all parameters have a counterterm expansion, the most efficient way to generate counterterm vertices of the theory is through an expansion of the bare vertices via (2.4). It is possible to add counterterm vertices by hand, or, as a third alternative, to induce counterterm vertices from bare ones, which are not included in the model, via counterterm expansion rules. The latter is used to handle 2-point counterterms and counterterms originating from the gauge-fixing function since both of these types have no corresponding tree-level Feynman rules.

Renormalization conditions
A standard set of renormalization conditions is implemented in Python as conditions, rather than solutions to conditions, which are solved upon request. As an advantage of solving equations, the form of vertex functions or conventions can change without breaking the system. REPT1L supports on-shell, MS, and momentum-subtraction conditions for general (mixing-)two-point functions. MS subtraction is implemented generically for n-point functions. We assume standard renormalization of the physical fields and masses from the complex poles of Dyson-resummed propagators and their residues, while we allow for several choices of renormalization conditions for the gauge-fixing function and for unphysical fields. In addition, we provide standard renormalization conditions for the SM couplings, e.g. the definition of α in the Thomson limit (TL) and in the G F scheme, which are implemented via self-energies 8 , and the N f -flavour scheme for α s in QCD 9 , which is implemented as a combined MS/momentum subtraction on vertex functions. All conditions are implemented in a model-independent way. Instead of the standard set of renormalization conditions already implemented, REPT1L can also handle alternative conditions properly set by the user.
Setting up renormalization conditions requires a RECOLA2 model file including counterterms. The derivation of model files is done as discussed in the previous section with enabled vertex counterterm expansion (see Fig. 1) and leaving the counterterm parameter unspecified. The renormalization conditions are derived analytically as FORM code. REPT1L uses RECOLA2 to generate the skeletons for processes. The result is written to a FORM file and evaluated, yielding vertex functions which are parsed to Python and processed with SymPy solving the conditions for the counterterm parameters. The procedure is visualized in Fig. 2. Multiple schemes for the very same counterterm parameters can be implemented by imposing different renormalization conditions. All schemes are exported to the RECOLA2 model file and, for a given parameter, a specific scheme can be selected before the process generation phase. For instance, this system can be used to allow the user to choose between different QCD and EW renormalization schemes within the same model file. The same system is used for dealing with singularities from light fermions. In general, particles can be tagged as light particles, which, when a particle is subject to on-shell renormalization, makes REPT1L to regularize the associated diagonal two-point function in three different setups, namely dimensional regularization, mass regularization, and keeping the full mass dependence. In a RECOLA2 session a suited regularization scheme for light particles is set automatically, depending on the choice of the mass value, unless the regularization for a particle is explicitly required in a specific scheme. In the case of unstable particles, i.e. massive particles with finite widths, REPT1L applies, by default, the Complex-Mass Scheme (CMS) as discussed in more detail in Section 5.2.

Computation of R 2 terms
The computation of R 2 uses the methods developed in Refs. [52][53][54] and follows the same computation flow as solving renormalization conditions which is depicted in Fig. 2. For renormalizable theories all existing R 2 terms can be computed. To this end, REPT1L can generate the skeletons at NLO for all vertex functions in the theory which are potentially UV divergent by power counting. FORM is used to construct each vertex function, replace tensor integrals by their pole parts and take the limit D → 4. The finite parts are identified as Feynman rules associated to the original vertices, which are precisely the R 2 terms. These steps are done in Python with the help of SymPy.
The computation of tensor coefficients is done in conventional dimensional regularization. Different regularizations will be supported in the future by exchanging the responsible FORMprocedure files. In view of EFTs, the power counting can be disabled, and specific vertex functions can be selected. Further, the R 2 extraction rules [52][53][54] have been extended to higher n-point functions and higher rank. 10

HAWK 2.0 interface to RECOLA2
In this section we describe the interface between HAWK 2.0 and RECOLA2 which allows for an automated computation of NLO EW and QCD corrections to observables in associated Higgs production with a vector boson or two jets. We start with the LO partonic channels and virtual corrections and conclude with the computation of the real corrections. The implementation has been realized in a model-independent way, allowing in the future, apart from the two presented BSM models, for predictions in alternative models.

Process definitions at LO and NLO with RECOLA2
In the case of associated Higgs production with a vector boson, also known as Higgs strahlung, we consider processes with an intermediate vector boson decaying leptonically as pp → HV → Hl + l − /Hl ± ν/Hνν.  Whenever possible, we optimize computations involving different quark generations. For instance, in (3.2) the processes involving the second generation are not computed explicitly, but the results for the first generation are employed instead. For the first generation of quarks the RECOLA2 library is used to generate the processes at tree and one-loop level. The second process class under consideration is Higgs production in association with two hard jets pp → Hjj, (3.3) also known as VBF. There are plenty of partonic channels and, again, we exploit optimizations with respect to the different quark generations. For the LO, virtual NLO EW, virtual NLO QCD, real emission EW, and real emission QCD contributions RECOLA2 generates 32 partonic channels each, with the real kinematic channels corresponding to the Born kinematic ones, with an additional gluon or photon. For the gluon-and photon-induced channels RECOLA2 generates 20 channels each. At the stage of the process definition the Higgs boson entering in (3.1) or (3.3) can be chosen freely 11 as long as it is supported by the RECOLA2 model file currently in use. For instance, in the case of the 2HDM the Higgs flavour can be set to H l , H h or H a (see Section 4), which is done in the HAWK 2.0 input file. In HAWK 2.0 the relevant parameters for process generation and computation are set by input files. This information is forwarded to RECOLA2, allowing to choose specific contributions. The selection works for individual corrections such as QCD or EW either virtual or real. For the results presented in this work we selected the pure electroweak corrections, including photon-induced corrections.

Infrared divergences
RECOLA2 provides the amplitudes for the partonic processes under consideration as well as the colour-correlated squared matrix elements needed for the Catani-Seymour dipole subtraction. In order to deal with IR singularities, an IR subtraction scheme needs to be employed. We adhere to the Catani-Seymour dipole subtraction [55] which is used in HAWK 2.0 and employ mass regularization for soft and collinear divergences, i.e. a small photon mass and small fermion masses are used wherever needed. From the point of view of the interface, dealing with EW dipoles is a matter of replacing certain Born amplitudes with the ones computed by RECOLA2. As for the QCD dipoles one needs in general colour-correlated matrix elements. For processes with only two partons, as it is the case for Higgs strahlung, the colour correlation is diagonal owing to colour conservation (see Eq. A1 in Ref. [55]) and again no colour-correlated matrix elements are required. For VBF we compute the colour-correlated matrix elements directly with RECOLA2, and use colour conservation to minimize the number of required computations. The dipoles are used as implemented in HAWK 2.0 and are not part of RECOLA2. For the QCD dipoles consider Refs. [55,56] and for EW dipoles see Refs. [57,58].

2HDM and HSESM model description
In this section, we sketch the definition of the scalar potential of the 2HDM and the HSESM. In both cases we restrict ourselves to a CP-conserving Z 2 -symmetric scalar potential, which in the case of the 2HDM is allowed to be softly broken. For a comprehensive introduction to the 2HDM we refer to Refs. [59,60] and for the HSESM to the original literature [61][62][63] and applications to LHC phenomenology in Refs. [64][65][66][67]. For the kinetic terms we refer to the conventions used in Ref. [64].

Fields and potential definition
Both models are simple extensions of the SM, only affecting the form and fields entering the scalar potential and for the 2HDM also the Yukawa interactions. In the case of the 2HDM we have two Higgs doublets, generically denoted as Φ i with i = 1, 2 and defined component-wise by with v i denoting the vevs. Under the constraint of CP conservation plus the Z 2 symmetry with five real couplings λ 1 . . . λ 5 , two real mass parameters m 2 1 and m 2 2 , and the soft Z 2 -breaking parameter m 2 12 . The HSESM scalar potential involves one Higgs doublet Φ and a singlet field S defined as Under the same constraints, the most general, renormalizable potential reads with all parameters being real.

Parameters in the physical basis
Both potentials are subject to spontaneous symmetry breaking which requires a rotation of fields to the mass eigenstates in order to identify the physical degrees of freedom. For the 2HDM there are five physical Higgs bosons H l , H h , H a , H ± and in the HSESM there are two neutral Higgs bosons H l and H h , intentionally identified with the same symbols as in the 2HDM. Besides the physical Higgs bosons, there are the three would-be Goldstone bosons G 0 and G ± in the 't Hooft-Feynman gauge. The mass eigenstates for the neutral Higgs-boson fields are obtained in both models by the transformation and α being fixed such that the mass matrix , with the potential V being either (4.2) or (4.4). The solution to (4.6) for symmetric 2 × 2 matrices is generically given by (see Ref. [59]) In the 2HDM there are additional mixings between charged and pseudo-scalar bosons and Goldstone bosons, which are diagonalized as follows The angle β is related to the vevs according to The Higgs sector is minimally coupled to the gauge bosons.
Collecting quadratic terms and identifying the masses one obtains the well-known tree-level relations where g and g denote the weak isospin and hypercharge gauge couplings, and M W and M Z the W-and Z-boson masses, respectively. For the 2HDM we identify v = v 2 1 + v 2 2 . Finally one employs the minimum conditions for the scalar potential which, in both models, read cos β sin β , (4.11) and the two mixing angles as c αβ := cos(α − β) (s αβ := sin(α − β)) and t β := tan(β), which is a natural choice for studying (almost) aligned scenarios. For the HSESM we use the neutral Higgsboson masses M H l (light Higgs boson), M H h (heavy Higgs boson) and the angles s α := sin(α) and t β := tan(β). To summarize, we transform the parameters from the generic basis to the physical one by choosing the following parameters as external ones 2HDM: where we have indicated that the vev is traded for gauge couplings and masses according to (4.9).

Yukawa interactions
The fermionic sector in the HSESM is the same as in the SM, whereas the 2HDM allows for a richer structure. In the general case of the 2HDM, fermions can couple to both Φ 1 and Φ 2 , leading to flavour-changing neutral currents (FCNC) already at tree level. Since FCNC processes are extremely rare in nature they highly constrain BSM models. In order to prevent tree-level FCNC, one imposes the Z 2 symmetry whereΦ 2 is the charge-conjugated Higgs doublet of Φ 2 . Neglecting flavour mixing, the coefficients are directly expressed by the fermion masses m d , m u and m l , and the mixing angle β, (4.14) Again, the vev v has been substituted using Eq. (4.9).
Type II: This is the MSSM-like scenario obtained by requiring odd parity for down-type quarks and leptons: d R → −d R , l R → −l R and even parity for up-type quarks. It follows that the down-type quarks and leptons couple to Φ 1 , while up-type quarks couple to Φ 2 . The corresponding Yukawa Lagrangian reads Neglecting flavour mixing, the coefficients are expressed by the fermion masses m d , m u and m l , and the mixing angle β, Type Y: This type, also referred to as lepton-specific 2HDM, is obtained by requiring odd parity only for leptons: l R → −l R .
Type X: This type, also referred to as flipped 2HDM, is obtained by requiring odd parity only for down-type quarks: In the analysis of this paper we focus on Type II, which is equivalent to Type I for massless leptons and quarks, except for the top quark. We remark that exactly one RECOLA2 model file can handle all Yukawa types, and switching between different Yukawa types is done by a simple function call.

Background-Field Method for extended Higgs sectors
The BFM is a powerful formulation for gauge theories which renders analytic calculations easier due to a simple structure of the Feynman rules and additional symmetry relations. The method was originally derived by DeWitt in Refs. [70,71] 12 and has since then been used in many applications. The additional symmetry relations emerge for gauge theories in combination with a suited gauge-fixing term and encode the invariance of the theory under so-called Background-Field gauge invariance. This property is particularly useful for the calculations of β functions [73] in higher orders and is also of interest in beyond flat space-time quantum field theory. The BFM can be used to calculate S-matrix elements, as constructed in Ref. [74], which, despite having to deal with many more Feynman rules, is in our implementation as efficient as the conventional formalism. Further, the BFM, which can be viewed as a different choice of gauge, allows for an alternative way of computing S-matrix elements and, thus, provides a powerful check of the consistency of the REPT1L/RECOLA2 tool chain. This is particularly useful for the validation of R 2 terms where mistakes are difficult to spot. In addition, we checked a few Background-Field Ward identities. We stress that the BFM can be used as a complementary method in RECOLA2 besides the usual formulation. Even though the use of the BFM in practical calculations is steered in precisely the same way as for model files in the conventional formulation, the internal machinery is different. In particular, the derivation of the Feynman rules and renormalization procedure requires special attention which is discussed in the following.

BFM action for extended Higgs sector
The results presented here are a simple generalization of Ref. [75], which deals with the BFM applied to the SM at one-loop order. The BFM splits fields in background and quantum fields and combines the new action with a special choice for the gauge-fixing function resulting in a manifest background-field gauge invariance for the effective action at the quantum level. This splitting separates the classical solutions of the field equations, represented by background fields, from the quantum excitation modes, represented by quantum fields. The Feynman rules are derived as usual, treating background and quantum fields on equal footing, which we have done with Feynrules. In principle, the splitting can be done for every field in the theory, however, as we are only interested in a background-field gauge-invariant action, it is sufficient to shift fields which enter the gauge-fixing function. Thus, we perform where W a (Ŵ a ) and B (B) are the SM quantum (background) gauge fields in the gauge eigenbasis with a = 1, 2, 3. The index i runs over all Higgs doublets Φ i in the theory under consideration, and S is a singlet field, absent in the 2HDM or SM. While the singlet field S does not appear explicitly in the gauge-fixing function [see (5.3)], the inclusion of S in the splitting (5.1) is necessary due to the mixing with the neutral component of a Higgs doublet. The components for the background-and quantum-field doublets are defined as 12 See Ref. [72] for a pedagogical introduction.
By convention, we keep the original vev of the Higgs doublet in the Higgs background-field doublet. The quantum gauge-fixing term has the traditional form. In the gauge eigenbasis it reads with generalized gauge-fixing functions and i running over all Higgs doublets. The covariant derivativeD is similar to the usual one, but with a background-field gauge connection instead of a quantum-field one. For a field α j in the adjoint representation it acts in the following waŷ with abc being the structure constants of SU(2). The form (5.3), (5.4) is invariant under background-field gauge transformations, which can be shown using the techniques presented in Ref. [72], but suitably generalized in the presence of spontaneous symmetry breaking. The construction of the ghost term follows the standard BRST quantization procedure. Once the symmetry transformations are defined on the fields, a valid ghost Lagrangian, leading to a BRST invariant action, is given by The fields in the gauge eigenbasis are rotated to the physical basis in the following way The BRST transformations on the gauge eigenbasis, expressed in terms of physical fields via (5.7), read Note that in contrast to the conventional formalism, the covariant derivatives entering the BRST transformations use the shifted gauge fields (5.1). For the Higgs doublets the BRST transformation rules can be defined at the level of components as follows The transformations for δ B ρ i and δ B η i are fixed by taking the real and imaginary part of the BRST transformation of the lower doublet component, respectively. In this way, if the ghost term is formulated directly in the physical basis, as it is done in Ref. [75], the Lagrangian is guaranteed to be hermitian.

Renormalization in the BFM
The renormalization in the BFM is performed in the same fashion as in the conventional formulation, except that only background fields need to be renormalized. REPT1L can distinguish between both types of fields by checking the field-type attribute. A field can be assigned to be a background and/or quantum field. In the conventional formalism, all fields play both roles and can thus appear in tree and loop amplitudes. In the presence of pure quantum fields, as it is the case in the BFM, the only contributing Feynman rules to tree and one-loop amplitudes are the ones with exactly none or two quantum fields.
Since we aim at the computation of S-matrix elements, an on-shell renormalization of physical fields is suited. However, fixing the field renormalization constants via on-shell conditions breaks background-field gauge invariance and, as a consequence, some Ward identities are not fulfilled. The reason is that demanding background-field gauge invariance requires, in particular, a uniform renormalization of all covariant derivatives in the theory which is only possible if the field renormalization constants of gauge fields are not independent parameters but chosen accordingly [75]. Since the theory is governed by BRST invariance, the breaking of the background-field Ward identities does not pose a problem, especially not for the renormalizability of the theory and the gauge independence of observables. Yet, we do not break the QED background-field Ward identity, which relates the fermion-fermion-photon vertex to fermionic self-energies [75] and can be used to fix the photon field renormalization constant or the counterterm δZ e . Requiring (5.16) for renormalized vertex functions yields the well-known one-loop relation in the BFM and Σ 1PI,T being an unrenormalized transverse 1PI mixing or self-energy. Note that all terms, except for the W self-energy, originate from vertex and box corrections, in particular, the term Σ 1PI,T AZ has just been introduced to match the divergence structure. Equation (5.19) is valid for the conventional formulation in the 't Hooft-Feynman gauge, but not in the BFM since mixing and self-energies, or, in general, vertex functions differ by gauge-dependent terms in both formulations. Since the parameter ∆r connects physical quantities it is necessarily gauge independent, which implies that both formulations differ merely by a reshuffling of gauge-dependent terms between the self-energy and vertex parts. We have determined the difference in the vertex corrections between the BFM and conventional formulation in the 't Hooft-Feynman gauge, and, as expected, it cancels against the difference in the W self-energy. For a model-independent evaluation in the BFM, the result can be expressed in the same form as (5.19), but with a modified vertex correction 13 which is valid only in the 't Hooft-Feynman gauge in the BFM. Another subtlety concerns the renormalization within the CMS. REPT1L automatically renormalizes unstable particles in the CMS following the general prescription of Refs. [79][80][81]. The corresponding on-shell renormalization conditions require scalar integrals to be analytically continued to complex squared momenta. This can be avoided by using an expansion around real momentum arguments, 14 which gives rise to gauge-dependent terms of higher perturbative orders. Thus, comparing the BFM to the conventional formalism leads to somewhat different results for finite widths. The effect can be traced back to the difference of full self-energies in both formulations, e.g. the difference in the W self-energy is given by 13 Note that Σ 1PI,BFM,T AZ (0) is zero in the BFM due to a Ward identity.
14 The expansion breaks down for IR-singular contributions resulting from virtual gluons or photons. This can be corrected by including additional terms (see Ref. [80]) which is automatically handled in REPT1L.
with the conventions for scalar integrals as in Ref. [51]. The gauge dependence drops out in the mass renormalization constant, i.e. δµ 2 W BFM = δµ 2 W in the CMS, because the self-energy is evaluated on the complex pole, i.e. for p 2 = µ 2 W . However, performing an expansion of the selfenergy around the real mass M 2 W results in differences of the order of O α 3 . For a comparison of both formulations it is useful to modify the expanded (exp) mass counterterm to match the conventional formalism in the following way with B 0 being defined as the derivative of B 0 with respect to p 2 . Note that the difference is of order O(α 3 ) and phenomenologically irrelevant. The renormalization of the tadpoles in the BFM is performed analogously to the conventional formulation. From a theoretical point of view the renormalization of tadpoles is not necessary, and the theory is well-defined just by including tadpole graphs everywhere. However, in practical calculations it is desirable to avoid unnecessary computations of graphs with explicit tadpoles if their contribution can be included indirectly by other means, e.g. via a suited renormalization. The renormalization of the tadpoles has to be done with care because a naive treatment of the tadpole counterterms can lead to spurious dependencies on the gauge-fixing procedure which ultimately spoil the gauge independence of the one-loop part of S-matrix elements. From the point of view of applicability, automation and gauge independence, we strongly recommend the FJ Tadpole Scheme, 15 which has been automated for arbitrary theories [23]. In contrast to other schemes, the FJ Tadpole Scheme is purely based on the field reparametrization invariance of quantum field theory (see Ref. [23]), which can be shown to be equivalent to not renormalizing the tadpoles at all, but with the benefit of not having to compute graphs with explicit tadpoles. Under the general assumption that the theory under consideration is expressed in the physical basis without tree-level mixings and restricting to the one-loop case, the FJ Tadpole Scheme is equivalent to the field redefinitionĤ for every physical (background-)fieldĤ i that develops a vev and with δtĤ i being the associated tadpole counterterm. By fixing δtĤ i to the tadpole graphs TĤ i via explicit tadpoles are cancelled and only tadpole counterterms to 1PI graphs remain. REPT1L can automatically derive all tadpole counterterms in the FJ Tadpole Scheme. In the FJ Tadpole Scheme the value of each counterterm needs to be independent of δtĤ i which can be verified analytically. 16 Additional checks concerning the tadpole renormalization can be performed on a process-by-process basis by including the tadpole graphs explicitly instead of renormalizing them. Finally, we note that RECOLA2 is able to use any tadpole counterterm scheme, but only the FJ Tadpole Scheme is fully automated.
6 Setup and benchmark points 6

.1 Input parameters
For the numerical analysis in the two Higgs-boson production processes we use the following values for the SM input parameters [84]: For the 2HDM we present updated and new results for the benchmark points in Tables 1 and 3 as proposed by the HXSWG [85]. The corresponding Higgs self-couplings λ i are given for convenience in Tables 2 and 4. For the HSESM we compiled a list of benchmark points in Table 5 featuring different hierarchies and being compatible with the limits given in Refs. [65,66]. 17 The results include the SM-like and heavy Higgs-boson production for both models. The computations were carried out in the 't Hooft-Feynman gauge both in the conventional formalism and in the BFM. In case of the 2HDM the matrix elements have undergone additional tests. Most notably, we have compared results obtained with RECOLA2 for Higgs decays into four fermions, which is closely related to the considered processes, to an independent calculation [86] based on FeynArts/FormCalc [7,8] for various channels. We found agreement to more than 7 digits for 3348 out of 3500 phase-space points in the virtual amplitude, none differing by more than 5 digits. We compared off-shell two-point functions for all distinct external states, i.e. scalars, fermions, and vector bosons, against an independent approach in QGRAF [87] and QGS, which is an extension of GraphShot [88]. Against the same setup we compared Higgs decays into scalars, fermions and vector bosons on amplitude level. In addition, we verified (on-shell) Slavnov-Taylor identities for two-point functions (see Eq. (4.16) and the following in Ref. [86]).

Cut setup
For the analysis of Higgs strahlung we consider the case of two charged muons in the final state, pp → Hµ + µ − + X. The muons are not recombined with collinear photons, and are assumed to be perfectly isolated, treated as bare muons as described in Ref. [21]. We use the cuts given in Ref. [89], i.e. we demand the muons to • have transverse momentum p T,l > 20 GeV for l = µ + , µ − , • be central with rapidity |y l | < 2.4 for l = µ + , µ − , • have a pair invariant mass m µµ of 75 GeV < m µµ < 105 GeV.
Further, we select boosted events with a • transverse momentum p T,µµ > 160 GeV.
For VBF we employ the cuts as suggested by the HXSWG in Ref. [85], i.e. we require two hard jets j i , i = 1, 2, emerging from partons i with  Table 1: 2HDM benchmark points in the alignment limit, i.e. s αβ → −1, c αβ → 0, taken from Ref. [92]. The parameter M sb depends on the other parameters and is given for convenience.  Table 2: Higgs self-couplings for the 2HDM benchmark points in the alignment limit. We omit the imaginary parts appearing in the CMS.
We present the results for hadronic cross sections at the centre-of-mass energy of 13 TeV using the NLO PDF set NNPDF2.3 with QED corrections [91].

Mixing angles at one-loop order
The prime vertices of interest in the processes studied in Section 7 are the H l VV and H h VV vertices. Thus, the relevant one-loop corrections require to renormalize α and β in the 2HDM   Table 4: Higgs self-couplings for the 2HDM benchmark points outside the alignment limit. We omit the imaginary parts appearing in the CMS.  and α, but not β, in the HSESM. We present the counterterms for the mixing angles in an MS scheme and two different on-shell schemes in the following: MS: The mixing angles α, β are renormalized using MS subtraction [23] for the vertices H l → τ + τ − , H a → τ + τ − , respectively, with β only being renormalized in the 2HDM. This is equivalent to using the identities with the relation for δα being valid in the 2HDM and the HSESM and the one for δβ only in the 2HDM. The origin of these relations can be traced back to the renormalizability of models in a minimal (symmetric) renormalization scheme. See Ref. [86] for the derivation of these and other UV-pole-part identities. The tadpole counterterms in (6.2) are treated in the FJ Tadpole Scheme (see Apps. A and B in Ref. [23]) and using the renormalization condition (5.24) for tadpoles. Estimating the size of higher-order contributions via the usual scale variations has been improved via a partial resummation including the renormalization-group (RG) running of parameters. 18 For the 2HDM this requires to solve a coupled system of differential equations, 3) The functions f α , f β and f M sb can be directly read off the pole parts of the corresponding counterterms. The counterterm δM sb was fixed from the vertex H h → H + H − in the MS scheme. Note that δM sb does not enter the considered processes at fixed one-loop order.
For the HSESM we keep β fixed, assuming no scale dependence, resulting in a simple differential equation for α, The (coupled) system has been solved to run the parameters from the reference scale µ 0 to µ = µ 0 /2 and µ = 2µ 0 . The results are presented in Tables 6 and 7 for the benchmark  points of Tables 1, 3, and 5 being defined at the typical scale of the process, µ 0 = 2M H l if not stated otherwise. 19 The cross sections are evaluated at the scales µ 0 /2, µ 0 , 2µ 0 , using the running parameters of c αβ , t β , M sb (s α ) at the corresponding scale as input parameters in the 2HDM (HSESM). The three different predictions for σ EW NLO normalized to the leading-order cross section σ LO (µ 0 ) at the central scale µ 0 and scale-dependent relative EW corrections are defined as Note that the tree-level matrix elements only depend on the scale through the running of parameters, whereas the one-loop matrix elements have an explicit scale dependence. As a shorthand notation for the relative corrections in the MS scheme we use with u and d being the upper and lower edges of the scale variation (see e.g. Table 8).
BP  Table 6: Running values for t β , c αβ and M sb in the 2HDM at the scales µ 0 /2 and 2µ 0 . The benchmark points are defined at the central scale µ 0 in Tables 1 and 3. The results for the alignment-limit scenarios are in the upper part of the table whereas the non-alignment scenarios are in the lower part. For BP22A the running β reaches π/2 for a scale greater than µ 0 /2, thus, t β becomes singular. In this particular scenario the steep running is caused by the Higgs selfcoupling λ 3 ≈ 8 (see Table 4) and can be stabilized by reducing its value. The running becomes stable only for values smaller than λ 3 0.5. Changing λ 3 to 0.  Table 7: Running values for s α in the HSESM at the scales µ 0 /2 and 2µ 0 . The benchmark points are defined at the central scale µ 0 in Table 5.
p * : The renormalized mixing angles α and β are defined to diagonalize radiatively corrected mass matrices which implies a scale and momentum dependence for the mixing angles. The scale dependence can be eliminated by a special choice for the momentum p 2 = (p * ) 2 at which the mixing two-point functions, and thus the running mixing angles, are evaluated. The original idea goes back to Ref. [95] (see also Ref. [96]) and has been applied to the HSESM in Ref. [94] and the 2HDM in Ref. [83]. In our conventions, the counterterms are defined as Note that for δβ alternatively the mixing energy with the charged Higgs and Goldstone boson can be used. As the mixing energies are gauge-dependent an additional intrinsic prescription is required to fix the gauge-independent parts. We choose the BFM with quantum gauge parameter ξ Q = 1, which corresponds to the gauge-fixing functions (5.3), (5.4). We remark that this is equivalent [97,98] to the self-energy in the pinch technique [99,100] and allows to extract a well-defined gauge-parameter-independent contribution to self-energies or, in general, vertex functions and hence counterterms in this scheme.
BFM: As an on-shell alternative to the p * scheme, the authors of Ref. [83] propose to use the on-shell scalar mixing energies defined within the pinch technique which has also been investigated in Ref. [101]. In our framework, this corresponds to with the mixing energies evaluated in the BFM with quantum gauge parameter ξ Q = 1.
In Ref. [83] it is argued that the use of the FJ Tadpole Scheme is essential for the consistency of on-shell schemes in combination with (5.24). There are, however, other options. A different tadpole counterterm scheme, such as the one of Ref. [51], yields different values and pole parts for counterterms, e.g. δα and δβ absorb tadpoles and become gauge dependent. Yet, the absorbed tadpoles drop out in momentum subtraction schemes [23] and do not spoil the gauge independence of the S-matrix. In tadpole counterterm schemes other than the FJ Tadpole Scheme special care has to be devoted to the formulation of renormalization conditions as they are necessarily gauge dependent. This situation is encountered in standard SM and MSSM on-shell renormalization schemes, where certain tadpole contributions to self-energies are left out, rendering the counterterms gauge dependent, but the S-matrix remains gauge independent. When employing gauge-fixing prescriptions in renormalization conditions, tadpoles can be handled naively in a favoured tadpole counterterm scheme if the same gauge is used in the renormalization and in the matrix-element evaluation. This is illustrated in App. C using the example of δα in the p * scheme. There, we also discuss the general case with arbitrary gauge-fixing functions, which is less trivial and cannot be done in the naive way due to the mismatch of the gauge prescription and the actual gauge-parameter choice. From a practical point of view the latter is only relevant if one is interested in verifying the gauge independence of the S-matrix in tadpole counterterm schemes other than the FJ Tadpole Scheme. We note that the use of MS schemes for the mixing angles in combination with alternative tadpole counterterm schemes can be made gauge independent by including finite tadpole terms 20 which is equivalent to the use of the FJ Tadpole Scheme.
The results for total cross sections in the BFM renormalization scheme in Section 7 were not computed directly, but were obtained from the results in the p * scheme using the following formulas, depending on the model (2HDM or HSESM) and on the produced Higgs flavour (H l or H h ) as follows  Table 8: Relative NLO corrections δ EW to SM-like Higgs-boson production in Higgs strahlung pp → H l µ − µ + in alignment scenarios in the 2HDM. The results in the MS scheme are given at the central scale µ 0 = 2M h = 250.18 GeV with scale uncertainties estimated including the RG running of parameters as given by (6.7). Both on-shell schemes agree within the integration error, and only results in the p * scheme are given. The SM EW correction is δ EW = −12.4%.
HSESM H h : δ BFM EW = δ p * EW + 2 sα cα δα p * − δα BFM Note that the formulas can be applied uniquely to the observables under consideration as these rely on the mixing-angle dependencies of the respective leading-order couplings.

Numerical results for total cross sections
In Table 8 we present updated results for the production of a SM-like Higgs boson in Higgs strahlung in the 2HDM in alignment scenarios. Non-alignment scenarios are given in Table 9. The corresponding SM correction is −12.4%. In Table 10 we provide the corresponding results for heavy Higgs-boson production in non-alignment scenarios. For the HSESM all considered scenarios are non-aligned. The results for light Higgs-boson production are given in Table 11, and the ones for heavy Higgs-boson production in Table 12. Note that for the benchmark points BP5 and BP6 with inverted hierarchy the heavy Higgs boson is SM-like with M H h = 125.09 GeV. For the benchmark points in the 2HDM the light Higgs boson is always identified as the SM Higgs boson. Finally, in Table 13 results for SM-like and heavy Higgs-boson production in VBF are presented for the 2HDM. The HSESM predictions for VBF are given in Table 14. The corresponding SM correction for SM-like Higgs-boson production in VBF amounts to −5.5%.

Discussion of the numerical results
In the following, we compare cross sections in different renormalization schemes and models for Higgs-boson production in Higgs strahlung. For VBF the picture is similar and not discussed in detail. In particular, for the MS scheme we collect some observations concerning large corrections. An analysis of the exact origin of these contributions would go beyond the scope of this paper.

MS scheme
We start with the MS scheme and SM-like Higgs production in the alignment limit of the 2HDM in Table 8. In a fixed-order calculation no scale dependence appears in the MS scheme, because the relevant counterterms δZ H h H l , δα and δβ entering the vertices H l ZZ and H l WW are screened by the factor c αβ /s αβ = 0 in the alignment limit. For the same reason, the on-shell schemes agree with the MS scheme at the central value. Yet, with the running of parameters, a small scale dependence is visible. For BP21B the correction is unstable for smaller scales, signalling a potential problem with the benchmark point (in fact, this scenario is close to the non-perturbative limit, see Table 2.) or with the MS scheme. In non-alignment scenarios the MS results for the 2HDM in Table 9 are almost all unstable and suffer from large scale dependencies, 21 which are reflected in the running parameters c αβ and t β in Table 6. For heavy Higgs-boson production in   the 2HDM (Table 10) no predictions in the MS scheme are presented as these scale uncertainties are even more enhanced due to ratios s αβ /c αβ entering the predictions. The situation for the MS renormalization in the HSESM for light (Table 11) and heavy (Table 12) Higgs-boson production is clearly more stable for the considered benchmark scenarios (see Table 5 for the λ i values). This is reflected in a reasonable running of the parameter s α in Table 7, except for BP5 and arguably for BP6. Due to the smaller running, we obtain results in the expected ballpark, with no artificially large corrections, even for the heavy Higgs-boson production near the alignment limit, where potential problems coming from the mixing energy would be enhanced by uncancelled finite parts. In the HSESM large scale uncertainties are observed in the MS scheme for light Higgs-boson production in BP6 and in particular for almost degenerate neutral Higgs bosons in BP5. Further, one observes that the MS scheme leads to larger deviations from the SM corrections, which, however, do not come with large scale uncertainties for the well-behaved benchmark points.
In   neutral scalar mixing energy, which is not subtracted in the MS scheme. This particular effect is enhanced for degenerate neutral Higgs bosons owing to the denominator structure in (6.2) which is not cancelled against the one coming from the on-shell off-diagonal field renormalization constants. Besides the top-quark contributions it is possible to induce moderate contributions coming from the Higgs potential by tuning λ 3 . This requires, however, large M 2 H h − M 2 H l with not too small s α , and it is not straightforward to tune the parameters in order to exceed the top-quark contribution without getting close to the non-perturbative limit |λ i | ∼ 4π. In the 2HDM, the reason for the large corrections in the MS scheme is more difficult to grasp, especially because in view of our observables we have to deal with the renormalization of β which is known to cause difficulties in the MSSM [102]. 22 The problem with β can be traced back to large contributions in the tadpoles. For α, the largest contributions cannot be explained by tadpoles nor by the top-quark contribution in the neutral scalar mixing energy. Here, we observe that the large contributions to the neutral scalar mixing energy are mediated through the charged and pseudo-scalar Higgs boson, which, eventually, exceed all other contributions. Since these large contributions are only found in the off-diagonal LSZ-factors they remain uncancelled in the MS scheme.

On-shell schemes
For the considered on-shell schemes none of the observed problems of the MS scheme is encountered because the large contributions in the mixing energy and the tadpoles are subtracted via δα and δβ, i.e. all terms involving the poles 1/(M 2 H h −M 2 H l ) and 1/M 2 Ha cancel in S-matrix elements. Further, in view of the size of the corrections the on-shell methods perform much better in the sense that the SM-like Higgs-boson production processes (see Tables 8, 9 for the 2HDM and  Tables 11, 12 for the HSESM) yield corrections which are close to the SM correction. In heavy Higgs-boson production (see Table 10 for the 2HDM and Table 12 for the HSESM) the results for on-shell renormalization schemes remain stable even for aligned 23 or degenerate scenarios. The difference between the p * and BFM schemes is tiny. It seems to us that the schemes are too similar for their difference to provide a qualitative estimate of higher orders. The difference between these schemes just results from the momentum dependence of the neutral scalar mixing energy, which turns out to be small and starts at the order O(M 2 H l − M 2 H h ). Note also, that the large contributions in the neutral scalar mixing energy were observed to have almost no momentum dependence. For VBF the computation has only been carried out in the on-shell schemes. For the 2HDM (Table 13) and HSESM (Table 14) the SM-like scenarios almost coincide with the SM predictions.

Distributions
We present distributions for the transverse momentum p T,H h and rapidity y H h of heavy Higgs bosons in Higgs strahlung and VBF. In addition, we show distributions in the rapidity y µ − of the muon µ − in Higgs strahlung and in the rapidity y j 1 of the hardest jet j 1 in VBF. We selected a typical subset of all benchmark points, namely the benchmark points BP3B1, BP43 Fig. 3 and Fig. 4, respectively, the ones for y H h in Higgs strahlung and VBF in Fig. 5 and Fig. 6, and those for y µ − and y j 1 in Fig. 7 and Fig. 8. 24 In the upper plots we show the LO and NLO EW differential cross section. In the lower plots the relative EW corrections δ EW are depicted. In order to isolate the genuine effects of the underlying model from the kinematic ones, we have computed the pure SM corrections with the SM Higgs-boson mass set to the heavy Higgs-boson mass M H h denoted as "SM" in the lower panels. The corresponding SM total EW cross sections are listed in Table 15.
In the following we focus on shape-distortion effects relative to the SM results. Starting with the distributions in Higgs strahlung, we observe quite large effects in the p T,H h distribution in Fig. 3 for BP3B1 and BP43 in the 2HDM, small effects for BP3 in the HSESM, and no effect in BP45 in the 2HDM, which perfectly reproduces the SM result. The situation changes for the distributions in the rapidities y H h and y µ − in Figs. 5 and 7. Here, the largest deviations from the SM are observed for BP43, where the relative EW corrections to the y H h distribution in the 2HDM are flatter than in the SM. For the y µ − curve the opposite tendency is observed, i.e. the SM correction is flatter. For BP3B1, BP45, and BP3 shape distortions relative to the SM appear at large rapidities, which are less important due to low statistics in those regions. Switching to the distributions for VBF in Figs. 4, 6, and 8, we observe a stronger trend towards SM-like results. The largest differences are observed for BP43 in the p T,H h and y j 1 distributions. For BP3B1 the effects for the same distributions are smaller but significant. For BP3 the shape distortion in the p T,H h distribution for VBF is not larger than the one for Higgs strahlung. In             general in the considered benchmark points for the HSESM the effects in VBF, but also in Higgs strahlung, are tiny compared to the ones observed in the 2HDM. The reason for the rather mild effects in the HSESM is due to the similiarity of the HVV vertices to the SM ones. In particular, in the HSESM all couplings of the light and heavy Higgs boson to gauge bosons or fermions are SM-like, but modulated with (−s α ) and c α , respectively. In the relative corrections these factors drop out, and the only difference due to the presence of an additional light Higgs boson and modified Higgs-boson couplings is small in the benchmark points under consideration (all λ i 1). Remarkably, even the corrections to the HVV vertices involving Higgs self-couplings (and thus all corrections) scale as the corresponding tree level with either (−s α ) or c α , respectively, in the (anti-)alignment limit. Furthermore, all mixing effects between H l and H h vanish in this limit. For these reasons the corrections cannot become enhanced with respect to the LO unless tree-level perturbativity is violated. In fact, in the HSESM the one-loop corrected H h VV vertices are exactly zero in the alignment limit. Note that these arguments apply to the whole phase-space region, thus, no significant shape-distortion effects are expected for the processes under consideration in the HSESM.
The 2HDM, on the other hand, exhibits non-decoupling effects in the alignment limit c αβ → 0, where the underlying vertices for heavy Higgs-boson production become loop induced. The largest corrections in BP43 are due to the non-decoupling term in the top Yukawa coupling 25 proportional to s αβ . In this case the Yukawa coupling is of the same size as the corresponding SM one, but with a different sign and further enhanced with respect to the LO by a factor of 1/c αβ , leading to a non-SM-like bosonic-fermionic interplay. Furthermore, the corrections in the 2HDM are very sensitive to the presence of new particles, especially the pseudo-scalar Higgs boson in the case of BP43. In general, the contributions involving Higgs self-couplings can be large since non-decoupling terms remain in the alignment limit giving rise to enhanced corrections with respect to the LO.

Conclusion
We reported progress towards fully automated one-loop computations in BSM models. The presented code RECOLA2 allows one to compute QCD and EW corrections for extensions of the SM for arbitrary processes. RECOLA2 can produce NLO corrections in general models, which requires the model file for each BSM model built in a specific format containing the ordinary, counterterm and R 2 Feynman rules. The model-file generation and the renormalization of general quantum-field-theoretic models is performed with the new tool REPT1L in a fully automated way, relying on nothing but the Feynman rules of the model in the UFO format. Once the renormalization conditions for the model are established, REPT1L performs the renormalization, computes the R 2 rational terms and builds the one-loop renormalized model files in the RECOLA2 format. We introduced the Background-Field Method as a complementary method in RECOLA2, which is useful for practical calculations and serves as a powerful validation method. We described the renormalization procedure in the Background-Field Method which is handled in RECOLA2 on equal footing with the usual formulation.
In summary, we realized the following generalizations with respect to RECOLA: • We developed a true model-independent amplitude provider, featuring a dynamic process generation in memory without the need for intermediate compilation.
• A generic interface has been developed supporting all methods available in RECOLA, but generalized to fit in the model-file approach. This includes the computation of amplitudes and squared amplitudes, the selection of specific polarizations and resonances, and the computation of interferences with different powers in new fundamental couplings. Furthermore, we provide spin-and colour-correlated squared matrix elements required in the Catani-Seymour dipole formalism. The latter methods are restricted to singlet, triplet and octet states of SU(3).
• RECOLA2 is limited to scalars, Dirac fermions and vector bosons. In the near future we will allow for Majorana fermions.
• We support Feynman rules with a general polynomial momentum dependence and allow for elementary interactions between more than four fields. Due to internal optimizations the number of fields per elementary interaction is restricted to at most 8.
• We generalized RECOLA2 to support the BFM as a complementary method. Furthermore, the R ξ -gauge can be used for massive vector bosons or, alternatively, non-linear gauges can be implemented.
• With REPT1L we have formed the basis for a fully automated generation of renormalized model files for RECOLA2. We provide a simple framework for the implementation of custom renormalization conditions. Presently available model files for RECOLA2 include the Z 2 -symmetric Two-Higgs-Doublet Model with all types of Yukawa interactions and the Higgs-Singlet extension of the Standard Model as well as models files with anomalous triple vector-boson and Higgs-vector-boson couplings.
The considered simple models do by far not exhaust the range of applicability of RECOLA2 and REPT1L, and further models will be implemented in the future.
As an application of the new tools we present first results for NLO electroweak corrections to vector-boson fusion and updated results for Higgs strahlung in the Two-Higgs-Doublet Model and the Higgs-Singlet extension of the Standard Model. We compared Higgs-production cross sections for different renormalization schemes in both models. We analysed the scale dependence in an MS renormalization scheme for the mixing angles, which has been improved including the renormalization-group running of parameters. We found unnaturally large corrections and scale uncertainties at one-loop order for the MS scheme, while the considered on-shell schemes remain well-behaved. These enhanced contributions can be related to uncancelled finite parts in the MS scheme and should be investigated in more detail in the future, since a proper estimation of higher-order uncertainties, as it can be done based on scale variation in MS schemes, is highly desirable. For the on-shell schemes, our results for the electroweak corrections to SM-like Higgsboson production are almost not distinguishable from the corresponding SM corrections for all considered benchmark points. Finally, we presented distributions for the production of heavy Higgs bosons. Here, interesting shape-distortion effects for the electroweak corrections at the level of several percent are observed in the 2HDM.

Appendices A Colour-flow vertices
In RECOLA2 the colour flow is constructed recursively. For a given off-shell current the outgoing colour configuration is determined from the incoming ones and the possible colour flows associated to the interaction vertex. As the UFO format does not incorporate the colour flow, we need to translate between the two representations. We implemented a dynamical system for computing the colour flow from the generators and structure constants, rather than substituting for known results. In the conventions of Ref. [10] the colour flow associated to a given colour structure with λ a being the Gell-Mann matrices. The computation then consists of eliminating the structure constants and the generators by solving (A.2) for the structure constants and using the (Fierz) completeness relation for the generators as follows Performing all contractions yields a sum of Kronecker deltas which represent the individual colour flows. For instance, the quartic gluon vertex of the SM reads with L 1 , L 2 , L 3 being Lorentz structures which, for the following discussion, are left unspecified. Focusing on the colour structure δ i 1 j 2 δ i 2 j 3 δ i 3 j 4 δ i 4 j 1 , we obtain for the two relevant contributions k,a 1 ,a 2 ,a 3 ,a 4 (∆ a 1 ) i 1 j 1 (∆ a 2 ) i 2 j 2 (∆ a 3 ) i 3 j 3 (∆ a 4 ) i 4 j 4 f ka 1 a 2 f ka 3 a 4 = 1 2 −δ i 1 j 2 δ i 2 j 3 δ i 3 j 4 δ i 4 j 1 + . . . , k,a 1 ,a 2 ,a 3 ,a 4 (∆ a 1 ) i 1 j 1 (∆ a 2 ) i 2 j 2 (∆ a 3 ) i 3 j 3 (∆ a 4 ) i 4 j 4 f ka 1 a 4 f ka 2 a 3 = 1 2 +δ i 1 j 2 δ i 2 j 3 δ i 3 j 4 δ i 4 j 1 + . . . . (A.5) 26 The ∆a generators are related to the conventional ones Ta, as used e.g. in Feynrules, via (∆a) i j = √ 2(Ta) i j with Tr {TaT b } = δ ab /2 and [Ta, T b ] = if abc Tc. Note that the structure constantsf abc in Ref. [10] are related to the ones in this paper viaf abc = √ 2f abc .

Combining this result with (A.4), results in the contribution
Thus, diagonalizing the vertex in colour-flow basis requires, in general, to redefine Lorentz structures and couplings.

B Off-shell currents
For a given Lorentz structure and a definite colour-flow state the BGR is derived from the Feynman rules by selecting one of the particles as the outgoing one, multiplying with the corresponding propagator and the incoming currents of the other particles. Since the structure of currents depends on the outgoing particle, one needs to derive the BGR for all distinct outgoing particles. Consider for instance the QED vertex e + e − γ. REPT1L constructs three different recursion relations with w i ,w j , w µ being either incoming or outgoing off-shell currents, depending on whether they are on the right-or left-hand side of (B.1). For many Feynman rules, the underlying BGR are formally the same if the couplings or masses of the particles are not further specified. Assuming that the colour flow has been factorized as explained in App. A, all fermion-fermion-vector rules, e.g. Ze → e or γe → e, can be mapped onto the same structures realizing that γ µ ω + and γ µ ω − form a suitable basis, Df αβ c 1 γ µ ω + + c 2 γ µ ω − δα × w µ ×w δ , with D f αβ , Df αβ , D V µν denoting generic propagators for fermions, anti-fermions and vector bosons, respectively. REPT1L has the ability to derive a minimal basis, dynamically, i.e. depending on the operators of the theory, without relying on the Lorentz basis in the UFO format. This is done in two steps. In the first step, all distinct BGR in the underlying theory are registered. In the next step the BGR are merged recursively until a minimal basis is obtained. The size of the BGR can be controlled by a parameter for the maximal number of allowed distinct generic couplings, and it is possible to allow for vanishing couplings to improve the merging. If a merge yields a BGR size larger than allowed, the merging is not accepted. Finally, all vertices are mapped to the minimal basis.