Road map through the desert with scalars

In the context of the gauge coupling unification, we present a comprehensive analysis of the extensions of the Standard Model with vector-like fermions and scalars. We find 145 models that satisfy the unification condition, which are distinguishable by the number of new particles in the spectrum and by their transformation properties under the gauge symmetry group of the Standard Model. For all models we derive lower bounds on the exotic fermion and scalar masses, stemming from the measurement of the strong gauge coupling scale dependence, from the heavy stable charged particle searches, and from the electroweak precision tests. We also discuss the potential of testing the unification scenarios at the future 100 TeV collider and in the proton decay experiments. We show that many models can already be excluded based on the current data, while many others will be entirely probed in the coming years.


Introduction
The idea of unification, i.e. the common high-scale origin of three fundamental interactions of the Standard Model (SM), has been driving intellectual efforts in particle physics since the mid 1970s [1][2][3][4][5]. The well known observation that the strength of the renormalized gauge couplings of the SM, which significantly differ at the scale of the electroweak (EW) interactions, seems to become similar at very high energies, may suggest the existence of a new unified description of fundamental forces known as a Grand Unified Theory (GUT).
It has also been known for decades that in the standalone SM the precise unification of the gauge couplings is not really achieved, as their GUT-scale values deviate from one another by several percent. This fact provided a strong motivation to consider various extensions of the SM in which the unification is guaranteed by the appropriate choice of the particle spectrum .
In [30] we proposed to use the concept of gauge coupling unification as an underlying organizing principle for the beyond-the-SM (BSM) model building and to systematically classify the BSM extensions that feature this property. In this way we intended to create a sort of road map that would guide the model building through the desert between the electroweak and the GUT scales. For starters, we looked for all possible extensions of the SM with two types of vector-like (VL) fermions in different representations of the SM gauge group. We found that the precise gauge coupling unification can be achieved for 13 distinct models only, as long as the masses of VL fermions are limited to 10 TeV. We also discussed complementarity of various experimental search strategies that can be employed to probe the BSM spectra favored by the unification requirement. We showed that all but one VL model can be excluded or tested by the current and future experiments.
The present study is a direct follow up of [30]. We extend the previous analysis by allowing a BSM theory to contain, beside VL fermions, also extra scalar fields and consider two broad categories of the SM extensions: with one type of fermions plus one type of scalars, and with two types of scalars in different representations. As we will see later, that increases significantly the number of possible unification scenarios as the impact of scalars on the renormalization group equations (RGEs) of the gauge couplings is milder than the one of fermions.
We also perform a detailed phenomenological analysis of all successful unification scenarios. We take into account a wide set of experimental collider searches, which include determination of the strong gauge coupling at different energy scales, measurements of the electroweak precision observables (EWPO), searches for heavy stable charged particles (HSCP), as well as the measurement of the proton lifetime at Super-and Hyper-Kamiokande. We update the analysis of [30] by incorporating the most recent experimental data and by adding experimental information that were not previously examined. The latter include projections for determination of the running strong and electroweak gauge couplings at a hypothetical future 100 TeV proton collider, which proves very efficient in testing the parameter space consistent with the unification condition.
The paper is structured as follows. In Sec. 2 we briefly revise our computational strategy of identifying BSM scenarios that allow for a precise unification of the SM gauge coupling constants. We also present a comprehensive list of all possible models that satisfy this condition. Sec. 3 is devoted to a detailed discussion of experimental searches that could test, and in many cases exclude, the parameter space of the models. We summarize our results and conclude in Sec. 4. Technical details of the analysis are collected in four appendices.

BSM models with gauge coupling unification
In constructing a set of generic anomaly-free extensions of the SM that allow for unification of the three gauge couplings of the SM, we strictly follow the approach introduced in [30]. First, we define the fundamental building blocks of our models, which are: • VL fermions F , which transform under the SU (3) C × SU (2) L × U (1) Y gauge group as

1)
• complex scalars S transforming as R S = (R 3S , R 2S , Y S ). (2.2) We denote by N F and N S number of copies of VL fermions and scalars that have the same quantum numbers. The left-and right-handed components of a VL pair are counted separately, so N F can only assume even values. The only exception is a fermion transforming in an adjoint representation of the non-abelian gauge symmetry group, in which case N F can be odd.
In [30] we discussed in details the SM extensions with two different VL fermion representations. In the present work we extend the previous analysis by incorporating scalar fields. We thus consider fermion plus scalar extensions of the SM, dubbed as 'FS' hereafter, which can be entirely parametrized by a set of properties and scalar plus scalar extensions, dubbed as 'SS', described by At this initial stage both N F and N S i are unconstrained and it is one of the goals of our analysis to determine their allowed values. The numerical procedure employed to establish which of the SM extensions with VL fermions and scalars could potentially lead to the precise unification of three SM gauge couplings at high energies was described in details in [30] and we briefly summarize it here. For a given BSM model, we initiate the renormalization group (RG) running of the SM gauge couplings at the scale associated with the top quark mass, M t , at which they assume the following values [32]: We then use the 2-loop SM RGEs given in [33] to run the system up from M t to the scale M 1 , at which the lightest of the BSM particles show up in the spectrum. To simplify the analysis we assume that all VL fermions (complex scalars) in the same representation R F (S) have a common mass. Above the scale M 1 we replace the SM RGEs with those corresponding to a generic BSM scenario [33], whose explicit forms adapted to the models considered in this study are given in Eqs. A.6-A.8 of Appendix A. Similarly, at the scale M 2 the impact of heavier BSM states is incorporated and the full system runs up. 1 The unification scale, M GUT , is then defined as the scale at which all three gauge couplings of the SM reach the same value, g GUT : Note that we use the SU (5) normalization, To quantify the precision of the gauge coupling unification we introduce a set of three mismatch parameters, 1,2,3 , defined as a deviation of one of the couplings, g k , from the unified value of the remaining two,   The putative unification scale is then calculated by imposing the requirement The last step is to set the upper bound on GUT , below which the gauge coupling unification will be defined as precise. Of course, any such choice is to certain extent arbitrary. For example, in the SM SM GUT = 7.3%, confirming the well known fact that the three gauge couplings do not really unify in this framework. In the minimal SUSY version of the SM, however, GUT can be much lower, reaching MSSM GUT = 1.1% if all sparticle masses are set at 1 TeV. We thus believe it is reasonable to define the precise gauge unification (PGU) by a condition GUT 1%. (2.9) Besides the PGU condition, that must be satisfied by any BSM model featuring the gauge coupling unification, there are also other requirements that need to be fulfilled by a model to make it a valid unification scenario. These extra requirements stem both from theoretical consistency of the theory and from the experimental data. We will now discuss them one by one and show that their inclusion significantly reduces the number of viable PGU models.
Perturbativity We will require that the running gauge couplings remain perturbative (g i 4π) up to the energy scale of 10 15 GeV, which is the proxy for the proton stability bound to be discussed in more details later on. This condition provides an upper bound on the possible dimensions of representations under which new fermions and scalar transform, and on the numbers of their copies N F and N S . In order to save computation time in the later stages of the analysis, it is convenient to discard all not perturbative models at the very beginning.
In the case of VL fermions, the corresponding perturbativity bounds were derived in [30] and we recall them here for completeness. In Table 1 we show the maximal allowed number of VL fermions, N Fmax , for which the SM gauge couplings g 2 and g 3 remain perturbative up to the unification scale. We set Y F = 0 (non-zero hypercharge could only reduce N Fmax ) and perform the analysis for different combinations of representations R 3F and R 2F . To be conservative and to avoid excluding viable scenarios, we fix the VL mass  at 10 TeV, which is the largest allowed mass scale we consider in this study. In Table 2 we report instead the maximal value of the hypercharge allowed by perturbativity of g 1 for a single pair of VL fermions (N F = 2). The corresponding results for complex scalars are presented in Table 3 and Table 4. Larger representations are allowed in this case due to the fact that scalar contributions to the gauge beta functions are by a factor of two smaller than the fermion ones, c.f. Eqs. A.9-A.11. We will see later that this property results in a larger number of the allowed PGU scenarios with at least one BSM scalar in comparison to the two-fermion case.
SU (5) unification Following [30], we assume that all BSM fields are embedded at the unification scale into the SU (5) multiplets. The decomposition of the irreducible SU (5) representations into irreducible representations of the SM gauge group is presented in Appendix B. Note that the representations of the dimension larger than 75 can be decomposed either to representations that have already appeared in the decomposition of lower dimensional representations, or to representations whose dimensions exceed the limits presented in Tables 1 and 3. Taking into account hypercharge perturbativity bounds from Tables 2  and 4, we are left with 24 fermion and 34 scalar where the superscript S indicates that a given representation is allowed for the scalars only.
In total there are 816 distinct perturbative models of type FS and 561 distinct perturbative models of type SS for which the possibility of precise gauge coupling unification needs to be investigated.
Proton decay Proton decay is a generic prediction of the GUT theories. In the framework of SU (5) the SM quarks and leptons belong to the same SU (5) multiplets. As a result, new heavy gauge bosons can mediate the interactions that violate both the baryon and lepton number conservation. The dominant contribution to the proton decay width stems from the dimension-6 operators of the type QQQL. While the exact form of these countributions is highly model-dependent, the approximate proton lifetime can be extracted from the unification scale M GUT and the value of the unified gauge coupling g GUT as [34] (2.10) Proton decay has been experimentally tested since the early 1990s by the underground water Cherenkov detector Super-Kamiokande (SK). The strongest lower bound on the proton lifetime comes from the decay channel p → e + π 0 and is given by [35] τ p > 1.6 × 10 34 years. (2.11) It is expected that in the not-so-far future the limit will be extended by at least one order of magnitude, up to ∼ 2 × 10 35 years [36], by a new generation detector Hyper-Kamiokande (HK). We apply the SK bound of Eq. 2.11 as an extra viability condition and discard all the models that do not satisfy it.

List of the allowed models
For each of 816 distinct models of type FS and 561 distinct models of type SS we scan over the number of BSM fermions and/or scalars, N F , N S , and their common masses, M 1 , M 2 . The upper bounds on the number of BSM particles are given in Tables 1 and 3, while the parameters M 1 and M 2 range from 0.25 TeV to 10 TeV. For each point in the 4-dimensional parameter space we determine GUT and M GUT from Eq. 2.8 and verify whether the PGU condition in Eq. 2.9 is satisfied. Finally, we check if the experimental bound on the proton lifetime, Eq. 2.11, is not exceeded.
The results of our numerical analysis are summarized in Table 5. In total, we found 22 different FS scenarios that allow for the PGU at the scale 10 15 − 10 18 GeV. The corresponding quantum numbers of fermions and scalar w.r.t. the SU (3) C × SU (2) L × U (1) Y gauge group are shown in columns 2 and 3. Many of the successful scenarios allow for various numbers of fermion and scalar copies, which are indicated in the fourth column of Table 5. If these two extra degrees of freedom are taken into account, one can identify 67 distinct FS-type PGU models, which we will denote as FSI (N F ,N S ) , with I= 1, . . . , 22.    (1, 1)  Table 5: FS and SS extensions of the SM, which allow for the PGU ( GUT ≤ 1%) at the unification scale 10 15 − 10 18 GeV and are not excluded by the SK measurement of the proton lifetime [35]. The BSM masses lie in the range 0.25 − 10 TeV. In columns 2 and 3 quantum numbers of the new particles w.r.t. the SM gauge group are given. Column 4 displays possible numbers of identical copies of fermions and scalars in each representation.
In the case of SS-type models, 17 distinct PGU scenarios were found. Note that the gauge coupling unification is not possible with the representations of a dimension higher than 8, even though these are allowed by perturbativity condition. Moreover, PGU can not be achieved with scalars in the adjoint representations of either SU (3) C or SU (2) L , contrarily to what is observed for VL fermions. Taking into account different numbers of scalar copies corresponding to the same representation, we can define 78 distinct models denoted as SSJ (N S 1 ,N S 2 ) , with J= 1, . . . , 17. They are all listed in Table 5. 2 The unification requirement does not only select specific representations of the SM gauge symmetry groups under which the BSM particles should transform, but it also determines the spectrum, i.e. a particular hierarchy between the mass parameters M 1 and M 2 . In this regard, the successful PGU scenarios can be divided in four categories, which will be referred to with the following labels hereafter: An illustration of different types of the PGU spectra is presented in Fig. 1. Each of the four panels shows a distribution of the mismatch parameter GUT as a function of the parameters M 1 and M 2 for a specific SSJ (N S 1 ,N S 2 ) scenario from Table 5. The PGU region, in which GUT < 1%, is marked in red. Different shades of yellow indicate departure from the precise unification condition as quantified by the increasing values of GUT , whose upper bound is given in the legend bar on the right hand side of the plot. Dashed black contours indicate the unification scale M GUT (in GeV). Types of spectra characteristic of a given PGU scenario are also indicated in the second column of Tables 6 and 7 in Appendix D.
In [30] we demonstrated that a particular hierarchy among the mass parameters required by the unification condition allows one to probe many of the PGU scenarios of the type FF up to at least 10 TeV thanks to complementarity of various experimental search strategies. In the present study we will assume the same approach to derive experimental bounds on the PGU models listed in Table 5.

Testing the PGU models at the colliders
In general, there are two broad categories of experimental searches that one may want to employ to probe the PGU models identified in Sec. 2. The first category encompasses those searches that provide the exclusion bounds which are model independent, in the sense that once the gauge properties of a BSM scenario are defined (in other words, once a particular model from Table 5 is chosen), the experimental limits would not depend on further  assumptions regarding the model properties, like the presence of other (non-gauge) interactions. Contrarily, the predictions based on the searches belonging to the second category strongly hinge on those additional assumptions. In particular, in Sec. 3.2 the experimental lower bounds on fermion and scalar masses will be derived under the assumption that the Yukawa interactions between the BSM and the SM particles are negligible. A similar analysis was performed in [30] for the FF-type PGU models. In this work we extend the previous study to the scenarios of FS and SS type. Additionally, we update some of the experimental searches used in [30] to incorporate the more recent results. Finally, we discuss the projections derived from putative measurements of the running weak and strong gauge coupling constants at the future 100 TeV collider.

Model-independent tests
Running of the strong gauge coupling The deviation of the renormalized strong gauge coupling constant from its SM running is determined by the representation and the number of BSM particles that transform non-trivially under the SU (3) C gauge group. At one-loop the effect is controlled by the size of corrections to the coefficient B 3 of the corresponding beta function, Eq. A.6, which reads Last year the ATLAS collaboration released a new determination of the running g 3 (Q), based on the measurements of the transverse energy-energy correlations (TEEC) and their associated azimuthal asymmetric (ATEEC) functions, using a data sample containing 139 fb −1 of proton-proton collisions at the center-of-mass energy of √ s = 13 TeV [37]. The new determination of g 3 (Q) shows a very good agreement with the SM predictions up to the energy scale of about 4 TeV. These results updated the previous ATLAS analysis of the same type, based on the √ s = 8 TeV data set [38], in which the running of the strong coupling constant was tested up to around 1 TeV; and the strong coupling determination based on the measurement of the double-differential inclusive jet cross section at √ s = 8 TeV by CMS [39], in which the consistency of the g 3 (Q) with the SM prediction was confirmed up to 1.5 TeV.
The experimental searches measuring the energy dependence of g 3 use as the main kinematical variable the leading jet transverse momentum p T . As a result, the maximal energy scale at which the running coupling can be evaluated is limited by the value of p T . According to the report [40], which summarizes most important properties of SM processes at the future 100 TeV proton collider, the reach in p T at such a machine could extend well above 20 TeV, assuming the total integrated luminosity of 1 ab −1 . To derive a (very approximate) projection for the future impact of the running strong coupling measurement on our models, we will assume that g 3 does not deviate from the SM prediction up to at least 20 TeV and will use this number as a rough approximation for the maximal energy scale at which it can be probed.
The impact of the strong coupling determination on the PGU scenarios, derived using the combination of ATLAS measurements with √ s = 8 TeV and √ s = 13 TeV, is illustrated in Fig. 2 for two exemplary cases that are either (a) excluded by the measurement of g 3 (up to 10 TeV), or (b) can be tested by this measurement in the future. The corresponding exclusion bounds are indicated as green solid lines, while the 100 TeV projections as green dotted lines. Dark gray shaded area corresponds to the PGU region of the parameter space, GUT 1%. The meaning of the remaining lines will be explained later on. The lower bounds on the masses of color particles in all the PGU models listed in Table 5 are also collected in column 3 of Tables 6 and 7 in Appendix D. Note that in our previous study [30] the CMS search [39] was employed, and thus the present exclusion bounds are significantly stronger. The corresponding lower bounds based on the 100 TeV projection are shown in column 9 of Tables 6 and 7. As it was noted in [30], the running of g 3 is particularly efficient in testing the hierarchical PGU spectra in which the color    (40,6) and (b) FS6 (10,2) . In dark gray the PGU region is indicated. The area below and left to the solid green line is excluded by the measurement of the running strong coupling constant by ATLAS [37]. The limits from the 13 TeV ATLAS R-hadrons search [41] are indicated as a dashed green line. The corresponding lepton-like HSCP search excludes the area left and below to a dashed red line. 100 TeV collider projections for the running g3 [40], running g2 [42], and for the EWP tests [43] are depicted as green, brown and red dotted lines, respectively. particle in much lighter than the non-color one, allowing one to exclude, or to probe in the future, many of the models characterized by the spectrum type H1 or H2.
Running of the electroweak gauge couplings A similar strategy could in principle be applied to the electroweak couplings g 1 and g 2 . While the running of neither of them has been experimentally measured yet, in [42] it was proposed to look at the high invariant mass spectrum of the Drell-Yan (DY) production p p → γ/Z * → l + l − measured by the LHC to derive the limits on the energy dependence of g 1 and g 2 . Since the available experimental data from LEP and the 7 TeV LHC run did not produce strong constraints, [42] provided predictions for 8, 14 and 100 TeV. By now, the actual neutral current DY cross section was measured by CMS [44] and ATLAS [45] at 8 TeV, and by CMS at 13 TeV [46] with 2.3 fb −1 of data, finding good agreement of the measured differential cross sections with the theoretical calculations. Therefore, we will treat the 8 TeV results of [42] as the actual limits.
The DY constraints on g 1 and g 2 at a given energy scale can be translated into the limits on one-loop BSM contributions to the corresponding beta functions, as a function of the BSM mass. Since the process p p → γ/Z * → l + l − is sensitive to both g 1 and g 2 , the limits in [42] are provided separately on ∆B 2 under the assumption that (a) Testable by g2 -model FS3 (3,27) (b) Testable by g2 -model FS16 (4,41) Figure 3: Summary of experimental bounds on masses M1 and M2 for scenarios (a) FS3 (3,27) and (b) FS16 (4,41) . The color code is the same as in Fig. 2.
there is no BSM contribution to the running of the hypercharge coupling (∆B 1 = 0), and on ∆B 1 assuming that ∆B 2 = 0. We applied those limits to all the PGU models listed in Table 5 and found that none of them is presently constrained by the DY production. We can use, however, the results of [42] to quantify the future reach of the g 2 measurement at the future collider with the centre-of-mass energy of 100 TeV and 3 ab −1 of data. The strongest exclusion bounds can be derived in this case from the charged current DY process p p → W * → lν, which also makes it possible to constrain ∆B 2 individually, with no further assumptions regarding the U (1) Y charges. The corresponding projected lower bounds on the BSM fermion and scalar masses in our PGU models are reported in column 8 of Tables 6 and 7.
In Fig. 3 we present two PGU scenarios whose parameter space could be entirely tested by the future measurement of g 2 in the charged current DY process. The corresponding lower bounds on fermion and scalar masses are indicated as dotted brown lines. One observes that the g 2 limits are complementary to the limits form the strong gauge coupling measurement and could provide a constraint on those models in which non-color particles are lighter than the color ones. The bounds, however, are significantly weaker due to the lower production cross section of the EW gauge bosons at hadron colliders.
Electroweak precision tests Bounds on the BSM sector can also be derived by looking at the DY processes below the BSM mass threshold and by testing deviations from the SM predictions in the EWPO [43]. The oblique parameters sensitive to the presence of states charged under the EW gauge symmetry are W and Y [47,48], related to the corresponding beta functions as [42] where ∆B 2,1 are given in Eq. 3.2.
In [43] the limits on W and Y were derived based on the data from LEP [49], and from the LHC 8 TeV measurements by ATLAS [45] and CMS [44]. We checked that they do not provide any bounds on the parameter space of the PGU scenarios from Table 5. However, as it was the case for the running coupling g 2 , the present experimental bounds can be improved at the projected 100 TeV machine with the luminosity 3 ab −1 by roughly two orders of magnitude [43]. The corresponding projections are summarized in column 7 of Tables 6 and 7. They are also indicated as dotted red lines in Figs. 2 and 3.
Both the EWPO and running g 2 bounds are particularly powerful in constraining the BSM scalars and fermions that transform non-trivially under SU (2) L , and thus they can probe the PGU spectra characterized by the non-color particle lighter than the color one. It is worth to emphasize that all three model-independent searches discussed in this subsection provide a set of complementary constraints that proves to be particularly efficient in testing the hierarchical spectra of the H1 and H2 type. As a matter of fact, a close inspection of Tables 6 and 7 reveals that these are exactly the PGU models that are excluded or testable by the model-independent collider searches.
Proton lifetime Although technically not being a collider search, the proton decay measurement can put strong constraints on the models with gauge coupling unification. As already mentioned in Sec. 2, it is expected that in the not-so-far future the present SK limit will be extended by at least one order of magnitude, by a new generation detector Hyper-Kamiokande [36]. In column 6 of Tables 6 and 7 we report the maximal predicted proton lifetime in years, τ max p , in the PGU region of the corresponding model. If this number is lower than the bound of Eq. 3.4, the model will be entirely tested by the HK measurement. This is exactly the case of scenarios FS1 (18,6) and SS7 (5,6) whose parameter space is presented in Fig. 4. The projected exclusion by HK is indicated as blue dashed lines, while a blue solid line denotes the present day SK bound.
The PGU scenarios in Fig. 4 are characterized by the spectra of type H3 and H0, correspondingly. Here both BSM particles can be very heavy and for this reason the PGU parameter space is difficult to probe with the collider searches. The proton lifetime measurement, however, can constrain scenarios with a relatively low unification scale independently on the BSM fermion and scalar masses, and thus can provide yet another complementary way of testing the parameter space consistent with the precise unification condition.

LHC searches for long-lived particles
In this study we make an assumption that the Yukawa interactions of the BSM fermions and scalars with the SM quarks and leptons are negligibly small, so that they can be treated as    (18,6) . The color code is the same as in Fig. 2. Blue solid line marks the exclusion by the proton decay measurement at Super-Kamiokande [35]. Blue dashed lines indicate the projected reach of Hyper-Kamiokande [36].
very long-lived, or quasi-stable particles which do not decay inside the detector. 3 Dedicated analyses looking for such heavy stable charged particles (HSCP) were performed both by the ATLAS and CMS collaborations, based on the observables related to ionization energy loss and time of flight, characteristic of heavy particles whose velocities are significantly smaller than the speed of light.
Color HSCPs It is usually assumed that if the lifetime of a color HSCP is longer than the typical hadronization time scale, it would form a colorless QCD bound state with the SM quarks and gluons, the so-called R-hadron. The most recent 95% confidence level (C.L.) upper bound on the color HSCP production cross section, which can be translated into a lower bound on the HSCP mass, were derived by ATLAS using a data sample corresponding to 36 fb −1 of proton-proton collisions at √ s = 13 TeV [41], and by CMS using 2.5 fb −1 of data at the same energy [50]. While in both analyses the HSCP under study is gluino (a supersymmetric partner of the gluon), the cross-section bounds can be easily applied to any color HSCP. The reason is that hadronic interactions of R-hadrons with detector material proceed through its lightquark constituents, and thus are to a good approximation independent on the original HSCP quantum numbers.
In hadron colliders a color particle can be pair produced at the leading order (LO) through gluon fusion and quark-antiquark annihilation. A detailed calculation of the corresponding p p → XX cross section for a BSM particle X in a generic representation R 3 of SU (3) C is presented in Appendix C. Comparing those results with the observed exclusion limits on the cross section reported by ATLAS [41], we can derive lower bounds on the masses of colored particles in our PGU models. Those bounds are shown in column 4 of Tables 6 and 7.
In Figs. 2-4 the bounds from the R-hadron searches are shown as green dashed lines, to emphasize their dependence on the negligible Yukawa coupling assumption. In the presence of non-zero Yukawa interactions, color BSM fermions and scalars would decay before a bound state is formed and the limits would not apply anymore. This is in striking contrast with the bounds from the running strong coupling constant, which are model-independent once the SU (3) C charges are fixed.
Note that at the moment the R-hadron limits are stronger (often significantly) than those from the running g 3 , and so it is due to them that many of the PGU scenarios can be excluded. We indicate this fact with an asterisk in column 4 of Tables 6 and 7.
EW-charged HSCPs If a HSCP is a SU (3) C singlet but it carries non-zero EW charges, it will be produced through Drell-Yan (DY) processes and will predominantly lose energy via ionization inside the detector. The corresponding LO cross section for a BSM particle with arbitrary SU (2) L and U (1) Y quantum numbers is calculated in Appendix C. The lower bounds on the BSM particle masses can be derived using the chargino-and tau-dedicated HSCP searches by ATLAS with 36 fb −1 of data [41].
The corresponding exclusion bounds are summarized in column 5 of Tables 6 and 7 and depicted in Figs. 2-4 as dashed red lines. At the moment they provide the only experimental way of testing at the colliders the PGU scenarios in which the non-color particle is lighter than the color one. In the future competing bounds, independent on values of the BSM Yukawa couplings, should be provided by the EW searches discussed in Sec. 3.1.

Synopsis of the results
In this section we performed a detailed analysis of various experimental signatures characteristic of the PGU scenarios identified in Sec. 2.
We found that out of 67(78) models of type FS(SS) that allow for the precise gauge coupling unification, 23 (35) are already excluded (up to 10 TeV) by the combination of the running strong coupling measurement and the HSCP searches. Those models are dubbed as "excluded" in column 10 of Tables 6 and 7.
We have also shown that many of the remaining scenarios, 23 in the FS case and 25 in the SS case, could be further probed at the future 100 TeV proton collider, or by the proton lifetime measurement at HK. Those models are dubbed as "testable" in column 10 of Tables 6 and 7. Note that they almost exclusively belong to the categories H1 and H2, with a few H0 and H3 exceptions that can be either tested by HK, or are border-line cases between different types of spectra.
Finally, there are 21(18) models of type FS(SS) that remain very challenging to probe experimentally. Those are predominantly of the H3 and H0 spectrum type, i.e. featuring the parameter space where both BSM particles can be heavy. Moreover, in many of these models the unification scale is very high, M GUT > 10 17 GeV, so there is not much hope to probe them in HK-like experiments either.
Nevertheless, it feels encouraging that the great majority of the PGU models identified in this study will be experimentally probe in the next decades.

Conclusions
In this study we performed a systematic analysis of the SM extensions with BSM scalars and fermions in different representations of the SM gauge group in the context of the gauge coupling unification. We analyzed all possible combinations of one representation of scalars and one representation of fermions, as well as of two different representations of scalars, where the number of fermions and scalars in each representation was limited by demanding perturbativity of the gauge couplings at the unification scale.
We found 145 different two-representation BSM extensions that allow for precise gauge coupling unification at the energies higher than 10 15 GeV, under the assumption that all BSM masses are in the range 0.25−10 TeV. The successful models are listed in Tables 6 and  7 of Appendix D, together with the information about the type of spectrum favored by the unification condition. The latter is particularly important in the context of experimental testability of the unification models. In particular, those PGU scenarios which feature hierarchical spectra, with one of the BSM particles significantly lighter than the other, are very susceptible to collider searches based on heavy stable charged particle signatures, as well as to direct measurements of the SM gauge coupling scale dependence.
We also found that complementarity of various experimental strategies plays a crucial role in probing the parameter space of the unification scenarios. Out of 145 possible models, 40% can already be excluded based on the present data, while 33% more may be probed in the future at the 100 TeV collider. This observation is in line with the findings of [30], where VL fermion extensions of the SM were analyzed.
We hope that the results of our study would add yet another small incentive in favor of building a more powerful proton collider.

A Group invariants and loop coefficients
Let us consider a symmetry group G. The quadratic Casimir operator for the representation R is defined as where t A are the generators of G in the representation R. The Dynkin index of a representation R is instead given by The two are related through the dimension of the adjoint representation d(Adj), For convenience, one can parametrize the group-theoretical factors through the weights (p, q) for an irreducible SU (3) representation R 3 , and through the highest weight for an SU (2) representation R 2 , General two-loop beta functions for a system of gauge couplings g i of a direct-product symmetry group G i × · · · read [33] where G i , R F and R S indicate contributions from gauge bosons, Weyl fermions, and complex scalars respectively. The sum is meant in both S(R F ) and S(R S ) over all fermion and scalar representations transforming nontrivially under G i . The two-loop beta functions for SM augmented with N F Weyl fermions in the representation (R F 3 , R F 2 , Y F ) and N S complex scalars in the representation (R S3 , R S2 , Y S ) can be straightforwardly derived from Eq. A.5, and read with the one-loop coefficients defined as The two-loop coefficients are instead given by

C Production cross section of the BSM fermions and scalars at the LHC
In this section we derive analytical formulae for the pair-production cross section of color and non-color fermions and scalars at the LHC. We further confront our results with the output of numerical calculations performed with MadGraph5 aMC@NLO [52].

C.1 Color particles
Pair production of heavy color particles at hadron colliders proceeds dominantly via strong interactions. At the LO, two classes of processes contribute to the parton-level production cross section: quark-antiquark annihilation and gluon fusion. The numerical value of the cross section is then entirely determined by three properties of the produced particle: its spin, its mass, and its transformation properties under the SU (3) C symmetry group, encoded in the representation R 3 .
Scalars The only process contributing to the quark-antiquark production of a color complex scalar S, qq → S S * , is depicted in Fig. 5, and its partonic cross section readŝ where β S = 1 − 4m 2 S /ŝ, α s = g 2 3 /(4π), m S is the mass of the scalar S, the Mandelstam variable readsŝ = (p 1 + p 2 ) 2 = (q 1 + q 2 ) 2 , and p i , q i denote the momenta of the initial and final state particles, respectively. "Fun" and "Adj" stand for the fundamental and the adjoint representation of SU (3) C , correspondingly.
On the other hand, a contribution to the scalar pair production stemming from the gluon fusion, g g → S S * , is less straightforward to calculate. The LO Feynman diagrams are depicted in Fig. 6 and the corresponding amplitudes read, from left to right, In the above T a are the generators of SU (3), f abc are the corresponding structure constants, and t = (p 1 −q 1 ) 2 and u = (p 1 −q 2 ) 2 are the Mandelstam variables. The objects M i capture the kinematics of the process and are given by where is a polarization vector of the gluon and the contraction of the Lorentz indices is understood. The calculation of the total g g → S S * amplitude can be simplified by employing the factorization property of the non-Abelian gauge theories, according to which the amplitude of the process can be expressed as a product of two terms, M = α=s,t,u,tu where G encodes the properties of the produced particle with respect to the non-Abelian gauge group, and M A involves the dynamical part corresponding to the pure Abelian process. To derive the explicit forms of both factors, let us first note that the total g g → S S * amplitude can be written as [53,54] M = α=s,t,u It is now easy to verify that the factors appearing in Eq. C.11 satisfy the following relations, From Eq. C.11 and Eq. C.15 it results that the group theoretical factor G and the Abelian amplitude M A in Eq. C.10 read Summing over the final and averaging over the initial degrees of freedom (colors and spins), one obtains the averaged squared amplitude and Finally, the partonic cross section for the process g g → S S * can be written aŝ The above formula reproduces the results of [54,55] for color triples and sextets, correspondingly.
Fermions The derivation of the total partonic cross sections for the processes qq → FF and g g → FF proceeds analogously to the scalar case. The corresponding Feynman diagrams can be straightforwardly obtained from those depicted in Figs. 5 and 6 by replacing the scalar external lines with the fermion ones (the only exception is the rightmost diagram in Fig. 6, which does not exists for the fermions). The factorization conditions Eq. C. 10 and Eq. C.15 are satisfied for the fermion case as well, although the kinematical factors T α assume a slightly different form. In the end, the total partonic cross sections read  21 reproduces the results of Ref. [56] for color triplets.
Hadronic cross section The hadronic cross sections σ(p p → S S * ) and σ(p p → FF ) are obtained by convoluting the corresponding partonic cross sections with the parton luminosties L AB ab as where τ =ŝ/s, τ 0 = 4m 2 S /s, s is the hadronic center of mass energy squared and f A a (x, µ) indicates a parton distribution function for the parton a in the hadron A with the factorization scale µ.
In Fig. 7 we show the hadronic cross section at the LHC with √ s = 13 TeV for the pair production of fermions (left panel) and scalars (right panel) as a function of their mass. We assume in both cases that the produced particle transforms as a SU (3) C triplet. Red dotted curve corresponds to the production via quark-antiquark annihilation, while the blue dashed one indicates the gluon fusion. One observes that for fermions lighter than ∼ 600 GeV the gluon fusion is a dominant contribution to the total cross section, while for the heavier ones the quark-antiquark annihilation becomes larger and the cross section scales as S(R 3 ). A similar pattern is observed in the scalar case, although the corresponding cross-over between the dominant contributions takes place for the scalar mass of around ∼ 2.5 TeV.
In Fig. 8 we present a comparison of the total hadronic cross sections for the pair production of fermions (left panel) and scalars (right panel) in different SU (3) C representations: triplet (dotted red), sextet (dashed blue), octet (solid green) and tenplet (longdashed purple). The most important observation, which agrees with the findings of [55], is that the cross sections for sextets and octets are almost the same, as a direct consequence of similar numerical values of the group theoretical factors C(R 3 ) and S(R 3 ).
Finally, in Fig. 9 we present a comparison of our analytical derivation of the tree-level hadronic sections based on Eq. C.1 and Eqs. C.20-C.23, corrected with the k-factor of 1.5 [57] to account for higher-order QCD corrections, with the results of the numerical simulation performed with MadGraph5 aMC@NLO. Once more, the left panel illustrates the case of the fermions and the right one of the scalars. SU (3) C triplets are indicated as solid red lines, while the SU (3) C octets as solid blue lines. Dotted and dashed grey lines denote the MadGraph5 aMC@NLO results, correspondingly. One can see that our results reproduce the numerical simulation quite precisely, in particular at the heavier end of the spectrum.

C.2 EW-charged particles
At the tree -level, pair production of heavy non-color particles at hadron colliders proceeds through the Drell-Yan processes depicted in Fig. 10. Averaging over the initial color and spin states, and summing over the final spin states, one finds In the following we present individual averaged squared contributions to Eq. C.26 for the scalars and fermions separately.
Fermions After the EWSB, the neutral current interactions of a generic fermion F with the photon A and the gauge boson Z are given by where couplings of the fermion with Z are entirely determined by its SU (2) L and U (1) Y quantum numbers as Here T 3F denotes the third component of the weak isospin, θ W is the Weinberg angle, and the electric charge Figure 10: Feynman diagram for the LO fermion and scalar pair production at hadron colliders via the Drell-Yan process.
The first term in Eq. C.26 resembles the matrix element for the process qq → FF with the gluon exchange in the s-channel. It reads where α = e 2 /(4π) and Q q indicates the charge of the SM quark. The s-channel exchange of the boson Z leads to where M Z denotes the mass of the Z boson, and Γ Z its total decay width. Finally, the interference term is given by Consequently, the partonic cross section for the pair production of colorless fermions readŝ The above formula reproduces the results of [58] for the SU (2) L singlets with an electric charge 1.
In the left panel of Fig. 11 we present the hadronic cross section at the LHC with √ s = 13 TeV for the DY pair production of EW-charged colorless fermions as a function of their mass. Solid green line indicates an SU (2) L singlet with the hypercharge Y F = 1, dashed blue line an SU (2) L doublet with Y F = 1 2 , and red dotted line an SU (2) L triplet with Y F = 1.
Scalars In the case of spin-zero particles, the relevant Lagrangian for the EW interactions reads where Q S denotes the scalar electric charge, Q S = T 3S + Y S . Analogously to the fermion case, the coupling between an EW-charged scalar and the Z gauge boson is defined as Likewise, three contributions to the averaged squared amplitude of the process qq → S S * are given by (C.37) Finally, the total partonic cross section for the DY production of the EW-charged scalar particles readŝ A comparison of the total hadronic cross sections for the EW pair production qq → S S * of colorless scalars in different representations of SU (2) L is presented in the right panel of Fig. 11. D List of the PGU models and experimental constraints Table 6: Scenarios with one VL fermion representation and one scalar representation that allow for the PGU ( GUT ≤ 1%) and the associated unification scale lies in the range 10 15 − 10 18 GeV. The VL masses vary between 0.25 TeV and 10 TeV. In column 2 type of the model spectrum is given, as defined in Eq. 2.12. In columns 3, 4 and 5 current lower bounds on the BSM masses (in TeV) are presented, as provided by the measurement of the running coupling g3 [37], R-hadrons [41] and charged HCSP [41] searches, respectively. In column 6 the maximal proton lifetime (in years) is quoted in the PGU region of a given model. Projected sensitivities to the EWPO [43] tests, the measurement of the running coupling g2 [42] and determination of g3 [40] at the future 100 TeV collider are collected in columns 7, 8 and 9. Finally, in column 10 we indicate whether a given model is excluded by the current experimental data, or if it will be entirely tested in the future. The experimental bound responsible for the exclusion (testability) of a model is indicated with an asterisk.