External leg corrections as an origin of large logarithms

The appearance of large logarithmic corrections is a well-known phenomenon in the presence of widely separated mass scales. In this work, we point out the existence of large Sudakov-like logarithmic contributions related to external-leg corrections of heavy scalar particles which cannot be resummed straight-forwardly using renormalisation group equations. Based on a toy model, we discuss in detail how these corrections appear in theories containing at least one light and one heavy particle that couple to each other with a potentially large trilinear coupling. We show how the occurrence of the large logarithms is related to infrared singularities. In addition to a discussion at the one-loop level, we also explicitly derive the two-loop corrections containing the large logarithms. We point out in this context the importance of choosing an on-shell-like renormalisation scheme. As exemplary applications, we present results for the two-loop external-leg corrections for the decay of a gluino into a scalar top quark and a top quark in the Minimal Supersymmetric extension of the Standard Model as well as for a heavy Higgs boson decay into two tau leptons in the singlet-extended Two-Higgs-Doublet Model.

The appearance of large logarithmic corrections is a well-known phenomenon in the presence of widely separated mass scales. In this work, we point out the existence of large Sudakov-like logarithmic contributions related to externalleg corrections of heavy scalar particles which cannot be resummed straightforwardly using renormalisation group equations. Based on a toy model, we discuss in detail how these corrections appear in theories containing at least one light and one heavy particle that couple to each other with a potentially large trilinear coupling. We show how the occurrence of the large logarithms is related to infrared singularities. In addition to a discussion at the oneloop level, we also explicitly derive the two-loop corrections containing the large logarithms. We point out in this context the importance of choosing an on-shell-like renormalisation scheme. As exemplary applications, we present results for the two-loop external-leg corrections for the decay of a gluino into a scalar top quark and a top quark in the Minimal Supersymmetric extension of the Standard Model as well as for a heavy Higgs boson decay into two tau leptons in the singlet-extended Two-Higgs-Doublet Model.

Introduction
Until now, only one scalar particle without a known substructure has been found: the Higgs boson discovered at the Large Hadron Collider (LHC) in 2012 [1,2]. In the Standard Model (SM), the Higgs boson arises as part of an SU (2) L doublet alongside with the neutral and charged "would-be" Goldstone fields. For suitable parameter choices, the presence of the Higgs potential triggers the spontaneous breakdown of the SU (2) L × U (1) Y gauge symmetry. As a direct consequence, the "would-be" Goldstone fields become the longitudinal degrees freedom of the gauge bosons W ± and Z which in this way acquire mass.
While the SM provides a phenomenological description of electroweak symmetry breaking, an understanding of the underlying dynamics is lacking so far. Furthermore, a number of experimental observations -e.g., the presence of Dark Matter or the existence of non-zero neutrino masses -and theoretical issues -e.g., the hierarchy problemcannot be explained within the framework of the SM and are, therefore, hints for the existence of physics beyond the SM (BSM). Indeed, BSM models with an extended Higgs sector, featuring additional Higgs bosons, provide a description of the available data that is at least comparable or even better than for the case of the SM, see e.g. Ref. [3] for a recent discussion of LHC results and Refs. [4,5] for global fits.
So far no direct evidence for BSM particles has been found via the searches at the LHC, which instead have placed increasingly strong lower bounds on the masses of potential BSM particles. While this does not rule out the possibility of light BSM particles with relatively small couplings to SM particles, in many BSM scenarios that are currently investigated it implies a rather large hierarchy between the electroweak (EW) scale of the SM and the scale of at least some of the BSM particles.
For the description of BSM scenarios where the masses of the new particles are so heavy that they are beyond the reach of present and future accelerators, effective field theory (EFT) frameworks -e.g., SMEFT [6,7] -offer a conceptually clean way to parameterise the effects of heavy BSM particles on low-scale physics. Moreover, they provide a way to resum large logarithmic corrections involving the EW and the BSM scales appearing in the theoretical predictions for physics at the EW scale, such as the prediction for the mass of the SM-like Higgs boson in supersymmetric theories, see e.g. Ref. [8] for a recent review.
By construction, these EFT frameworks, where the heavy particles are integrated out, do not allow the description of the dynamics of heavy BSM particles. Such a description, however, is of interest since the LHC high-luminosity run and also future colliders like an e + e − linear collider or the FCC offer a significant potential for discovering new particles.
Precise theoretical predictions for the production and decay of new heavy particles are important for determining the viable parameter space of the considered BSM scenarios and for assessing the sensitivity for discoveries and for discriminating between different possible realisations of BSM physics. In those predictions, which have been obtained in the literature for a variety of models of BSM physics, the appearance of large logarithmic corrections potentially spoiling the reliability of the perturbative expansion is a re-occurring issue. For describing the dynamics of heavy BSM particles, it is, how-ever, obviously not possible to resum these logarithmic corrections by integrating out the heavy particles.
A well-known example in this context are the fermionic decays of a heavy BSM Higgs boson. In order to avoid large QCD corrections, the fermion mass entering via the Yukawa coupling of the Higgs boson to a pair of fermions should be evolved to the heavy Higgs mass scale [9,10]. Large logarithms develop, however, not only in this kind of QCD corrections, but also in the form of electroweak Sudakov logarithms (see e.g. Refs. [11,12]) which are related to the exchange of W , Z or light Higgs bosons. In principle, these electroweak Sudakov logarithms can be resummed using soft-collinear effective theory (SCET) [13][14][15]. A specific SCET approach for resumming large logarithms in the decay of BSM particles to SM particles has been developed in Ref. [16] and applied to various example models in Refs. [17][18][19].
In the present paper, we carry out a detailed investigation of the occurrence of large logarithmic corrections in processes with external scalar BSM particles. We point out that large logarithms can arise via external leg corrections involving an interaction between two heavy scalar particles and one light scalar particle. The size of these corrections can be further enhanced 1 by large trilinear couplings between the involved scalars. Focusing first on a simple toy model, we explicitly show the origin of these logarithms at the one-loop level. Moreover, we point out how these logarithms are related to infrared divergencies which can be cured by light scalar radiation or by resumming the mass corrections of the light scalar. While in principle one would expect that it should be possible to resum these Sudakov-like logarithms with the methods of e.g. Ref. [16] (the corresponding EFT realisation to the best of our knowledge has not been obtained in the literature so far), we follow a more direct approach in the present paper by explicitly calculating the two-loop corrections involving the large logarithms. In this context, we compare different renormalisation schemes for the involved parameters. As a result, we stress the importance of choosing on-shell-like schemes in order to avoid corrections that are enhanced by powers of the heavy BSM scale over the EW scale. In addition to these conceptual studies, we discuss various exemplary applications. First, we discuss the decay of a gluino into a top and a scalar top quark in the Minimal Supersymmetric extension of the Standard Model (MSSM). As a second example, we study heavy Higgs decays to two tau leptons in the next-to-minmal Two-Higgs-Doublet Model (N2HDM). For each of these examples, we explicitly derive the quantum corrections involving the large logarithms at the one-loop and the two-loop order.
Our paper is organised as follows. In Section 2, we introduce our toy model. The occurrence of large logarithms for this toy model is then discussed in detail in Section 3. The various exemplary applications are presented in Section 4. Section 5 contains the conclusions. Details on the two-loop calculations carried out in this paper are given in Appendices A and B.

Toy model
In order to illustrate our discussion as clearly as possible, we start by focussing on a toy model containing three real singlet scalars φ 1,2,3 that are coupled to a Dirac fermion χ. This model offers a simple setting for discussing the issue of large trilinear-couplingenhanced logarithms arising via external scalar leg corrections, and the different mass configurations in which they can appear. Moreover, it can be directly mapped to many BSM models (as done in Section 4 for some examples).
We endow this model with a Z 2 symmetry under which the scalars and fermions transform as We assume that the Lagrangian parameters are chosen such that this Z 2 symmetry is not broken spontaneously. Therefore, only φ 3 acquires a vacuum expectation value (VEV), which we denote as v 3 . We can, however, redefine the couplings in the Lagrangian in order to absorb this VEV v 3 (and thus work with a field φ 3 having φ 3 = 0) and we can furthermore assume without loss of generality that the scalar mass matrix is diagonal. The Lagrangian of this toy model can then be written as where the trilinear couplings A ijk will be of special interest in the discussion below. Note that since the fermions are not charged under the Z 2 symmetry, only φ 3 can couple to them. Furthermore, because the three scalars are real, the trilinear, quartic, and Yukawa couplings in the Lagrangian are also real. We are mostly interested in the situation in which φ 1 is much lighter than φ 2 and φ 3 , while φ 2 and φ 3 are approximately mass-degenerate. Accordingly, we will consider in the following the double limit in which m 1 m 2 ∼ m 3 and concurrently m 2 → m 3 . From a phenomenological point of view, φ 2 and φ 3 can be regarded as heavy BSM fields whereas φ 1 represents a SM field with a mass around the electroweak scale.

Scalar external leg corrections as an origin of large logarithms
One important ingredient in obtaining predictions for physical observables is to ensure that the external particles have the correct on-shell properties as required by the LSZ formalism [20]. For unstable particles that can mix with each other this can be achieved by employing a (in general non-unitary) matrix, denoted byẐ in the following (see e.g. Refs. [21][22][23]), which relates the (renormalised) one-particle irreducible vertex functions involving the external loop-corrected mass eigenstate fields φ physical a to the (renormalised) vertex functions involving the internal lowest-order fields φ j , The (aj) element of theẐ matrix can be written as a combination of the usual LSZ factor for the case with non-vanishing mixing and a term accounting for the mixing between the states i and j,Ẑ where no summation over repeated indices is implied. The terms on the right-hand side are given byẐ The superscript " " denotes a derivative with respect to the external momentum squared. The elements of the propagator matrix, ∆ ii (p 2 ), ∆ ij (p 2 ), . . . , are obtained from inverting the matrix involving the renormalised self-energies of the fields φ i , φ j , . . . , denoted aŝ Σ ii (p 2 ),Σ ij (p 2 ), . . . , where the diagonal and off-diagonal entries read p 2 − m 2 i +Σ ii (p 2 ) andΣ ij (p 2 ), respectively. Here m i is the tree-level mass of φ i . The effective self-energŷ Σ eff ii (p 2 ) is composed of the usual self-energyΣ ii (p 2 ) and mixing contributions. For the example of the case where the three particles i, j, k can mix with each other it readŝ The quantities in Eq. (5) are evaluated at the complex propagator pole, p 2 = M 2 a , which is associated with φ i .
We note in this context that the elements of theẐ matrix given above contain higherorder contributions that are crucial for the description of the resonant mixing of two or more unstable particles that are nearly mass-degenerate [22,23]. In the following discussion focusing on large logarithmic contributions arising from external-leg corrections of heavy scalar particles we will always perform a strict perturbative expansion, keeping only the terms contributing to the order -one-or two-loop -at which we are working. In this way a mixing of orders in perturbation theory is avoided, so that there is no associated residual dependence on unphysical gauge or field-renormalisation contributions (as discussed, e.g., in Refs. [24,25]). Thus, we do not provide here a treatment of the resonance-type behaviour of the nearly mass-degenerate heavy fields φ 2 and φ 3 that would in general be expected to occur (this refers in particular to parameter regions ϕ 3χ χ Figure 1.: φ 3 →χχ decay process at tree level. where their mass difference is smaller than the sum of their total widths). We note that in the considered toy model the mixing between the fields φ 2 and φ 3 is forbidden by the imposed Z 2 symmetry, but we will study a model with non-vanishing scalar mixing in Section 4.1.3. The treatment of the resonance-type behaviour, as outlined in [22,23], can be carried out in addition to the analysis of large logarithmic contributions that will be presented in the following.

Large trilinear-coupling-enhanced logarithms appearing in scalar external leg corrections
The external leg corrections can become a source of large logarithmic contributions. We discuss this here in the context of the toy model introduced in Section 2. The discussion is, however, straightforwardly transferable to other models (as we will demonstrate in Section 4).
As an example process, we investigate the φ 3 →χχ decay process, which is shown at tree level in Fig. 1. In more realistic models, this process could for example correspond to the decay of a heavy Higgs boson to two SM fermions (see also Section 4.2). The corresponding process for the other heavy field, φ 2 →χχ, is forbidden as a consequence of the imposed Z 2 symmetry (which also eliminates interference effects in the production and decay of φ 2 and φ 3 ). We note, however, that both the decay processes of φ 2 and φ 3 need to be taken into account in order to ensure the cancellation of infrared divergences in the limit where m 2 = m 3 and m 1 = 0, see the discussion below.
At tree-level, the decay width for this process reads Once radiative corrections are taken into account, one must compute additional contributions to this process, as illustrated in Fig. 2. We denote the relative one-and two-loop corrections to the decay width by ∆Γ (1) φ 3 →χχ and ∆Γ (2) φ 3 →χχ respectively, so that Among the radiative corrections, one must take into account one-particle-reducible diagrams corresponding to field renormalisation contributions of the external particlesor equivalently, to the LSZ factors on the external legs (see above).
Here, we focus on large logarithmic corrections that are proportional 2 to powers of the trilinear couplings A ijk . These arise only via the LSZ factor of the external φ 3 leg. Taking into account only terms proportional to at least two powers of trilinear couplings, we obtain for the one-loop corrections to the φ 3 →χχ decay process 3 where k ≡ (4π) −2 is used to indicate the loop order, B 0 is the usual Passarino-Veltmann function (we recall its definition in Eq. (81)), and the ellipsis denotes terms that are not proportional to at least two powers of trilinear couplings. As already noted in Section 2, we concentrate on the situation in which m 1 m 2 ∼ m 3 , and hence the second term in Eq. (9) is of particular interest because it is infrared divergent in the double limit of m 1 → 0 and m 3 → m 2 . At the one-loop level no other terms proportional to A 2 123 appear that could cancel this divergence. For this mass hierarchy, we can distinguish two different cases, In both cases, logarithms appear which become large in the limit m 1 → 0 or → 0. The suppression by 1/m 2 3 can be compensated by the prefactor A 2 123 if A 123 ∼ m 3 , a situation which can easily be realised in many BSM theories. If one of the scalar fields would obtain a vev, similar terms proportional to two powers of the vev times a scalar quartic coupling would appear. However, such terms are expected to be smaller than genuine BSM trilinear terms because they are of order vev times an O(1) numberthe size of the quartic coupling being limited by unitarity -while the natural scale for BSM trilinear couplings is the BSM mass scale itself. Note on the other hand that the magnitude of trilinear couplings can also be constrained by unitarity (at finite energies), see e.g. Ref. [28], or by other considerations like vacuum stability, as is for instance the case with X t in the MSSM, see for instance Ref. [29].
The described logarithms do not involve the renormalisation scale and can, therefore, not be avoided by an appropriate scale choice. Also a straightforward EFT approach, in which φ 3 and φ 2 are integrated out, would not circumvent the issue, since we consider here a process in which φ 3 appears as an external particle. Consequently, the appearing logarithms can not straightforwardly be resummed using renormalisation group equations. 4 We also want to remark that this type of divergence, caused by a light scalar close to or in the massless limit, is similar to the so-called "Goldstone boson catastrophe," which has been discussed in the literature -see in particular Refs. [30][31][32][33][34][35].

The infrared limit
As stated above, the external leg corrections to the φ 3 → χχ process develop infrared divergences in the limit m 1 → 0 or → 0, and this obviously raises the question of how to address them. In the following, we discuss two methods to cancel out the IR divergences: firstly the inclusion of real radiation, as prescribed by the KLN theorem [36,37], and secondly a resummation of φ 1 contributions. 5

Inclusion of real corrections
In the case where m 2 = m 3 and m 1 = 0, the radiation of a real φ 1 scalar becomes kinematically allowed, and the appearance of the IR divergence in the φ 3 → χχ decay width can be understood as being due to considering an IR-unsafe observable. Following the KLN theorem, in order to obtain an IR-safe observable, one should include also the soft part of the Bremsstrahlung process φ 2 → χχφ 1 , shown in Fig. 3, which cannot be experimentally resolved from the φ 3 → χχ decay. Reintroducing a mass m 1 as an IR 4 As already mentioned in the introduction, we expect the SCET approach worked out in Ref. [16] to provide a way to resum these logarithms. While the approach seems to be suitable, the needed concrete EFT has not been worked out yet. 5 For completeness, we should mention another possible approach, following Ref. [38], that is to include the width Γ 3 of the φ 3 scalar, and to evaluate the process at the complex pole, M 2 3 = m 2 3 − iΓ 3 m 3 . However, this approach for curing the IR divergence works only at the one-loop order if no non-zero width is taken into account for the internal particles. Moreover, in some models or scenarios, the width of φ 3 can vanish. regulator, we obtain where E denotes the energy resolution of a detector that would be used to study this process. We can observe that the second term in the brackets in the lower equation exactly cancels the divergence caused by d dp 2 B(p 2 , m 2 1 , m 2 3 )| p 2 =m 2 3 in the limit m 1 → 0c.f. Eq. (10) -and we find where the ellipsis denotes IR-safe terms, not involving A 123 . We note that the inclusion of the soft real radiation not only cancels the IR-divergent logarithm, but also changes the finite part of the result (c.f. the −1 piece in Eq. (12)). If we instead set m 1 = 0 and regulate the IR divergence by (= m 2 3 − m 2 2 ), the real radiation diagram in Fig. 3 evaluates to cancelling the IR divergent term in Eq. (11), and giving for the terms of order (A 123 ) 2 in the sumΓ(φ 3 → χχ) + Γ (0) (φ 2 → χχφ 1 ) soft the same expression as in Eq. (13).
Finally, we note that while including only the soft part of the real φ 1 radiation process introduces, as mentioned already, a dependence on a new parameter E , this can be avoided by computing also the hard part of the real radiation process. We will do this in the following in our numerical analysis, where we perform the necessary three-body phase space integration numerically. Figure 4.: Examplary Feynman diagrams appearing in the resummation of φ 1 contributions. The hashed dots denote the one-loop φ 1 self-energyΣ 11 .

Resummation of φ 1 contributions
Inspired by one of the proposed solutions to the "Goldstone boson catastrophe," a second possibility to address these problematic IR-divergent terms is to resum the contributions from φ 1 , following Refs. [30,33,35]. While the resummation approach was devised for the IR divergences from Goldstone bosons in the effective potential (and its derivatives), it can also be applied for a more general situation, such as the one here. The reason for this is that scenarios with large mass hierarchies lead to significant corrections to the light-scalar mass(es). This means that the perturbative expansion is not adequately organised in this case, because all diagrams with (subloop) self-energy insertions yield large contributions, and a resummation is necessary to capture the relevant contributions of higher-order diagrams with self-energy insertions. We note that the approach of setting the φ 1 (the Goldstone-like state in our toy model) on shell, as done in Ref. [33], would not work here as there is no tree-level contribution involving φ 1 that would produce a one-loop shift in the calculation of the φ 3 →χχ process in such a way as to cancel the IR divergences.
We show in Fig. 4 some examples of the diagrams that are resummed by applying this procedure. At the cost of mixing orders in perturbation theory, this resummation generates a finite mass for φ 1 and therefore avoids divergences in the derivatives of the B 0 function in Eq. (9). For m 1 = 0, the contributions to the φ 1 mass read whereΣ (1) 11 (p 2 ) denotes the one-loop φ 1 self-energy at external squared momentum p 2 , and the loop function A 0 (m 2 ) is defined in Eq. (81) below. We recall also that the quartic couplings λ 1122 and λ 1133 are defined in the Lagrangian in Eq. (2). The (A 123 ) 2 terms in the one-loop decay width are then obtained aŝ which is free of IR divergences. We should point out here that employing this resummation means that one interprets the occurrence of the IR divergence in the one-loop Figure 5.: φ 3 self energy at the two-loop level.
corrections to the φ 3 decay as being due to the breakdown of the perturbative expansion explained above, rather than a lack of inclusiveness of the observable we are computing. As we will see in the following section, at two loops, there are additional IR-divergent contributions that are not contained in the class of diagrams that are incorporated by the resummation as indicated in Fig. 4, and the resummation approach therefore would have to be extended for applications beyond one-loop order to include further classes of diagrams. It should furthermore be noted that this approach of curing IR divergences would have to be modified or may not be realised if a symmetry enforces the considered particle to be strictly massless (like e.g. in the case of the photon).

Scalar external leg corrections at the two-loop level
We now investigate the external leg corrections to the φ 3 → χχ decay process at the two-loop order. For this discussion, we set m 2 2 = m 2 3 = m 2 and m 2 1 = . Following the discussion at the beginning of Section 3, the external-leg corrections to the φ 3 → χχ decay width can be expanded up to two loops aŝ whereΣ 33 is the renormalised φ 3 self energy. In particular, the two-loop contributions are given by To leading powers in the trilinear coupling A 123 , the one-and genuine two-loop contri-butions to the φ 3 self-energy (see Fig. 5) read in terms of MS-renormalised parameterŝ The superscript "genuine" for the two-loop self-energy indicates that we consider therein only pure two-loop contributions. Definitions 6 of the two-loop T 11234 and T 12345 loop integrals as well as expressions for their derivatives (obtained in the MS scheme) can be found in Appendices A and B. The two-loop integrals contain potentially dangerous terms of the form 1/ , 1/ √ , ln 2 /m 2 , and ln /m 2 , which in turn give rise to potentially even more problematic terms of the form 1/ and 1/ √ (as well as ln 2 /m 2 , ln /m 2 ) in the first derivatives that enter the external leg corrections. Indeed, we have up to order where the results in Appendix B for the derivatives of the T integrals have been applied, and the notation lnx ≡ ln x/Q 2 (Q being the renormalisation scale) has been used. Applying Eq. (17), we obtain a pure MS expression for the dominant external-leg contributions in powers of A 123 to the φ 3 → χχ decay width, which readŝ We will use this pure MS expression for comparison in our numerical analysis below, but otherwise for obtaining our results we employ an on-shell renormalisation scheme -at least for the scalar masses. This implies that finite contributions from subloop renormalisation need to be included, which for the parameters in our result are of the formΣ where and m 2 2 , the one-loop counterterm for A 123 , and the one-loop field renormalisation counterterm for φ 3 , respectively.
We now derive the dominant contributions to these counterterms in powers of A 123 . It should be noted that because the expressions forΣ are MS -renormalised quantities, only the finite parts of the counterterms need to be considered here (as the UV divergences have already been cancelled). Adopting an on-shell renormalisation scheme for the scalar masses, the mass counterterms are Next, renormalising the φ 3 field in the MS scheme, we simply have It should be noted that the MS field renormalisation drops out in the sum of the vertex correction and the LSZ factor. Therefore the prescription that has been chosen for the field renormalisation is insignificant. Turning finally to the trilinear coupling counterterm δ (1) A 123 , there are several possible choices for the renormalisation of A 123 . We consider here the following options.
Taking now the derivative ofΣ with respect to p 2 , we can observe that the 1/ and 1/ √ divergent terms contained in these two terms exactly cancel with the ones in the derivative ofΣ . Indeed, the first term of Eq. (25) is which yields an exact cancellation of the 1/ and 1/ √ terms in the expression for d dp 2 T 11234 (p 2 , , , m 2 , m 2 , m 2 ). The second term of Eq. (25) is which yields a cancellation of the 1/ √ term in d dp 2 T 11234 (p 2 , m 2 , m 2 , , m 2 , ). Accordingly, we obtain d dp 2 T 11234 (p 2 , , , m 2 , m 2 , m 2 ) We note that the renormalisation-scale dependence has dropped out in this expression. Thus, for the two-loop external-leg corrections at leading order in A 123 , i.e.
O((A 123 ) 4 ), the incorporation of OS counterterms for the involved masses is already sufficient to obtain a result that is independent of the renormalisation scale. This fact can be explained by the finiteness of the A 123 counterterm at leading order in powers of A 123 . Indeed at this order, vertex corrections -in this toy model -can only involve scalar triangle diagrams, resulting in the scalar integral C 0 (defined in Eq. (82)) that is UV-finite, while field renormalisation contributions only involve derivatives of the B 0 function, which are also UV-finite (this is further illustrated in Eq. (31) below). In other words, even in the MS scheme, the contributions to the trilinear coupling A 123 are scale-independent at leading powers in A 123 .
The derivative of the total two-loop O((A 123 ) 4 ) contribution to the φ 3 self-energy hence readŝ Taking all pieces together, we can write the dominant contributions in powers of A 123 to the φ 3 → χχ decay width including external leg corrections up to two loops asΓ (ii) A second possible choice is to renormalise A 123 on-shell: for this we require that the on-shell-renormalised loop-corrected φ 2 → φ 1 φ 3 amplitude with on-shell momenta shall remain equal to the tree-level result. With this condition, we obtain and the one-loop scalar integral C 0 is defined in Eq. (82) below. Since this result contains only derivatives of the B 0 integral, it is manifestly UV-finite, in accordance with the above discussion concerning the scale-independence of A 123 (to leading powers in A 123 ). This result for δ (1) OS A 123 yields additional terms inΣ (2), subloop 33 . Accordingly, one finds for the derivative of the total two-loop φ 3 self-energŷ This yields for the decay widtĥ which as discussed above is independent of the choice of the renormalisation scale.
(iii) As a third possible choice, we explore whether it is possible to choose the finite part of the A 123 counterterm in such a way as to cancel the ln 2 term inΓ(φ 3 → χχ)we will denote this version of the counterterm δ (1) no-log-sq A 123 . We should emphasise that this only shuffles the ln 2 term away from the calculation of the φ 3 → χχ crosssection and into the extraction of A 123 from a physical observable (e.g. related to the φ 3 → φ 1 φ 2 process) and thus it does not remove the potentially large logarithmic term entirely. However, if A 123 is only considered as a user-given input, this choice may be convenient. The piece of the subloop renormalisation related to δ (1) A 123 is given by In this prescription we require that the term is cancelled. This yields Inserting this back intoΣ (2) 33 (m 2 ), we obtain and in turn for the decay width, It should be noted in this context that in general it is not possible to choose the counterterm δ (1) A 123 such that both the ln 2 and ln terms are cancelled.

Numerical analysis
In this Section we present numerical investigations of the external leg corrections at one and two loops, in both mass configurations specified in Eqs. (10) and (11). We start by showing in Fig. 6 the relative modification of the φ 3 → χχ decay width due to the inclusion of one-loop external φ 3 leg corrections as a function of the mass of the light scalar φ 1 . First, the red curve shows the one-loop result in the absence of any treatment of the infrared divergence, and therefore it diverges logarithmically if m 1 approaches zero, as expected from Eq. (10). Next, we compare the impact of the different treatments of the IR problem discussed above. The blue curves are obtained by including soft real radiation of a φ 1 scalar, displaying several values of the detector energy resolution E . While the IR divergence is cured independently of the numerical value of E , the one-loop contribution toΓ(φ 3 → χχ) is sensitive to it: for m 1 = 1 MeV, one finds a one-loop contribution of approximately −35% for E = 1 GeV, compared to −18% for E = 20 GeV. If the hard part of the real radiation is also included (performing the phase-space integration numerically), one finds the magenta dotted curve, which is independent of the value of E . In this case the one-loop contribution amounts to about −7% of the tree-level result. Finally, the two green curves indicate the result obtained after resummation of φ 1 contributions, for two possible choices of the renormalisation scale, Q = m 3 and Q = 2m 3 , entering the loop functions in ∆m 2 1 , see Eq. (15). One can see that for small values of m 1 the resummation procedure also cures the IR divergence in the external leg corrections. The difference between the curves -for m 1 = 1 MeV, we find a one-loop effect of about −12% for the curve with Q = m 3 , while it is only   −2.3% for Q = 2m 3 -can be interpreted as an indication of the theoretical uncertainty of this prescription. It can be associated in particular with diagrams containing higherorder subloop insertions of φ 1 that are not included in the resummation performed here. This large theoretical uncertainty of the resummation procedure indicates a limitation of applying this approach for the purpose of obtaining reliable higher-order predictions for the φ 3 → χχ decay width. While the different predictions for the results including real radiation can be understood as corresponding to different experimental situations, i.e. they refer to different physical observables, there is no clear interpretation of which observable the resummed prediction for the decay width should be compared toespecially as varying the renormalisation scale Q leads to significantly modified results. It is furthermore apparent that the curve for Q = m 3 exhibits a spike around m 1 42 GeV arising from the fact that the sum m 2 1 + ∆m 2 1 , appearing in logarithms of the loop functions, vanishes for this particular value of m 1 . Since the resummation procedure is meant to be applied to the region where m 1 approaches zero, this numerical artifact at a relatively large non-zero value of m 1 has no particular physical significance.
Next, in Fig. 7 we investigate the impact of the two-loop external leg corrections for the different possible choices of renormalisation of the trilinear coupling A 123 and of the scalar masses that were discussed in the previous section. We set the mass of the heavy scalars to be m 2 = m 3 = 1 TeV and take A 123 = 3 TeV. In the left plot, the one-loop (blue curve) and two-loop (red curves) results for the φ 3 → χχ decay width are shown relative to the tree-level result. For the two-loop results the OS scheme has been adopted for the scalar masses, and A 123 is renormalised in either the MS (dashed), the OS (solid), or the "no-log-sq" (dotted) scheme. The two-loop corrections employing either an OS or MS renormalisation for the trilinear coupling A 123 turn out to be numerically significant. For m 1 = 10 −3 GeV they amount to +25.4% and +21.6% in the OS and MS schemes, respectively, as compared to −73.0% at one-loop order, enteringΓ(φ 3 → χχ) with the opposite sign than the one-loop corrections. On the other hand, in the "no-log-sq" scheme -devised to eliminate the ln 2 term inΓ(φ 3 → χχ) -the two-loop corrections are noticeably smaller (only −7% for m 1 = 10 −3 GeV) and have the same sign as their one-loop counterparts. It should be noted, however, that the red curves cannot be compared directly to each other, since they correspond to different physical situations as a consequence of the different renormalisation schemes used (we have not incorporated the shifts in the numerical values of A 123 that are induced by the different renormalisation schemes).
In the right plot of Fig. 7 the two-loop corrections are shown relative to the one-loop result for the considered three different choices of the renormalisation of A 123 as well as for a complete MS renormalisation (of the masses and of A 123 ), using two different values of the renormalisation scale Q. While the three choices for the renormalisation of A 123 together with an on-shell renormalisation of the masses give rise to the moderate twoloop effects that are shown in more detail in the left plot, one can see that the two-loop corrections in the full MS scheme yield unphysically large effects. Those corrections are many orders of magnitude larger than at one loop and furthermore exhibit a huge scale dependence. These unphysically large effects arise as a consequence of the severe nature of the terms enhanced by 1/m 1 and 1/m 2 1 in the derivatives of the T 11234 integrals, see the expressions in Appendix B. Thus, our analytical and numerical results highlight the importance of adopting an on-shell renormalisation for the masses entering the computation at the one-loop order, in order to avoid unphysical power-enhanced terms at the two-loop level. This situation is similar e.g. to the known occurrence of non-decoupling effects in the context of Higgs-mass calculations in supersymmetric models that is encountered for scenarios with heavy gluinos if the parameters in the scalar top sector are renormalised in the DR scheme, see for instance Refs. [42][43][44][45][46].
As a final step of our numerical analysis of the toy model we show in Fig. 8 our results for the other mass configuration where m 1 = 0 and there is a mass splitting = m 2 3 −m 2 2 . In our one-loop analysis we had seen that the mass splitting = m 2 3 − m 2 2 acts as an IR regulator for this mass configuration. The results presented in Fig. 8 are obtained using the analytical expressions for the derivatives of two-loop self-energies given in Appendix B.4, and adopting an OS scheme for all scalar masses. Because both the twoloop integral T 11234 (p 2 , m 2 1 , m 2 1 , m 2 + , m 2 + , m 2 ) and its derivative with respect to p 2 are always IR divergent for m 1 = 0 (see Eq. (94) below), does not suffice to regulate all IR divergences at two loops for this mass configuration. We therefore introduce an IR regulator mass squared of m 2 IR = 10 GeV 2 for φ 1 in order to treat the otherwise divergent two-loop integral. We emphasise that this IR divergence that is caused by the occurrence of a squared massless propagator in T 11234 is different from the IR divergences that we had encountered in our one-loop analysis and also in the two-loop analysis for the   other mass configuration. This can been seen from the fact that the divergence appears already in the self-energy and not only in its derivative. In fact, this is precisely the situation encountered in the so-called "Goldstone boson catastrophe" discussed above. In the prediction for a physical observable these IR divergences are expected to cancel with the corresponding IR divergences in the real radiation contributions at the two-loop level in an analogous way to the one-loop case that we had discussed above. While the sum of those contributions should be free of logarithms involving m IR , it still contains logarithms involving the small but non-zero quantity that can be numerically large. On the lefthand side of Fig. 8, we plot, as a function of the ratio m 2 /m 3 , the relative deviation of the φ 3 → χχ decay width from its tree-level prediction at the one-loop (blue and green curves) as well as at the two-loop level (red curve). We choose here a parameter point for which m 3 = 500 GeV, m χ = 200 GeV, A 123 = 3m 3 , λ 1122 = 1, and λ 1133 = 1.2, while the other trilinear couplings are set to zero for simplicity. As expected, we find the external leg corrections to be divergent when m 2 /m 3 → 1 at one-and two-loop order. The two green curves show the one-loop result incorporating a resummation of φ 1 contributions for Q = m 3 (dashed green curve) and Q = 2m 3 (dotted green curve). As for the other mass configuration the resummation regulates the IR divergence, but again we find that the resummation procedure exhibits a large dependence on the renormalisation scale Q. This difference originates from the quite significant difference in ∆m 2 1 between the two cases. It stems from the dependence on Q in ∆m 2 1 , or in other words from higherorder contributions from φ 1 self-energy insertions. On the right-hand side of Fig. 8, we illustrate the size of the two-loop correction relative to the one-loop result. The more rapid divergence of the two-loop corrections is a consequence of the ln 2 term.
To summarize our study of the toy model, we have identified -at the one-and two-loop level -the potentially dangerous logarithmic and squared-logarithmic contributions arising from external-leg contributions to a process involving a heavy scalar as an external state. We have shown that in the IR limit, in which these terms become divergent, the inclusion of real radiation cures the IR divergences and yields a finite and well-defined result. However, should the experimental resolution allow the separation of the real-radiation contributions (this could happen if either of the mass parameters that we defined as has a non-zero, but relatively small, value), large logarithmic and squared-logarithmic terms remain in the computed decay width. We have also illustrated how a resummation of higher-order contributions involving the light scalar can in principle alleviate the IR problem, although this method is in the present calculation plagued by large theoretical uncertainties arising from the problem of relating the prediction to a well-defined physical observable. Furthermore, we have discussed several possible choices of renormalisation schemes for the parameters entering this calculation, and their impact on the size of the remaining large logarithms. In this context, we pointed out the importance of renormalising the involved masses on-shell in order to avoid unphysically large power-enhanced corrections. At the same time, we found that, within the class of corrections that we considered, the renormalisation scale dependence of the contributions to A 123 vanishes, and the differences between the discussed options for renormalising this coupling are numerically moderate. While the two-loop effects remain smaller than those at one-loop order, their impact can become significant for large gt L,R t Figure 9.: Decay of a gluino into a stop and a top quark mass hierarchies, i.e. close to the IR limit. This is in particular the case if the MS or OS schemes are adopted for the renormalisation of A 123 , while as expected a scheme that was specifically devised to absorb squared-logarithmic terms yields the smallest effects at two-loop order.

Applications
In this Section, we discuss several concrete BSM models where large logarithms appear in the external leg corrections in an analogous way to the toy model as described in Section 3.

MSSM
As a first example, we consider corrections to processes involving external scalar top quarks, i.e. the superpartners of the top quark, in the MSSM. One such process is the decay of a gluino into a stop and a top quark as shown Fig. 9. The following discussion is straightforwardly transferable to other processes involving scalar top quarks on an external leg. 7 In the limit of vanishing electroweak gauge couplings, the stop mass matrix is given by where mt L and mt R are the stop soft SUSY-breaking mass parameters, m t is the topquark mass and X t ≡ A t − µ/ tan β (A t is the trilinear stop coupling, µ is the Higgsino mass parameter, and tan β ≡ v 2 /v 1 is the ratio of the vacuum expectation values of the two Higgs doublets). Here (and also for the rest of this Section), we assume all parameters to be real for simplicity. 7 The same type of corrections appears in the widely used OS scheme of the trilinear stop coupling

Y t terms
First, we aim at calculating all corrections to the external stop leg that are leading in powers of the trilinear coupling Y t ≡ A t + µ tan β, which couples the heavy Higgs bosons H, A, and H ± to the stops. We assume the heavy Higgs bosons to be much heavier than the electroweak scale but still lighter than the stops. In order to simplify our discussion, we work in the approximation of vanishing electroweak gauge couplings. We neglect all terms that are proportional to the electroweak scale by setting the vacuum expectation value v 2 ≡ v 2 1 +v 2 2 to zero. In this approximation the stops do not mix, since the off-diagonal entry of the stop mass matrix is zero (i.e., m t = 0 in this limit). Consequently, the left-and right-handed stop chirality eigenstates (t L andt R ) are also mass eigenstates.
In this approximation, all relevant couplings between the heavy Higgs bosons and the stops are given by where h t is the MSSM top-Yukawa coupling, and c β ≡ cos β. Note that the couplings of the charged Higgs bosons involve onlyt R . As a consequence of the coupling structure, the off-diagonal stop self-energy is zero, The diagonal left-and right-handed stop self-energies receive contributions from the H and A bosons (see the left plot of Fig. 10). The diagonal right-handed stop self-energy receives in addition a contribution from the charged Higgs boson (see the right plot of Fig. 10). The leading contribution in powers of Y t to the stop self-energies at the one-loop level is given byΣ where we set mt L = mt R = M SUSY , and m A is the mass of the heavy Higgs bosons. For the leading genuine two-loop corrections, we obtain The subloop renormalization is given bŷ where δ (1) m 2 A is the one-loop mass counterterm of the heavy Higgses, δ (1) m 2 t L ,b R ,t R are the one-loop mass counterterms of the stops and sbottoms, δ (1) Zt L,R are the one-loop field renormalisation constants of the stops, and δ (1) (h t c β Y t ) is the one-loop counterterm of the Higgs-stop-stop interaction.
Renormalising all counterterms in the OS scheme (apart from the field renomalisation), we obtain with ReΣ (1) The counterterm for the coupling h t c β Y t is fixed by demanding that the amplitude for the processt LtR → H up to the one-loop level should coincide with the tree-level result (taking into account contributions at leading order in Y t ). All genuine vertex corrections to this process vanish at this order, such that only the external leg corrections enter the definition of the counterterm. Using instead for the renormalisation of Y t e.g. the DR scheme would not significantly change the numerical size of the two-loop corrections (see the discussion in Section 3.3). We employ the DR scheme for all field renormalisation constants. As already mentioned above, this choice is insignificant, as the DR field renormalisation constants drop out in the sum of the vertex corrections and the LSZ factor. Using these expressions, we obtain the following results for the two-loop stop selfenergies (Σ (2) =Σ (2,genuine) +Σ (2,subloop) ), where we dropped the "Re" in the counterterm expressions of Eqs. (45a) and (45d) since all involved loop functions are real for the considered mass hierarchy. The leading corrections in powers of Y t to the gluino decay width are then obtained by taking into account the external stop LSZ factor (see Section 3). For an externalt L we obtain where on the right-hand side of the second equality we performed an expansion in m A /M SUSY and Y t ≡ Y t /M SUSY . For an externalt R we obtain Note that even though the one-loop corrections to the decay widths witht L andt R differ from each other, the logarithmic terms in the two-loop corrections are identical. We show an exemplary numerical evaluation for these corrections in Fig. 11. For this example, we fix Y t = √ 6, m A = 500 GeV, t β ≡ tan β = 2, and mg = 1.2M SUSY . We vary M SUSY and show the relative deviation from the tree-level results for the gluino decay widths including the one-(blue curves) and two-loop (red curves) corrections in leading powers of Y t . For the decay intot L and a top quark (solid curves), the one-loop corrections lower the tree-level decay width by up to ∼ 4% in the considered parameter region. We find the two-loop corrections to have a significantly smaller effect. The inclusion of the one-loop corrections has a larger effect for the decay intot R and a top quark (dashed curves) reaching values of up to ∼ 8%. This is a consequence of the additional factor of two appearing in the one-loop corrections of Eq. (49) with respect  to Eq. (48). Also for the gluino decay intot R and a top quark, the two-loop corrections have much smaller effect of 0.5% in this case. The example of the corrections in the MSSM that are enhanced by powers of the trilinear coupling Y t is a realisation of the situation of the toy model in Section 3.3 for m 2 2 = m 2 3 = m 2 and m 2 1 = . In the MSSM example the mass m 1 of the light scalar in the toy model corresponds to the common mass scale of the MSSM Higgs bosons H, A, and H ± , which has been set to ∼ 500 GeV in the present example. The heavy scale m of the toy model corresponds to mt L = mt R = M SUSY . If a gluino is detected in future searches in the decay modesg → t +t L,R , there should be significant experimental sensitivity for distinguishing the cases with and without additional H, A, or H ± radiation for this process. Thus, the appropriate theoretical description would not be the IR limit discussed above but rather the case illustrated here where the prediction for a measurable physical observable contains potentially large logarithms containing the ratio of the widely separated scales M SUSY and m A .
Regarding the encountered numerical effects, we find these potentially large logarithms to have a sizeable impact at the one-loop level. At the two-loop level, the size of the corrections is moderate for the considered example suggesting that a fixed-order treatment is sufficient and that no resummation of large logarithmic corrections is needed.

X t terms -case 1
Next, we look at the corrections leading in powers of X t with X t = A t − µ/t β . As for corrections leading in Y t , we work in the unbroken phase of the theory (i.e., v = 0). In this limit, X t only appears in the couplings of the light CP-even Higgs boson h, the neutral Goldstone boson G, and the charged Goldstone bosons G ± to stops and sbottoms, As a direct consequence of this coupling structure, we obtain as for the case above The limit v = 0 also implies that h, G, and G ± are massless. This situation corresponds to the second mass configuration discussed in Section 3, where the two heavy particles t L andt R are separated by a small mass difference , while the light particle is massless. Accodingly, we take mt L = M SUSY , m 2 t R = M 2 SUSY + and then consider the limit → 0. The stop self-energies obtain corrections proportional to X 2 t from diagrams analogous to the ones in Fig. 10 with H, A, and H ± replaced by h, G, and G ± , respectively. These corrections are given byΣ (1) For the subloop renomalisation, we obtain Here, we also introduce mass counterterms for the light scalars. In contrast to the broken phase of the MSSM, in which the mass counterterms of the light Higgs boson h (and also the Goldstone bosons) are dependent quantities that can be expressed in terms of the counterterms of the independent parameters m A , tan β and m Z , in the unbroken phase, independent mass counterterms can be introduced (by introducing counterterms for the bi-linear potential parameters). It should also be noted that in the unbroken phase the Goldstone bosons are physical particles with a finite mass. We fix all counterterms in the OS scheme. The counterterm for the coupling combination h t s β X t is fixed by demanding that the one-loop matrix element for the h →t LtR process is equal to its tree-level value. This yields It is important to remark that the counterterm . This is in contrast to the broken phase of the MSSM in which the loop corrections to the h tree-level mass are of O(v 2 ). 8 Since we assume that the physical masses of h, G, and 8 In the broken phase, the O(M 2 SUSY ) contributions of the hh self energy are cancelled as a consequence of the renormalisation procedure.
G ± are zero, we need to cancel the O(M 2 SUSY ) one-loop corrections by choosing the mass counterterm appropriately.
Summing the genuine and the subloop two-loop contributions, we get The leading X t corrections to the gluino decay widths up to the two-loop level are then given by (in the limit → 0 and m IR /M SUSY → 0), and where X t ≡ X t /M SUSY . Taking into account all two-loop corrections in the considered approximation, all infrared divergences ∝ 1/m IR or ∝ 1/m 2 IR cancel. In order to cancel the remaining ln m 2 IR terms, additional real light scalar radiation contributions would have to be taken into account, as explained in Section 3.3. This real radiation contribution would not introduce any additional large logarithms. Therefore, we set m 2 IR = for our numerical results presented below.
We assess the numerical size of these corrections in Fig. 12 showing the relative deviation from the tree-level results for the gluino decay widths including the one-(blue curves) and two-loop (red curves) corrections in leading powers of X t as a function of mt R /mt L . For the gluino decay intot L and a top quark (solid curves) the effect of the loop corrections is relatively modest (∼ 5%) for |1 − mt R /mt L | 0.05. If the two stops are nearly mass degenerate -corresponding to ∼ 0 -the external leg corrections are numerically large. As discussed in Section 3.2, the real radiation contributions should be included if the additional final state particles cannot be resolved experimentally. Close to mt L ∼ mt R also the two-loop corrections have a sizeable impact, reducing the overall size of the loop corrections by up to ∼ 10%. For the gluino decay intot R and a top quark (dashed curves) the overall size of the one-loop corrections is enhanced because of the additional factor of two in the one-loop part of Eq. (58) with respect to Eq. (59). As for the previous MSSM example we find that large logarithms occur in the prediction for the decay width of the gluino, which apart from the parameter region where real radiation contributions should be included are theoretically well under control in the two-loop fixed-order result.

X t terms -case 2
As a final MSSM example, we show that large logarithms also appear in the broken phase of the theory (i.e., if v = 0). Here, we assume that mt L = mt R . In the broken phase, the off-diagonal terms of the stop mass matrix (see Eq. (39)) are non-zero. After g → t +t L,R leading X t terms (case 1) g → t +t R @ 2L Figure 12.: The gluino decay widthsg → t +t L,R as a function of mt R /mt L , where the final state t +t L is displayed by the solid curves and t +t R by the dashed curves. The one-loop results including all corrections in leading powers of X t are shown relative to the tree-level decay width by the blue curves, while the corresponding two-loop results are displayed by the red curves.
diagonalisation, we denote the mass eigenstates byt 1 andt 2 with mt 2 ≥ mt 1 . Their couplings involving X t are given by 9 In contrast to the discussion in Section 4.1.2, for the case considered here the couplings are not the only source of the X t dependence. This is caused by the fact that also the stop masses depend on X t . The difference between the two stop masses is given by m 2 t 2 − m 2 t 1 = 2m t X t . Consequently, the stops become mass-degenerate in the limit v/M SUSY → 0, and large logarithms appear in analogy to Section 4.1.2.
We restrict ourselves here to a discussion at the one-loop level. The stop self-energy contributions containing h, G and G ± yield at leading order in X t of X t are then given bŷ It is important to remark that here the infrared divergences related to the occurrence of the massless Higgs boson h do not cancel (see below regarding the Goldstone contributions). This issue can be addressed by one of the strategies discussed in Section 3.2: either external h radiation could be included summing over different final states, or the higher-order mass corrections for the Higgs boson h could be included in analogy to the Goldstone boson resummation.
We focus here on the inclusion of real radiation. Having a close look at the coupling structure in Eqs. (60a) and (60f) and at the diagonal stop self energies in Eq. (61a), it is easy to identify the diagrams with an internal h boson and the same internal stop as the external stops as origin of the infrared divergence. Following the discussion in Section 3.2, we obtain for the gluino decay into a top quark, a stop, and the Higgs boson h in lowest order The infrared-divergent term exactly cancels with the corresponding term in Eq. (62) rendering the combinationΓg →t+t 1,2 + Γ (0) g→t+t 1,2 +h free of infrared divergences. It should be noted that the inclusion of real radiation of the neutral and charged Goldstone bosons would not yield additional infrared divergences. Accordingly, the virtual contributions of the neutral and charged Goldstone bosons also do not give rise to large logarithms (see Eq. (14) with = 2m t X t and m 3 = M SUSY ).
Taking into account h radiation, we show the numerical result for the decay width of the gluino into to the lighter stop and a top quark in Fig. 13. We choose X t = √ 6 and t β = 10. Three different values are chosen for the detector resolution: E = 1 GeV (red), E = 10 GeV (green), and E = 100 GeV (blue). The relative difference between the one-loop corrected and the tree-level decay width increases logarithmically for increasing g → t +t 1 (+h) leading X t terms (case 2) E = 1 GeV E = 10 GeV E = 100 GeV Figure 13.: The gluino decay widthg → t +t 1 (+h) at one-loop order including the corrections in leading powers of X t is shown relative to the tree-level decay width as a function of M SUSY . The result is shown for different detector resolutions: E = 1 GeV (red), E = 10 GeV (green), and E = 100 GeV (blue).
M SUSY , reaching values of up to ∼ −40%. Varying the detector resolution by two orders of magnitude corresponds to a variation of the decay width of ∼ 10%.
As an alternative to including real radiation diagrams, one could consider to take into account terms proportional to the electroweak gauge couplings. This would render the Higgs and Goldstone boson masses non-zero. While in this case no infrared divergences would appear, large logarithms would still occur in the limit v/M SUSY → 0, which corresponds to taking simultaneously the limits m 1 → 0 and m 2 → m 3 in Eq. (11).
Similarly to the examples discussed above, we expect the two-loop corrections to be much smaller than the one-loop corrections evaluated here and thus a fixed-order treatment should be sufficient for obtaining a percent-level precision.

N2HDM
We now turn to a second example of a frequently discussed BSM model in which potentially large logarithms arising from wave-function normalisation contributions may appear, namely the N2HDM, i.e. a real-singlet extension of the Two-Higgs-Doublet Model. We consider once again corrections to processes involving a heavy scalar on an external leg. Specifically, we discuss the decay process h 3 → τ + τ − , where h 3 is the heaviest of the three CP-even mass eigenstates (which we will assume to be mostly doublet-like).

The model
We start with a CP-invariant tree-level potential written in terms of two SU (2) L doublets Φ 1 , Φ 2 of hypercharge 1/2 and a singlet Φ S as For this model, we derived a FeynArts model file using SARAH [53][54][55].
In the following we will use the shorthand notations for combinations of the trilinear couplings that appear in the calculated self-energies, where tan β denotes the ratio of the vacuum expectation values of the two Higgs doublets.

4.2.
2. X 4 a external leg corrections to the h 3 → τ + τ − process Restricting ourselves to contributions involving only powers of X a -defined in terms of Lagrangian trilinear couplings in Eq. (68) -we find for the contributions to the h 3 self-energy at one-loop order where α 3 is the third CP-even mixing angle (see e.g. Eq. (11) in Ref. [56] for its definition). At two-loop order we find We adopt an on-shell renormalisation scheme for the masses and the trilinear coupling, together with an MS renormalisation of the h 3 field (which drops out in the sum of the vertex corrections and the LSZ factor). We then obtain and that involve heavy states (we drop terms involving only light states or involving c 2 2β ). In order to investigate the potentially large logarithms arising from the external leg correction in the prediction for the h 3 → τ + τ − decay width we consider a mass hierarchy where the first two CP-even mass eigenstates as well as the Goldstone bosons are light, 10 with m 2 At the two-loop level, for reasons of clarity we study the terms with different powers of s α 3 separately. Starting with terms of order s 4 α 3 , we havê At order s 2 α 3 , we find Finally, the terms without s α 3 are given bŷ Finally, the loop-corrected decay width for the h 3 → τ + τ − process, incorporating the external leg contributions, is found via the relation For the sake of brevity, we do not present the complete result, which can be straightforwardly obtained from the expressions in the previous equations.

Numerical results
We now turn to some numerical examples in order to illustrate the possible size of the external leg corrections computed in the previous section. In Figs. 14 and 15 we show the relative size of the one-loop (blue curve) and two-loop (red curve) external-leg corrections to the h 3 → τ + τ − decay width with respect to the tree-level result as a function of the heavy mass scale m. For these plots a light mass scale of √ = 50 GeV and a rather large trilinear coupling X a = 3m have been chosen. For s α 3 (and tan β) we consider two different scenarios: for the first one -inspired by one of the benchmark points in Ref. [57] (see table 5 therein) -we take s α 3 = 0.94. This choice ensures that the factor c α 3 , entering with the same power as X a , is not too small so that an overall suppression of this class of contributions is avoided. As can be observed in Fig. 14, this scenario indeed gives rise to sizeable contributions at the one-loop and the two-loop level. For m = 1 TeV, we find effects of −8.2% at the one-loop and −0.84% at the two-loop level.
The corrections grow to −27.1% and −8.0% at one-loop and two-loop order, respectively, for m = 100 TeV. Values of s α 3 closer to unity are however more easily reconcilable with the experimental data on the Higgs signal at about 125 GeV, and we therefore consider a further scenario with s α 3 = 0.99, making use of the public code ScannerS [58]. As shown in Fig. 15, in this case the size of the external leg corrections relative to the tree-level prediction for the h 3 → τ + τ − decay width is (as expected) much smaller than for the first scenario. At one-loop order, the corrections increase from approximately −1.4% for m = 1 TeV to −4.6% for m = 100 TeV, while at the two-loop level the corrections grow from −0.025% to −0.24% in the same range of m. The large relative increase of the two-loop corrections for increasing m is due to the appearance of ln 2 m 2 / terms at the two-loop level.
We briefly comment in this context on the constraints that have been applied for the above two scenarios. The analysis of Ref. [57] was carried out for an N2HDM without trilinear couplings, and the same is also true for the implementation of the N2HDM in the code ScannerS. Concerning the theoretical and experimental constraints that the considered parameter points have to fulfill, we note that theoretical properties such as boundedness-from-below or perturbative unitarity (in the high-energy limit) are not significantly altered by the inclusion of trilinear couplings in the N2HDM scenario that we are considering here. On the other hand, the trilinear couplings can modify the production cross-sections and decay widths, such that a more thorough analysis would be needed in particular for the scenario with s α 3 = 0.94 in order to assess its compatibility with the latest experimental results.
In comparison to the MSSM scenarios discussed in the previous sections, we have found that the relative size of the two-loop corrections is larger in the considered N2HDM scenario. Nevertheless we have found also for the case of the N2HDM that the perturbative expansion remains well-behaved even if the mass ratio entering the logarithmic contributions in the external leg corrections becomes relatively large. In particular, the two-loop corrections to the h 3 → τ + τ − decay width remain much smaller than their one-loop counterparts. Thus, also for the N2HDM case it appears safe to rely on a fixed-order analysis, and there does not seem to be a pressing need for developing a SCET-like resummation of the logarithmic contributions.

Conclusions
If a new heavy BSM particle (or more than one) is discovered, the characterisation of its properties will be of foremost importance in order to unravel the underlying structure of the observed physics beyond the SM. A crucial ingredient in this context will be precise theoretical predictions to which the experimental measurements can be compared.
One re-occurring issue when calculating loop corrections for heavy BSM particles is the appearance of potentially large logarithmic corrections. In this work, we have pointed out the existence of large Sudakov-like logarithmic contributions related to externalleg corrections of heavy scalar particles in scenarios where in addition at least one light scalar particle is present, giving rise to the possibility of large trilinear couplings between the scalars and a significant hierarchy between the masses of the different scalars.
Working in a simple toy model, containing one light and two heavy scalars coupled to each other with a potentially large trilinear coupling, we discussed in detail how the occurrence of the large logarithmic corrections is related to a limit in which infrared singularities emerge. We showed how these singularities can be regulated by including the radiation of the light scalar particle or by resumming higher-order contributions involving the light scalar. On the other hand, if for the decay of a heavy scalar BSM particle the final states with and without the radiation of the light scalar particle can be experimentally distinguished from each other, potentially large logarithmic corrections involving the ratio of the two widely separated mass scales occur in the prediction for the decay width of the heavy particle.
In order to assess the size of these logarithmic corrections at higher orders, we derived the leading two-loop external-leg corrections in the considered toy model. In this context, we compared different choices of renormalisation schemes for the parameters entering our computation (masses and trilinear coupling). As a result, we found that the choice of an MS renormalisation for the involved masses leads to unphysically large corrections that are enhanced by powers of the BSM scale over the electroweak scale. We showed that those huge effects are absent if the involved masses are renormalised in the OS scheme. In the OS scheme, we found the numerical size of the two-loop corrections to be moderate. This implies that a resummation of the Sudakov-like logarithms, which should in principle be possible in SCET-like frameworks, is not expected to be mandatory for obtaining theoretical predictions of sufficient accuracy.
As a further illustration of the qualitative features that we had found in our analysis of the toy model we then considered specific examples of models that are popular for studying possible BSM phenomenology. In the MSSM, we investigated corrections related to external scalar top quarks as appearing e.g. for the decay of a gluino into a top quark and its scalar superpartner. As a second example, we discussed the decay of a heavy Higgs boson into two tau leptons in the N2HDM. For both models, we found large one-loop corrections but only moderate effects at the two-loop level. Thus, the analyses in the specific models confirm our results for the toy model indicating that the class of logarithmic contributions that we have investigated here is not expected to spoil the perturbative expansion. The corrections that we have obtained at the two-loop order turned out to be much smaller than their one-loop counterparts, suggesting that a fixed-order analysis at this level should be sufficient for obtaining reliable theoretical predictions. with the N2HDM model file. We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306. H.B. acknowledges support by the Alexander von Humboldt foundation.

A. Definition of loop functions
In this Section, we list the definitions of the one-and two-loop integrals that have been used in this paper. We first introduce the shorthand notation where UV = (4 − d)/2, and µ is the regularisation scale (related to the renormalisation scale by the relation Q 2 = 4πe −γ E µ 2 ). We will also use throughout the following discussion.
We would like to draw the reader's attention to our conventions for denoting loop integrals and loop functions: in the following, we refer to the (unrenormalised) complete loop integrals with bold capitals, such as A, B, T, while we use normal script -A 0 , B 0 , T -for the finite parts of these integrals, defining the loop functions that appear in our computations.
First, the one-loop functions A 0 (x) and B 0 (p 2 , x, y) are defined as the UV-finite parts of the integrals as while the function C 0 is equal to the finite integral . (82)

B. Derivatives of two-loop functions
We present here our derivation of expressions for the derivatives of the two-loop selfenergy diagrams T 11234 and T 12345 -shown in Fig. 5 -which contribute to external leg corrections of heavy particles and contain potentially large terms as a consequence of the mass hierarchy between the involved particles.

B.1. Setup of the calculation and intermediate results
Our calculation makes use of the results of Refs. [40,41,59], in which T 11234 and T 12345 correspond to V and M , respectively. We start with the expressions of derivatives of the two-loop self-energy basis integrals with respect to the external momentum as well as mass arguments -see in particular Eqs. We require a number of limiting cases of the loop functions entering the relations for the derivatives: first, with one single mass scale we need where once again we have retained the external momentum argument explicitly (here p 2 = x). Next, we also use limits with two distinct mass scales, T 134 (0, x, y) = − I(0, x, y) T 234 (x, y, x, y) = − S(x, y, y, x) T 2234 (s, x, x, 0, 0) = T (s, x, 0, 0) T 2234 (x, x, x, y, y) = T (x, x, y, y) T 1234 (x, 0, y, x, y) = U (x, y, 0, y, x) Here Li 2 denotes the dilogarithm, and we recall that ζ(2) = Li 2 (1) = π 2 /6. In addition to these limits, we also had to derive several expansions of loop functions in powers of , such as The expansions of the B 0 and T 134 integrals can be obtained directly from the known analytical expressions of these functions -the former can be reproduced with Package-X [60], while the latter can be verified to correspond up to order to Eqs. (3.15) and (3.21) in Ref. [32]. For the T 1234 and T 12345 integrals, we can use Eqs. (3.20), (3.21), (3.22), and (3.32) in Ref. [41] to derive differential equations in the variable . The solutions of those differential equations yield the desired integrals (as we know their expressions for = 0).  We compare values obtained using the program TSIL [59] with those from our approximate expansion, at different orders in . Note that as the integral T 12345 is UV finite, there is no Q dependence in the corresponding loop function.

B.2. Example derivation
In order to illustrate the procedure of our calculations, we present here how we derived the expansion of T 134 (x, x, ) in powers of . Starting from Eq. (5.3) in Ref. [41], which gives the first derivative 11 of the integral I with respect to its first mass argument, and substituting the mass arguments with the two variables x and y, we find This is simply a differential equation in the function T 134 , for which we also know the boundary condition at y = 0 One can then straightforwardly solve this equation, and while the resulting function is fairly complicated, it can be expanded to arbitrary order in powers of y = x. We should emphasise here that for the case of the T 134 integral one could in principle obtain the same result by directly expanding the analytical expression that is known for general mass assignments. However, we discuss the case of this integral as it provides a simple example of our setup to derive expansions of the more complicated integrals T 1234 and T 12345 -for which analytical results are not known in general. Finally, for the derivatives of T 11234 and T 12345 , we make use of the relations in Sec. IV of Ref. [41].  11 Ref. [41] provides complete results for the derivatives of a basis of two-loop self-energy integrals with respect to external momentum and mass arguments. Solving this system of differential equations allows a numerical evaluation of those basis integrals. This method is implemented in the public tool TSIL [59].
We present in Table 2 some example values of these three derivatives, multiplied by m 4 in order to obtain a dimensionless number, for the set of mass parameters (a) in Eq. (89). We compare the results from TSIL (left column) to the values obtained with the approximate expressions of Eq. (92) (right column). We note that because TSIL provides results for the integrals V and M -corresponding to T 11234 and T 12345 in our notation -we need to take the derivatives with respect to p 2 numerically. We choose for this a step size of δ = 10 −3 GeV 2 , although we have verified that the numerical derivatives are stable and depend only very slightly on δ. Given that we derived approximate results only to order O( 0 ), we choose a set of mass inputs with small , for which we expect the higher-order terms to be moderate. Indeed, we find good agreement between the numerical results of TSIL and the approximate ones, with differences of only 0.06%, 3 · 10 −4 %, and 2%, respectively.