The triple Higgs coupling: A new probe of low-scale seesaw models

The measure of the triple Higgs coupling is one of the major goals of the high-luminosity run of the CERN Large Hadron Collider (HL-LHC) as well as the future colliders, either leptonic such as the International Linear Collider (ILC) or hadronic such as the 100 TeV Future Circular Collider in hadron-hadron mode (FCC-hh). We have recently proposed this observable as a test of neutrino mass generating mechanisms in a regime where heavy sterile neutrino masses are hard to be probed otherwise. We present in this article a study of the one-loop corrected triple Higgs coupling in the inverse seesaw model, taking into account all relevant constraints on the model. This is the first study of the impact on the triple Higgs coupling of heavy neutrinos in a realistic, renormalizable neutrino mass model. We obtain deviations from the Standard Model as large as to $\sim +30\%$ that are at the current limit of the HL-LHC sensitivity, but would be clearly visible at the ILC or at the FCC-hh.


Introduction
The CERN Large Hadron Collider (LHC) was home to one of the biggest discoveries in particle physics with the observation of a Higgs boson with a mass of around 125 GeV in 2012 [1,2], thanks to the data collected in Run 1 at 7 and 8 TeV. The Higgs boson is the remnant of the electroweak symmetry-breaking (EWSB) mechanism [3][4][5][6] that generates the masses of the other fundamental particles and unitarizes the scattering of weak bosons [7,8]. The Run 2 data collected in 2015 and 2016 at a center-of-mass energy of 13 TeV still displays a compatibility of this Higgs boson with the Standard Model (SM) hypothesis; nevertheless we know that the SM cannot be the ultimate theory. In particular the observation of neutrino oscillations, confirmed in 1998 at Super-Kamiokande [9], implies that neutrinos are massive, which cannot be explained in the SM framework and thus calls for an extension of the SM. One of the simplest possibilities to explain the non-zero neutrino masses and mixing is to add fermionic gauge singlets that will play the role of right-handed neutrinos. The addition of these heavy sterile neutrinos leads to the type I seesaw model and its various extensions [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. A very recent study summarizes the possible direct detection possibilities and indirect tests for heavy sterile neutrinos at lepton-lepton, proton-proton and lepton-proton colliders [26], see also references therein.
In a recent article [27] we have presented the triple Higgs coupling λ HHH as a new observable to test neutrino mass generating mechanisms in a regime of mass difficult to probe otherwise. The measure of λ HHH is one of the main goals of the high-luminosity run of the LHC (HL-LHC) as well as of the future colliders, such as the electron-positron International Linear Collider (ILC) [28] or the Future Circular Collider in hadron-hadron mode (FCC-hh), a potential 100 TeV pp collider (for the Higgs studies see reviews in refs. [29][30][31]). It would be a direct probe of the shape of the scalar potential that triggers EWSB. Any deviation of this coupling from the SM prediction is then welcomed to unravel new physics. In ref. [27] the study of neutrino effects on λ HHH was done in the context of a simplified model with the SM plus one heavy Dirac neutrino. It was found that effects as large as +30% at one-loop could be obtained, at the limit of the currently foreseen ∼ 35 % sensitivity that the HL-LHC will have to the SM triple Higgs coupling, when combining ATLAS and CMS data [32], but clearly measurable at the ILC [33] or the FCC [34]. A comprehensive study in a realistic and renormalizable model of neutrino masses was still left to be done.
In this article, we fill the gap and present the first analysis of Majorana neutrino effects on λ HHH . We work within the inverse seesaw (ISS) model [17][18][19], a renormalizable low-scale seesaw model generating neutrino masses. After taking into account all relevant constraints, we obtain effects that can be as large as a ∼ +30% increase of λ HHH , similar to the effects that we found in our previous article [27] using a simplified model. In the case of the ISS model, more heavy neutrinos are present, enhancing the effects as we expected, but the constraints on the model are stronger, reducing the end-effect back to the simplified model expectations. This can be clearly measurable at the ILC and at the FCC-hh and is at the limit of the currently foreseen sensitivity of the HL-LHC.
This article is organized as follows. In Section 2 we present the ISS model as well as the theoretical and experimental constraints that we consider. We give the technical details of our calculation in Section 3 and present the numerical analysis of the ISS one-loop corrections to λ HHH in Section 4. A short conclusion is given in Section 5. We present the details of the parameterization adopted for the light neutrino mass matrix in Appendix A and the analytical expressions of the one-loop corrections involving the neutrinos are collected in Appendix B.

Model and constraints
While our calculation and the analytical results presented in Section 3 are applicable to all models with extra fermionic gauge singlets and Majorana neutrinos like the type I seesaw [10][11][12][13][14][15][16] or the linear seesaw [22][23][24][25], we will focus in this work on the inverse seesaw (ISS) model for illustrative purposes. After introducing the model and the different parameterizations used to reproduce neutrino oscillations data, we will present the theoretical and experimental constraints considered in our study.

The inverse seesaw model
One particular variant of the type I seesaw is the ISS model [17][18][19] which has very interesting characteristics leading to a rich phenomenology. In the ISS model the suppression mechanism that guarantees the smallness of neutrino masses is the introduction of a slight breaking of lepton number in the singlet sector (composed of right-handed neutrinos ν R and new gauge singlets X with opposite lepton number), in the form of a small Majorana mass µ X for the X singlets, compared to the electroweak scale v ∼ 246 GeV. This allows for large Yukawa couplings compatible with a low (TeV or even lower) mass for the seesaw mediators, contrary to the seesaw model of type I for example, where the mediators have a mass of the order of the GUT scale or the Yukawa couplings are very small.
In the inverse seesaw, the additional terms to the SM Lagrangian are where Φ is the Higgs field and Φ = ıσ 2 Φ * , i, j = 1 . . . 3, Y ν and M R are complex matrices and µ X is a complex symmetric matrix whose norm is taken to be small since lepton number is assumed to be nearly conserved. In this work, we do not consider a possible Majorana mass term for the right-handed neutrinos ν R since this extra parameter is not relevant to our study. It would only induce negligible corrections to the heavy neutrino masses and the observable that we consider conserves lepton number. Assuming 3 pairs of ν R and X, the 9 × 9 neutrino mass matrix reads after electroweak symmetry breaking in the basis (ν C L , ν R , X), with the 3 × 3 Dirac mass matrix given by m D = Y ν Φ . M ISS being complex and symmetric, we can use the Takagi factorization to write where U ν is a 9 × 9 unitary matrix. A specificity of the ISS model is the presence of a nearly conserved lepton number. The light neutrino masses are then suppressed by the small lepton number breaking parameter µ X and the heavy Majorana neutrinos, which have nearly degenerate masses, form pseudo-Dirac pairs. This can clearly be seen if we consider only one generation. In the inverse seesaw limit µ X m D , M R , we have one light neutrino ν and two heavy neutrinos N 1 , N 2 with masses . (2.5) With three generations, M ISS can be diagonalized by block to give the light neutrino mass matrix, at leading order in the seesaw expansion parameter The next order terms are given in Appendix A. This 3 × 3 complex symmetric mass matrix is diagonalized by using a unitary matrix identified with the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U PMNS [35,36]: with m n 1 , m n 2 and m n 3 the masses of the three light neutrinos.
In order to reproduce low-energy neutrino data, different parameterizations can be introduced. Working in the basis where M R is diagonal with entries M i , neutrino oscillations are generated by off-diagonal terms in m D and µ X . In a first parameterization, we can reconstruct m D as a function of neutrino oscillation data and high energy parameters. This leads to a Casas-Ibarra parameterization [37] adapted to the inverse seesaw is a complex orthogonal matrix that can be expressed as with c i = cos θ i , s i = sin θ i , θ i being arbitrary complex angles. The other possibility is to use the µ X -parameterization that was introduced in ref. [38], Both parameterizations are based on eq.(2.6) where only the leading order term in the seesaw expansion is considered. While this is sufficient in most of the parameter space, these formulas fail to reproduce low-energy neutrino data when the active-sterile mixing becomes very large. Indeed, a large active-sterile mixing corresponds to a large seesaw expansion parameter m D M −1 R , which makes the next order terms presented in eq.(A.1) relevant. Including the next order terms in the seesaw expansion in the µ X -parameterization gives which allows to better reproduce neutrino oscillation data. The complete derivation of this formula is given in Appendix A. Finally, we need to specify the couplings between SM particles and the new fields that are relevant for our calculation of the corrections to the triple Higgs coupling λ HHH . Following ref. [21], we introduce the B and C matrices defined as where V L is the unitary matrix that diagonalizes the charged lepton mass matrix M charged according to with V R another unitary matrix. In the Feynman-'t Hooft gauge and in the mass basis, the relevant interaction terms in the Lagrangian are where g 2 is the SU(2) coupling constant, θ W is the weak mixing angle and P L , P R are respectively (1 − γ 5 )/2 and (1 + γ 5 )/2.

Constraints on the ISS model
Strong experimental and theoretical constraints on the parameter space of the model have to be considered, in particular on the size of the active-sterile mixing. Our use of the modified Casas-Ibarra or µ X -parameterization allows to reproduce neutrino oscillation data. In our numerical study, we explicitly check the agreement with the neutrino masses and mixing obtained in the global fit NuFIT 3.0 [39]. The light neutrino masses are also chosen to agree with the Planck result [40] 3 i=1 m n i < 0.23 eV . (2.17) The mixing between the active and sterile neutrinos will also induce deviations from unitarity in the 3 × 3 sub-matrixŨ PMNS of the full 9 × 9 mixing matrix U ν , that controls the mixing between the light neutrinos [41,42]. Using a polar decomposition, this square complex matrix can be expressed asŨ where η is a Hermitian matrix that encodes the deviations from unitarity. We have included the following constraints from a recent fit [43] to electroweak precision observables, tests of CKM unitarity and tests of lepton universality, In the presence of a large active-sterile mixing, the off-diagonal entries in the neutrino Yukawa couplings Y ν might also induce large branching ratios for lepton flavor violating (LFV) decays.
We have implemented the analytical expressions from ref. [21] for the LFV radiative decays and the LFV three-body decays. The corresponding experimental upper limits on the LFV radiative decays [44,45] are at 90% C.L. while the upper limits on LFV three-body decays [46,47] are at 90% C.L. We will also require in our study that Yukawa couplings are perturbative since the complex angles of the R matrix in the Casas-Ibarra parameterization or the use of Y ν as an input parameter in the µ X -parameterization can lead to arbitrarily large entries in Y ν . We will ensure the perturbativity of the Yukawa couplings by requiring for i, j = 1 . . . 3. Since the decay width of heavy neutrinos grows like m 3 n when m n m H , we also require that their decay width verifies, for i = 4 . . . 9, in order for the quantum state to be a definite particle. The formulae used to calculate the heavy neutrino widths are taken from ref. [48].

Framework of the calculation
Our calculation is done in the Feynman-'t Hooft gauge and we use the Lagrangian of eq.(2.16) for the neutrino interactions. The SM scalar potential is written as with the Higgs field Φ given by H stands for the Higgs boson, G 0 the neutral Goldstone boson, G ± the charged Goldstone bosons and v 246 GeV is the vacuum expectation value (vev) of the Higgs field. We can define the Higgs tadpole t H , the Higgs mass M H and the triple Higgs coupling λ HHH as follows, This helps to redefine the triple Higgs coupling using t H , M H and v as input parameters, At tree-level, t H = 0 and we recover the usual definition of the tree-level triple Higgs coupling, For the one-loop corrections to the triple Higgs coupling, our set of input parameters that need to be renormalized in the on-shell (OS) scheme will be the following: We use the following relations to define the Higgs vev v and the weak angle as well as e 2 = 4πα. We require that we have no tadpoles at one loop: H being the one-loop un-renormalized contributions to t H . For the other parameters we introduce their counter-terms as follows, The full renormalized one-loop triple Higgs coupling is finally and λ HHH (q * H ) stands for the un-renormalized one-loop contributions to the process H * → HH with the momentum q * H for the off-shell Higgs boson H * . For the numerical analysis carried in the next section, we define the deviation induced by the BSM contribution ∆ BSM as where λ 1r,SM HHH stands for the renormalized one-loop SM contribution without the light neutrinos.
Introducing the notation Σ XY for the self-energy of the process X → Y , we use the usual OS conditions for M W , M Z and M H , For the electric charge e we use the following condition to be independent from the light fermion masses [49,50], (3.14) For the Higgs field renormalization we have The neutrino interactions induce changes in the W and Z self-energies as well as in the Higgs tadpole, self-energy and self-couplings. We display in Fig. 1 the Feynman diagrams for the neutrino contributions to the W , Z and Higgs bosons self-energies, the Higgs tadpole and the one-loop un-renormalized triple Higgs coupling. We also collect in Appendix B the analytical expressions of the neutrino contributions to δM W , δM Z , δt H , Σ HH and λ (1) HHH . They were obtained using FeynArts 2.7 [51] and FormCalc 7.5 [52], in which we have implemented our own Model File for the ISS model. The scalar and tensor loop functions [53,54] have been evaluated with LoopTools 2.13 [52,55,56]. We have checked numerically that the UV divergences cancel in the final result and that the renormalized one-loop triple Higgs coupling does not depend on the choice of the renormalization scale.

Numerical results
We present in this section the phenomenological study of the one-loop corrected triple Higgs coupling and the dependence of the corrections induced by the heavy neutrinos on the relevant input parameters of the ISS model. The SM parameters are taken from the Particle Data Group (PDG) [57] (with the exception of the SM Higgs boson mass) and read as The up-, down-and strange-quark masses are also taken from the PDG, but their impact on the calculation is negligible so that we do not list them here. The lightest neutrino mass  to comply with cosmological constraints as stated in eq.(2.17). We have explicitly checked that choosing a smaller mass for n 1 does not qualitatively modify our results and would only induce negligible numerical corrections to our final conclusions. We chose the normal ordering for the neutrino masses and the light neutrino mixing parameters are taken from NuFIT 3.0 [39], with δ CP = 0. Since the contributions of the light neutrinos are negligible and flavor constraints do not play an important role in our final conclusion, we do not expect our conclusion to change if we consider the inverted ordering. In our study, we will focus on two choices for the off-shell Higgs momentum q * H = 500 GeV and q * H = 2500 GeV. These choices follow from the behavior of the BSM corrections that exhibit a similar dependence on q * H between the ISS model and the simplified Dirac 3+1 model that was studied in ref. [27]. In particular, the maximal negative deviation was obtained for q * H = 500 GeV while the maximal positive deviation was obtained for large off-shell Higgs momenta. To facilitate the comparison between the Majorana ISS case and the simplified Dirac case we take the same fixed values of q * H as in ref. [27] in all the scans.

Casas-Ibarra parameterization
In order to get an insight into the parameter space of the ISS model we perform a scan in a Casas-Ibarra parameterization, see eq.(2.8). The goal is to get an idea of the corrections that are obtained in this parameterization and the impact of the constraints on the model. We perform a random scan using a flat prior on the three real rotation angles θ 1/2/3 of the orthogonal matrix R and a logarithmic prior on both the lepton number violating term µ X so that the Majorana mass term is µ X ≡ µ X I 3 , and the mass term M R so that the matrix M R is We take all mass and rotation matrices to be real in order to avoid generating CP violation. We use 180 000 randomly generated points in the following parameter range, 7.00 × 10 −4 eV ≤ µ X ≤ 8.26 × 10 4 eV.
The range choice for the parameter µ X follows µ , see eq.(2.4). Heavy neutrino masses below 200 GeV are better probed with direct searches at colliders [42] (see also ref. [26] and references therein), thus we do not take M R < 0.2 TeV. The result of our scan is displayed in Fig. 2 (upper row) in the M R − µ X plane. The topright corner (in yellow) of the parameter space is excluded by theory constraints, essentially the perturbativity of the neutrino Yukawa couplings. The region in light blue is excluded by EWPO, while the region in purple is excluded both by EWPO and theory constraints. The dashed black line displays the limit coming from neutrino oscillations. This comes from a breakdown of the leading-order Casas-Ibarra parameterization when the active-sterile mixing is too large, as is evidenced by the flat behavior in M R . In variants of the type I seesaw including the inverse seesaw, the active-sterile mixing is proportional to the seesaw expansion parameter m D M −1 R . However, in the Casas-Ibarra parametrization, m D grows linearly with M R , see eqs.(2.8)-(2.9). As a consequence, m D M −1 R appears constant in M R but increases when µ X decreases according to The breakdown happens for µ X ∼ 3 eV, which in turns roughly corresponds to m D /M R ∼ 0.1 when taking m ν = m n 3 . It is worth noting that this value can be predicted from eq.(A.1), where next-order corrections to the light neutrino mass matrix appear at O(m 2 D /M 2 R ), and from the current error on ∆m 2 being at the percent level [39]. The most stringent experimental constraint comes from LFV observables as displayed by the solid vermilion line. The top-left corner (in green) is allowed by all the constraints.
This scan has to be compared to the map of ∆ BSM displayed in Fig. 2 (lower row), fixing the off-shell Higgs momentum at q * H = 2500 GeV. The parameter space passing all the constraints only contains corrections up to ∼ +1%. The most interesting regions are in vermilion, blue, yellow and purple where ∆ BSM reaches +15%, +25%, +35% and more than +35%, respectively. In order to enter these regions, it is needed to escape the LFV constraints as much as possible. Following ref. [38] we will investigate this region using the µ X -parameterization and start with the case of degenerate heavy neutrinos.
Higgs coupling we want to escape these constraints and we require for example (Y ν Y † ν ) 12 = 0 since decays that involve a µ−e transition usually give the strongest constraints. This leads to either a diagonal Yukawa matrix or a Yukawa texture as defined in ref. [38], with degenerate heavy neutrinos, M R ∝ I 3 .
We investigate in this sub-section the case of the degenerate heavy neutrinos in a µ X -parameterization with the texture Y (1) τ µ taken from ref. [38] and defined below, We display in Fig. 3 (left) the two-dimensional scan in the plane (M R , |Y ν |) where M R represents the common scaling factor of the 3 × 3 diagonal mass matrix M R . The off-shell Higgs momentum is again fixed at q * H = 2500 GeV. A large part of the parameter space is excluded and a maximum of ∆ BSM ∼ +5% can be reached at M R 13 TeV. When compared to the Casas-Ibarra scan, this is the expected order of magnitude for the correction when entering the vermilion region which is excluded by LFV observables only. For large M R the most important constraint is the neutrino width (2.31). For lower M R the constraints are driven by the violation of the unitarity of the 3 × 3 matrixŨ PMNS controlling the mixing between the light neutrinos.
To get an insight into the behavior of the contour lines in Fig. 3 (left) we display a onedimensional plot of the neutrino correction ∆ BSM at a given M R = 10 TeV, as a function of the Yukawa scaling factor |Y ν |, in Fig. 3 (right). The correction is negligible for low Yukawa scaling factors, then rises to a maximum at ∼ +60% at |Y ν | 2.5 before dropping rapidly and eventually becoming negative for large Yukawa scaling factors.
From this behavior we devise the following approximate formula to reproduce ∆ BSM at M R > 3 TeV, The numerical coefficients are found to be universal in term of the parameters of the model and only depend on the kinematics of the off-shell Higgs boson, for the case of the three textures of ref. [38] as well as for the case of a diagonal texture. The dependence of the numerical coefficients on the kinematics of the off-shell Higgs boson is expected, as when compared to the full calculation they would result from the loop functions depending on q * H , see Appendix B. It is expected that eq.(4.6) be valid for the whole class of textures introduced in ref. [38]. At a given M R > 3 TeV, the approximate formula in eq.(4.6) is driven at low |Y ν | by the positive contribution and by the negative contribution at high |Y ν |, the latter falling more rapidly than the positive increase at low |Y ν |. This reproduces the behavior seen in Fig. 3 (right) where the result of the fit is also displayed. We can also reproduce the contour  lines for high M R in Fig. 3 (left) as seen from the green contour lines coming from the fit, that agree to a very good extent with the full contour lines for M R > 3 TeV.
The approximate formula in eq.(4.6) implies that the best way to maximize the neutrino effects on the triple Higgs coupling would be to maximize the ratio . The Yukawa couplings being real and limited by perturbativity requirements, this leads to the choice of a diagonal texture, Y ν ∝ I 3 . This will be considered in the next sub-section, but with the condition of degenerate heavy neutrinos being relaxed. In such a way the constraints on the non-unitarity of the matrixŨ PMNS are softened and the blue region of the Casas-Ibarra scan of Fig. 2, excluded by EWPO as well as by LFV observables, moves down.

Hierarchical heavy neutrinos
The analysis carried in the previous sub-section has lead us to consider a diagonal Yukawa matrix, Y ν = |Y ν |I 3 . In order to reduce as much as possible the impact of unitarity constraints on η for the matrixŨ PMNS , we chose hierarchical heavy neutrinos with M R = diag(M R 1 , M R 2 , M R 3 ) and we still work in the µ X -parameterization. More specifically, for illustrative purpose within this class of parameters, we chose with M R being a rescaling factor that is varied between 200 GeV and 20 TeV. This ensures that all the diagonal constraints of eq.(2.19) have the same impact on our study. This specific choice maximizes the individual contribution of each heavy neutrino, which in turns will maximize the one-loop correction to the triple Higgs coupling originating from the leptonic sector. Other choices for the heavy neutrino masses will only reduce the allowed maximum value of the triple Higgs coupling deviation from the SM. The result of the parameter scan in the M R − |Y ν | plane is displayed in Fig. 4. On the left-hand side, we display the map of ∆ BSM for an off-shell Higgs momentum q * H = 500 GeV. As already expected by the analysis in the simplified model of ref. [27], the heavy neutrino corrections are negative, and they reach a minimum of ∼ −8%, close to the minimum that was obtained in the simplified model. The most interesting results are displayed in the right-hand side of Fig. 4, for q * H = 2500 GeV. The corrections can now reach a maximum of ∼ +30%, similar to what has been obtained in the case of a simplified model. The corrections are generically bigger in the ISS model than in the simplified model, but the constraints are also stronger, reducing the heavy neutrino corrections back to the maximum obtained in the simplified model. This also confirms in a realistic, renormalizable, low-scale seesaw model that heavy Majorana neutrinos can induce sizable deviation to the triple Higgs coupling.
As a further test of our approximate formula for the heavy neutrino corrections, the green lines in Fig. 4 are the approximate contour lines obtained using eq.(4.6) but rescaled with a common factor γ = 0.51. This type of rescaling was expected as now the heavy neutrino mass matrix is not proportional to the identity matrix anymore. Once again we obtain a very good approximation for M R > 3 TeV, and in particular in the region allowed by the constraints. This approximate formula thus describes well the behavior of ∆ BSM in the allowed region of the parameter space, for q * H = 2500 GeV. We end this section with a comparison with the currently expected sensitivity to the triple Higgs coupling at the HL-LHC and at the future planned colliders. The sensitivities to the SM triple Higgs coupling are defined by its measure extracted from the Higgs pair production yields. As stated for example in refs. [58,59], a precision of ∼ 50% on the total cross section leads to a precision of ∼ 50% on the SM triple Higgs coupling. The sensitivity for the HL-LHC follows from ref. [32] (see also ref. [60]), scaled by a factor of 1/ √ 2 to account for both ATLAS and CMS accumulated data, while the sensitivity for the future colliders follow from refs. [33,34]. For the FCC-hh we do the same as for the HL-LHC to account for both ATLAS and CMS accumulated data 1 , as well as for the fact that the analysis in ref. [34] is only done for one search channel; we expect the sensitivity to improve when more search channels are taken into account. We display in Fig. 5 the maximally allowed deviation ∆ BSM (in percent), in black solid line, as a function of the heavy neutrino rescaling factor M R (in TeV). This is compared to the sensitivities to the SM prediction for λ HHH in the case of the More specifically and using current experimental constraints, the ILC at a center-of-mass energy √ s = 500 GeV could probe heavy neutrino masses in the range 8.5 < M R < 10.5 TeV, at 1 TeV with 5 ab −1 of data this extends to the range 5 < M R < 17.5 TeV. The FCC-hh collider could extend the analysis to a bigger range 3.3 < M R < 20 TeV. Indirect searches and in particular EWPO could probe heavy neutrinos with masses in the multi-TeV range and future improvements are expected, especially at future e + e − colliders [63]. Improved constraints on the EWPO would tend to shift the left-hand part of the black curve in Fig. 5 towards the right. This makes the triple Higgs coupling a new, viable and attractive observable to test low-scale seesaw mechanisms that will be complementary to improved EWPO measurements.

Conclusions
We have investigated in this article the one-loop effects of heavy neutrinos on the triple Higgs coupling in the framework of an inverse seesaw model, that is a realistic, renormalizable model accounting for the masses and mixings of the light neutrinos. After having presented the model and its constraints, both theoretical and experimental, in Section 2, we have given the technical details of the one-loop calculation in Section 3. We have presented in Section 4 our numerical investigation of the model. After having performed a scan in a Casas-Ibarra parameterization of the neutrino input parameters we have found that a µ Xparameterization is more suitable to get the maximal effects on the triple Higgs coupling and we have obtained a deviation as high as ∼ +30% for the class of parameters in which the 3 × 3 heavy neutrino mass matrix M R is diagonal and hierarchical while the 3 × 3 neutrino Yukawa texture is proportional to the identity matrix. This confirms our expectations coming from the simplified model analysis, and establishes the triple Higgs coupling as a viable, new observable to probe heavy neutrino mass regimes that are hard to probe otherwise, as this deviation is at the current limit for the expected sensitivity at the HL-LHC but clearly visible at the ILC and at the FCC-hh. Heavy neutrinos can also give rise to new diagrams that contribute to the complete HH production cross section and need to be evaluated. We leave this for future projects.
A Next order corrections in the seesaw expansion parameter to the µ Xparameterization Following the method of ref. [66], we can diagonalize M ISS by block to an arbitrary order in the seesaw expansion parameter m D M −1 R . This gives for the 3 × 3 light neutrino mass matrix in agreement with previous results [67]. This can be written in a symmetric form If m D is invertible, we can then express µ X as a function of M light and the other blocks of M ISS , The light neutrino mass matrix is diagonalized by using the unitary PMNS matrix. Using eq.(2.7) to rewrite M light in eq.(A.3), we get a formula for the µ X -parameterization that includes the effect of sub-leading terms in the seesaw expansion, It is easy to see that if we were to consider only the leading order term in the seesaw expansion, we would recover eq.(45) from ref. [38]. Interestingly, our results would not be modified by the addition of an extra mass term µ R ν C R ν R . The neutrino mass matrix would then be where taking ||µ R || ||m D ||, ||M R || corresponds to the inverse seesaw limit while taking ||µ R || ≥ ||M R || leads to the extended seesaw limit. In both cases, the next order corrections to M light are given by eq.(A.1) in the limit where ||µ X M −1 R µ R || ||M R ||. Thus, eq.(2.12) would remain unchanged 2 .

B Analytic expressions of the new ISS contributions
We give in this appendix all the analytic formulae of the new ISS contributions involved in the calculation of the renormalized one-loop triple Higgs coupling presented in Section 3. The SM contributions, denoted with a SM, can be found in ref. [49] and will not be reproduced in this appendix.

B.1 Counter-terms
By convention all loop integrals in this sub-section are to be understood as their real part only. We use the conventions of LoopTools 2.13 [52,55,56] for the scalar integrals and the tensor coefficients.