Constraints on coloured scalars from global fits

We consider a simple extension of the electroweak theory, incorporating one $SU(2)_L$ doublet of colour-octet scalars with Yukawa couplings satisfying the principle of minimal flavour violation. Using the HEPfit package, we perform a global fit to the available data, including all relevant theoretical constraints, and extract the current bounds on the model parameters. Coloured scalars with masses below 1.05 TeV are already excluded, provided they are not fermiophobic. The mass splittings among the different (charged and CP-even and CP-odd neutral) scalars are restricted to be smaller than 20 GeV. Moreover, for scalar masses smaller than 1.5 TeV, the Yukawa coupling of the coloured scalar multiplet to the top quark cannot exceed the one of the SM Higgs doublet by more than 80%. These conclusions are quite generic and apply in more general frameworks (without fine tunings). The theoretical requirements of perturbative unitarity and vacuum stability enforce relevant constraints on the quartic scalar potential parameters that are not yet experimentally tested.


Introduction
The main focus of the first LHC run was on the search for the Standard Model (SM) Higgs boson and culminated in its spectacular discovery in 2012 [1,2]. After the conclusion of the run 2, the properties of this new scalar are much better known, and they seem to agree very well with the SM predictions. In combination with direct LHC searches for additional spin-0 particles, they put strong bounds on extensions of the scalar sector of the SM, as for example the Two-Higgs-Doublet model [3]. On the other hand, new physics beyond the Standard Model (BSM) is needed to fix the shortcomings of the current theoretical framework, and many BSM scenarios could manifest themselves through the appearance of additional scalars at the LHC. Most of these SM extensions contain scalar fields which transform as singlets under the SU (3) C group of QCD and, therefore, do not carry colour charge. But coloured scalars are less constrained than one might naively think: they do not mix with the SM Higgs and contribute to its signal strengths only at the loop level. However, recent LHC searches for scalar resonances are expected to yield strong constraints on their properties.
Here we want to study a simple extension of the SM scalar sector with one SU (2) L doublet of scalar fields that transforms as an octet under the colour SU (3) C group. In some sense, one can think of it as a Two-Higgs-Doublet model in which one of the doublets is coloured. Scalars of this type with masses around the electroweak scale could emerge naturally from SU (5) or SO(10) unification theories at high energy scales [4][5][6][7]. As this model was first proposed by Manohar and Wise [8], we will refer to it as the Manohar-Wise (MW) model in the following. Several phenomenological aspects of this dynamical scenario have been studied in the literature , yet a comprehensive analysis including all relevant up-to-date constraints has not been performed.
In this article we will present such an analysis. After defining the MW model in Section 2, we briefly describe in Section 3 the HEPfit package [30] that we employ to perform the global fits. This section explains in detail the various constraints on the parameter space of the model, which come from theoretical considerations (Section 3.1), Higgs measurements (Section 3.2), direct LHC searches (Section 3.3) and flavour observables (Section 3.4). Putting all pieces together, we show the results of our fits in Section 4. Finally, Section 5 contains a summary of our findings and a brief outlook on future prospects.

The scalar colour-octet model
The MW model extends the SM with one electroweak doublet of colour-octet scalar fields with hypercharge Y = 1 2 . Since it has colour charge, the new scalar multiplet does not mix with the SM Higgs doublet. Furthermore, the coloured particles cannot acquire a vacuum expectation value (vev) because colour must be conserved. Therefore, only the SM Higgs boson will acquire a vev which will minimise the most general potential that can be built with this scalar sector: where Φ = (φ + , φ 0 ) T is the usual SM doublet, the traces are taken in colour space, and i and j denote SU (2) L indices. The additional (8, 2) 1/2 scalar fields S A = (S A,+ , S A,0 ) T are contained in the multiplet S = S A T A with T A the generators of the SU (3) C group. All potential parameters are real except ν 3 , ν 4 and ν 5 , but performing a phase rotation we can always take ν 3 to be also real. From Eq. (2.1) we can derive the masses of the physical neutral octet scalar (m R ), the neutral octet pseudoscalar (m I ) and the charged octet (m S ± ), which are split by the Higgs vev, φ 0 = v/ √ 2: The kinetic term (the factor two gives the correct canonical normalisation for the fields) generates the interaction of the octet scalars with the gauge bosons through the covariant derivative with G µ = G A µ T A the octet gluon field. Thus, these interactions are determined by the gauge symmetry and do not introduce additional free parameters.
The coloured scalars can also couple to the quarks through the Yukawa interaction. In order to guarantee the suppression of unwanted flavour-changing neutral currents, which are extremely suppressed experimentally, we will assume the principle of Minimal Flavour Violation (MFV) [31,32], which is based on the hypothesis that all Yukawa matrices are proportional to the same flavour structures that break the SU (3) Q L ⊗ SU (3) u R ⊗ SU (3) d R symmetry in the SM. This is in fact one of the main motivations of the MW model because the only scalar representations that can couple to quarks and be compatible with this principle are the colour octet or singlet electroweak doublets [8]. With this assumption, the Yukawa couplings of the coloured scalars take the form Here, i and j are family indices, Y q = √ 2M q /v (q = u, d) denote the up and down SM Yukawa matrices and the tilde in the S field indicates charge conjugation. The proportionality constants η U and η D are, in general, complex parameters.
Looking at Eqs. (2.1) and (2.5), we observe that the MW model contains 18 more parameters than the SM, 14 of which are real while the other 4 are imaginary phases. In order to simplify the phenomenological analysis and to reduce the total number of free parameters we will only work in the CP-conserving limit. This assumption removes the imaginary parts of ν 4 , ν 5 , η U and η D and we end up with only 14 new free parameters.

Fit constraints
Our statistical data analysis will be based on a global Bayesian fit. We make use of the public HEPfit package [30], which is interfaced with the Bayesian Analysis Toolkit [33]. This code has been already applied to several BSM analyses, including the Two-Higgs-Doublet model [3,[34][35][36] and the Georgi-Machacek model [37]. We have adapted the code, including the additional routines needed to study the MW model. These routines are also public and can be extended in future works to incorporate additional observables. In our fit we have only included observables that have been directly calculated with this model. Constraints obtained from existing bounds on higher-dimension operators 1 are omitted, since we prefer to directly include in the fit the observables used to derive those constraints. Indeed, HEPfit has also been proven to be extremely reliable to constrain higher-dimension operators [39][40][41]. One of the key features of HEPfit is its independence of other codes at runtime, which provides a very fast framework for statistical analyses.
Since we will try to use in our fit all the available information, we did not have any previous constraints on the MW parameters. Therefore we decided to use a uniform distribution as a prior for the 14 free parameters of the MW model. The ranges adopted are shown in Table 1. The range taken for η U follows from the assumption of a perturbative top Yukawa coupling, while larger values of η D are possible since m b m t . The dependence of our results on the priors used turns out to be small, as long as these priors are reasonable. For instance, increasing the range of m 2 S up to 2 2 TeV 2 leads to the same limits, but we do not gain any further information since many of the direct-searches analyses included do not go beyond 1.5 TeV. For ν n , µ n and η U the posterior probabilities are basically zero in the limits of the ranges considered and the same constraints are always found if this (reasonable) condition is satisfied. The only exception is η D which could not be constrained within the range considered. However, higher values of η D will always bring stronger limits on the other parameters. The parameter space is constrained by the theoretical and experimental requirements that we now discuss. Parameters (0.4 2 , 1.5 2 ) TeV 2 (-10, 10) (-10, 10) (-5, 5) (-20, 20) Table 1. Prior values of the fitted parameters.

Theoretical constraints
Since we are interested in the LHC phenomenology, we will assume that the MW model is well defined up to the few TeV scale. As a consequence, the renormalisation group (RG) evolution should be stable between the electroweak scale and 5 TeV. By RG stability we mean the absence of Landau poles as well as the fact that the scalar potential V MW should be bounded from below at any of the considered scales. The RG equations (without Yukawa and gauge couplings) were taken from Ref. [20] and a set of necessary positivity bounds was derived in Ref. [25]. Another strong constraint on the quartic couplings of V MW is the unitarity requirement that the two-to-two scattering processes of the scalar particles should not have a probability larger than 1. This condition is usually expressed in terms of the partial wave amplitudes a j . Using their LO (a (0) j ) and NLO (a (1) j ) expressions [25], perturbative unitarity can be imposed at LO, NLO and at the so-called NLO+ approximation [42] that includes the square of the NLO correction and therefore contains some, but not all, NNLO terms [25]: Note also that the LO and NLO expressions for the partial wave amplitudes are only available in the large-s approximation and should not be applied below a certain scale µ u , which we choose to be µ u = 1.5 TeV as the current limits on the mass of the colour-octet scalars are around 1 TeV [28]. On top of these unitarity bounds, we also impose perturbative behaviour of the quantum corrections and allow only for scenarios in which the NLO contributions to the partial wave amplitudes are smaller in magnitude than the LO term.
In the fits, we run up to 3 or 5 TeV each parameter set, sampled at the electroweak scale, and control at each intermediate step if the parameters still comply with positivity and (above µ u ) perturbative unitarity. If this is not the case, the running is stopped and the corresponding cut-off scale is returned as an output.

Higgs constraints
In the MW model there are additional contributions to some of the production and decay channels of the SM Higgs boson, due to the presence of the colour-octet scalars. At the LHC, the SM Higgs boson is produced mainly through gluon fusion (ggF), vector boson fusion (VBF), associated production with a tt pair (tth) and associated production with vector bosons (Wh/Zh). The LO contribution to the first process is a one-loop amplitude already in the SM. As the coloured scalars couple to gluons and to the SM Higgs boson, they will contribute to the gluon-fusion production mode at the same order than the SM. Similarly, the decay of the SM Higgs boson to photons is also a one-loop process in the SM to which the new scalars contribute. The LHC data for the Higgs physics are parametrized in terms of the Higgs signal strengths, which are defined as the measured cross section times branching ratio for a given production and decay Higgs channel, in units of the SM prediction. Table 2 compiles the experimental papers from which we have taken the values of the different Higgs signal strengths. Comparing the measured values of these observables with their theoretical predictions [22], we will try to constrain the parameters ν 1 , ν 2 and ν 3 , which are directly related with the mass splittings of the new scalars in Eq. (2.2).

Constraints from direct searches
In order to find constraints from direct experimental searches of additional scalars, we will compare the theoretical prediction of the production cross-section times branching ratio of a given process, σ · B, with the experimental upper limits of the ATLAS and CMS collaborations. The channels that we have included in the fit are shown in Table 3. Since decays of the colour-octet scalars into purely electroweak bosons are forbidden, we have only considered their decays into gluons and quarks. Furthermore, in order to be sure that the decay modes studied are indeed the dominant ones, we impose for each search that the mass difference between the coloured scalar we are analysing and the lightest member of the S multiplet is smaller than the W ± mass [28]. In this way, we forbid kinematically the decays into another coloured scalar and a gauge boson, which can become the dominant decay modes for some configurations of the parameter space. As we will show later, this assumption is well justified, given the constraints that perturbative unitarity enforce on the mass splittings of the new scalars.  Table 3. Data from direct searches included in the fit. CMS8 means CMS at 8 TeV, the rest is at 13 TeV.
Since HEPfit cannot generate events by itself, the theoretical predictions for σ · B must be provided with other tools. We have used MadGraph [74] to create tables with the values of σ · B for different input choices of the parameters on which these observables depend: η U , η D , (ν 4 + ν 5 ) and m S . These tables are read by HEPfit, which performs a linear interpolation to obtain the value of the observables at any point. The error of the linear interpolation is estimated to be 10% or less for log(σ · B).
The experimental data are also provided in the form of tables, which compile the values of the 95% upper limits on σ·B, as a function of the resonance mass. In these tables, HEPfit also performs a linear interpolation if needed.
In order to compare the theoretical results with the experimental data, we define the ratio of the theoretical prediction over the experimental upper limit. To this ratio we will assign a Gaussian distribution (restricted to positive values) centered at 0 such that the value 1 is excluded with a 95% probability.

Flavour constraints
Regarding flavour transitions, we have included in the analysis the B 0 s −B 0 s mixing and the decay B 0 s → µ + µ − , which are the most constraining observables. The expressions for the relevant Wilson coefficients were taken from Ref. [21], where a complete one-loop calculation was performed. Only the charged scalars contribute to these processes at the one-loop level, and the corresponding amplitudes involve the parameters η U and η D . The leading η U contributions are proportional to the top-quark mass, while the η D terms are weighted by the bottom-quark mass. Therefore, the transition amplitudes are mainly governed by the parameter η U .
These observables strongly depend on the numerical values of the relevant quark mixings. We cannot use the standard determinations of the CKM entries, since they assume the validity of the SM. In order to obtain non-contaminated values, we perform a specific CKM fit using HEPfit, which only includes observables that are not affected by the coloured scalars. The inputs of this CKM fit, shown in Table 4, are taken from the PDG 2020 [75] except |V cs |, for which we have adopted the more recent result of Ref. [76]. The value of |V ud | has been considerably shifted with respect to the one quoted in the PDG 2018 [77]. Combined with the reduction of its uncertainty, this generates a violation of unitarity in the first row of the CKM matrix at the 3σ level. Since there is no full consensus in the community, we have adopted a conservative attitude, increasing the error of |V ud | by a factor 2.4 so that CKM unitarity is recovered at 2σ, without modifying any central values. The results of this CKM fit, shown in Table 5, are then used as inputs in our analysis of the MW parameter space.  Table 6. Results for the fit of the oblique parameters S, T and U , excluding the information from R b .

Electroweak precision data
In this work, we are studying the most general CP-conserving MW model, so we have not imposed any additional assumptions such as Custodial Symmetry. Therefore, the coloured scalars generate non-vanishing contributions to the oblique parameters [78,79] that can be used to constrain the parameter space of the model. The standard bounds on S, T and U [75] include in the global fit the ratio R b = Γ(Z → bb)/Γ(Z → hadrons), which incorporates quantum corrections to Γ(Z → bb) that are enhanced by the large value of the top quark mass [80,81]. Since the ratio R b receives also sizeable corrections from the additional scalars, we cannot use these bounds. Similarly to what has been done before for the CKM matrix, we first determine the oblique parameters, performing a combined fit with HEPfit of electroweak precision observables without R b [3], and then use the resulting values as inputs in our analysis of the MW model. The values of the oblique parameters that we have obtained from this fit are shown in Table 6.
The theoretical predictions for S, T and U in the MW model have been taken from Ref. [11]. The oblique parameters depend on the masses of the new scalars and, therefore, provide important constraints on m S , ν 1 , ν 2 and ν 3 . Once, the oblique parameters have been fixed, we also include R b in the fit to the MW parameter space. The theoretical prediction of this observable, including the QCD corrections, has been taken from Ref. [12]. The ratio R b will basically constrain η U because, like for the flavour observables, the new physics contribution comes from one-loop diagrams involving virtual charged scalars, coupled to heavy quarks.

Theoretical constraints
As mentioned before, we have imposed that this model should satisfy perturbative unitarity and be well defined up to, at least, 3 TeV. These conditions generate bounds on all the parameters of the scalar potential that, in some cases, can be stronger than the ones obtained from experimental data. For instance, the quartic couplings are not constrained experimentally but the theoretical requirements imply that none of them can be, in modulus, higher than 6. The theoretical constraints can be translated to physical observables like the mass differences, which depend on these parameters. In Fig. 1 we show the theoretical constraints on the scalar mass differences for two choices of the UV scale, 3 TeV and 5 TeV. Obviously, requiring the model to be valid up to a higher scale results in stronger bounds. The figure exhibits also the importance of perturbative unitarity, showing how the limits get weakened if the scale above which we impose it is chosen to be higher. The large impact of these theoretical constraints can be better appreciated in Fig. 8, which displays the bounds they impose on the µ n and ν n potential parameters, in the form of two-dimensional correlated plots, compared with the final results from the global fit, including all constraints, that we will later discuss.

Experimental constraints
As we have shown before, using theoretical arguments we can constrain the parameter space of all the new quartic couplings that appear in the scalar potential. In this section we will use experimental measurements in order to constrain the masses of the coloured scalars and their Yukawa couplings. But first of all let us consider the observables that constrain also the parameters of the potential and let us compare the results with the ones obtained from theoretical requirements.
In Fig. 2 we show the constraints on ν 1 , ν 2 and ν 3 emerging from the Higgs signal strengths and compare them with the theoretically allowed regions, imposing RG stability up to 3 TeV. Clearly the theoretical restrictions are much stronger than the Higgs constraints so, a priori, one could think that the latter will be irrelevant in the global fit. However, the addition of these experimental observables will have some effect on the global fit, as we will show later. The fact that the constraints from the Higgs signal strengths, alone, are so weak is an expected behaviour, given that the coloured scalars only start to contribute to them at the one-loop level.
The oblique parameters constrain the scalar mass splittings, bounding the combinations of quartic couplings ν 2 ± 2ν 3 that appear in Eq. (2.2). This is clearly seen in the narrow allowed regions (orange) displayed in Fig. 3. At large values of ν 2,3 , the oblique parameters impose that ν 2 ≈ ±2ν 3 in order to reduce the splitting between the charged scalar and either the CP-odd or CP-even neutral scalars. The oblique bounds are relevant even when we consider the theoretical constraints up to 5 TeV, but they are specially important when we only require RG stability up to 3 TeV. Indeed, combining both constraints we can obtain a harder restriction on these parameters for the 3 TeV case.
The absolute mass scale m S is constrained by the electroweak ratio R b , which is also  Figure 3. Allowed two-dimensional regions of ν 2 and ν 3 from the oblique parameters (orange, 95% probability) and from theoretical constraints (green) with the NLO+ approximation. The darker and lighter green areas correspond to requiring RG stability up to 3 and 5 TeV, respectively. sensitive to the Yukawa alignment parameter η U . The m S -η U plane is also constrained by the two flavour observables that we have considered: B 0 s −B 0 s mixing and Br(B s → µ + µ − ). The left panel of Fig. 4 displays the regions allowed by each of these observables, individually, at 95% probability. The mass difference between B 0 s andB 0 s turns out to provide the strongest bounds, although they are quite similar to the ones from Br(B s → µ + µ − ). Furthermore, since the three observables exhibit a similar dependence on these parameters, the combined constraints obtained with the three observables together are not much better than the ones given by just ∆M Bs or Br(B s → µ + µ − ).
The m S -η U plane is also constrained by the direct searches, as shown in the right panel of Fig. 4. In this figure we have only included searches with at least one top quark in the final state because those are the channels generating the most interesting constraints. Since the figure does not include any observable constraining the scalar mass differences, we have set these mass splittings to zero. Doing so, we are able to totally exclude values of m S (equal in this case to the physical scalar masses) smaller than 1.1 TeV, for any value of the other parameters, with a probability of 95%. The plot shows also the additional constraints obtained on this plane from flavour observables.
In Fig. 5 we display the analogous constraints emerging from direct searches in channels without top quarks in the final state. Relevant bounds are only obtained at small values of η U or for high values of η D . This is to be expected because those are the parameter regions where the scalar decays into top quarks are suppressed and, therefore, the branching ratios for the other channels are increased. The decay channels to top quarks are of course the dominant ones outside these particular regions of the parameter space, and decay modes Similar features are found in the m S -η D plane, where the bound m S > 1.1 TeV is also obtained with a 95% probability from direct searches in channels with top-quark production. This could be, in principle, a bit surprising because one could naively expect that for very small values of η U those searches should not generate any constraint. This is true for the channels in which a neutral scalar decays to a tt pair, but not for a charged scalar decaying to tb, a process which also depends on η D . Indeed, in the right panel of Fig. 4, one observes that the lower bound on m S decreases when η U approaches zero. This is because all channels with neutral scalars become irrelevant in that limit, but the information from charged-scalar channels is still good enough to generate a quite strong constraint. Therefore, as long as η U and η D are not both extremely close to zero, we obtain good constraints on m S . Since the MW model is motivated by MFV, the particular region where the two quark Yukawa couplings are both zero does not seem to have much theoretical interest.

All constraints
Once we have analysed the constraints emerging from individual observables, we can combine all of them into a global fit. The first interesting result is that we are able to find 95% DS no t-quark 99.7% DS no t-quark 4 2 lower limits for the physical masses of our scalars and constraints on their mass splittings, as shown in Fig. 6. Note that the mass differences between the different scalars are now restricted to be smaller than 25 GeV with a 99.7% probability, if we impose RG stability and perturbative unitarity up to 3 TeV, and to be smaller than 20 GeV, with the same probability, when we go up to 5 TeV. This justifies a posteriori our approximation of not including in the analysis of direct searches the decays of the coloured scalars into another scalar plus a weak boson.
Another interesting result is shown in Fig. 7, where we plot again the constraints on the m S -η U and m S -η D planes. We can clearly see how η D is not constrained. On the other hand, |η U | cannot take values higher than 1.8, for masses of the scalars smaller than 1.5 TeV, within a probability of 95%. The constraint on the scalar mass scale is mainly coming from direct searches including decays to top quarks. This is why for values of η U close to zero the lower bounds on m S are weaker. However, as mentioned before, even for very small values of |η U | the constraints remain strong. We have also performed some global fits in which we set η U = 0 or η D = 0. For the first case we found that the limits on the scalar masses are roughly 100 GeV smaller, while for the second we found the same result than for the global fit in which we vary both parameters. Therefore, even in the hypothetical case that some symmetry would force one of the Yukawa couplings to vanish, one would still find important constraints on the scalar masses.
The two-dimensional correlations among the quartic potential parameters are displayed in Fig. 8, which compares the limits emerging from the global fit with those obtained with only theoretical constraints, imposing RG stability and perturbative unitarity up to 3 and 5 TeV. The allowed regions from theoretical constraints must be understood as having 100% Theo 3TeV Theo 5TeV 95% All (3TeV) 99.7% All (3TeV) 95% All (5TeV) 99.7% All (5TeV) m SI probability because they correspond to a discrete condition: a given point of the parameter space is either allowed or not. This introduces a small difference when comparing the theoretical and global limits, since the latter ones refer to a slightly smaller probability. Taking this into account, we can still see that ν 1 and ν 2 are slightly more constrained in the global fit, if RG stability is only imposed up to 3 TeV. Nevertheless, for the other observables the allowed regions remain almost invariant when including the experimental information. This is not a surprise because the measured observables depend on ν 1 , ν 2 and ν 3 . The improvements introduced by the global fit are a consequence of adding the oblique parameters and, specially, the Higgs signal strengths. In order to check this, we performed a fit with only theory and the Higgs signal strengths and the results for ν 1 and ν 2 were basically identical to the ones of the global fit. Indeed, although the Higgs signal strengths alone do not produce strong limits, combining them with the theoretical constraints, which reduce the allowed range of ν 3 , results in very good bounds on ν 1 and ν 2 .
If RG stability is imposed up to 5 TeV, the current experimental measurements have 95% All (3TeV) 99.7% All (3TeV) 2 1 0 1 2 U 800 1000 1200 1400 m S (GeV) 20 10 0 10 20 D Figure 7. Allowed regions on the m S -η U (left) and m S -η D (right) planes from the global fit, at 95% and 99.7% probability, imposing RG stability up to 3 TeV and perturbative unitarity at NLO+. Quite similar results are obtained, requiring RG stability up to only 1 TeV. a quite small impact in the allowed regions in Fig. 8, which are mainly governed by the much stronger theoretical constraints.
Finally, in Table 7 we present the marginalised allowed ranges for the parameters of the potential, from the global fit, with a probability of 95%. As can be seen there, none of the quartic couplings can be, in modulus, higher than 5 within this probability.

Summary
In this work, we have presented the first global fit of the MW model. We have combined the more relevant theoretical and experimental constraints, including flavour and electroweak precision observables, Higgs signal strengths and direct searches. The theoretical constraints have a quite large impact on the scalar parameter space, specially when RG stability is imposed up to 5 TeV, providing bounds on all the scalar potential parameters. The current experimental information on the Higgs signal strengths and the oblique parameters gives an improved sensitivity to ν 1 , ν 2 and ν 3 , which becomes relatively more important if RG stability is only imposed up to 3 TeV. This allows for more stringent constraints on the scalar mass splittings.
The flavour observables and R b constrain the m S -η U plane, but all of them in the same direction. The strongest limits come from ∆M Bs and Br(B s → µ + µ − ), which combined, Theo 3TeV 99.7% All (3TeV) Theo 5TeV 99.7% All (5TeV) require that |η U | < 1.8 (5.1) for scalar masses smaller than 1.5 TeV, within a probability of 95%. A quite strong bound on the absolute mass scale m S emerges from the LHC data on direct searches for new scalars. The more sensitive channels, which involve the production of top quarks in the final state, imply that the masses of all coloured scalars must satisfy the bound m S ± , m S R , m S I > 1.05 TeV (5.2) for any value of the other parameters, with a 95% probability. The global fit also restricts the scalar mass splittings to be smaller than 20 GeV, with a 95% probability, even when RG stability is only imposed up to 3 TeV. As shown in Ref. [28], even for tiny values of |η U | > 10 −7 , one still finds a strong lower bound on the scalar masses, provided |η D | > 10 −5 . These bounds can only be avoided for  Table 7. Allowed ranges for the quartic potential parameters from the global fit, at 95% probability.
fermiophobic scalars with η U ≈ η D ≈ 0, such that their decay into a fermion-antifermion pair is highly suppressed. In that case, they would have a completely different experimental signature, since they would either decay into a lighter coloured scalar and a gauge boson or would behave as strongly-interacting long-lived particles. Although fermiophobic scalars are not compelling from the MFV point of view of the MW model, their phenomenology could be interesting by itself, but requires a more specific analysis that we postpone to future works. The forthcoming Run3 of the LHC and its subsequent high-luminosity phase will provide much larger data samples, substantially increasing the sensitivity to coloured scalar particles. Mass scales up to 1.3 TeV seem to be reachable with the expected luminosity of 3 ab −1 . The constraints on those potential parameters not related to scalar masses will remain, however, largely dependent on theoretical considerations, unless a real discovery of a colourful scalar state emerges from the data.