Radiative Effects in the Scalar Sector of Vector Leptoquark Models

Gauge models with massive vector leptoquarks at the TeV scale provide a successful framework for addressing the B-physics anomalies. Among them, the 4321 model has been considered as the low-energy limit of some complete theories of flavor. In this work, we study the renormalization group evolution of this model, laying particular emphasis on the scalar sector. We find that, despite the asymptotic freedom of the gauge couplings, Landau poles can arise at relatively low scale due to the fast running of quartic couplings. Moreover, we discuss the possibility of radiative electroweak symmetry breaking, and characterize the fine-tuning associated with the hierarchy between the electroweak scale and the additional TeV-scale scalars. Finally, the idea of scalar fields unification is explored, motivated by ultraviolet embeddings of the 4321 model.


Introduction
Recently, there has been a resurgence of interest in the study of leptoquark fields in light of discrepancies in semi-leptonic B-decays, collectively known as B-physics anomalies [1][2][3][4][5][6][7][8][9][10], which challenge the assumed lepton flavor universal nature of fundamental interactions. In particular, the vector leptoquark U 1 with Standard Model (SM) quantum numbers (3, 1) 2/3 , originally proposed in Ref. [11][12][13], has been identified as the only single-mediator solution [14,15]. The search for an ultraviolet (UV) completion of the U 1 simplified model naturally leads to variations on the Pati-Salam (PS) gauge group, G PS ≡ SU (4) × SU (2) L × SU (2) R [16], as it features quark-lepton unification. The U 1 arises as a massive gauge boson after spontaneous symmetry breaking of the PS group to the SM. However, the original PS model predicts a leptoquark with flavor-universal couplings which has to be very heavy in order to avoid the strict bounds from processes involving the light generations.
The minimal variation containing the U 1 leptoquark at TeV scale while being phenomenological viable is G 4321 ≡ SU (4) × SU (3) × SU (2) L × U (1) X [17,18], where SU (2) L is identified with SM weak isospin, and both color and hypercharge are diagonal subgroups (for examples of other variations, see Refs. [19][20][21][22]). Models based on G 4321 are called 4321 models and generally fall into two categories based on the charge assignments of the SM matter fields: i) flavor universal 4321 [18,23], where all would-be SM fermions are singlets under SU (4) and have the usual SM charges and ii) flavor non-universal 4321 [24][25][26][27][28][29][30], where the would-be third generation quarks and leptons are unified into SU (4) multiplets. This non-universality in charge assignments translates into non-universality of the couplings of the new heavy vector boson, among them the U 1 , to the different generations. The flavor non-universal models can also be viewed as the low-energy limit of more fundamental theories that envision larger symmetry groups broken at various scales and aspiring to dynamically explain the SM Yukawa hierarchies (see e.g. the three-site PS model [31], PS variations in five dimensions [32,33] or the twin PS model [34]). In the following, we will refer to UV models reducing to the 4321 model at low-energy as UV embeddings of 4321.
For both types of 4321 model, the spontaneous symmetry breaking to the SM requires the presence of new scalars at the TeV scale. The purpose of this work is to investigate the Renormalization Group (RG) evolution of the parameters in the extended scalar sector of the 4321 model. The system of RG equations together with boundary conditions is in general under-determined, since the number of parameters is larger than the number of formal and phenomenological bounds, even in the most simplified version of the model. To study this system, we employ different benchmark scenarios for the boundary conditions. We consider two types of scenarios: i) agnostic towards further assumptions in the UV, and ii) inspired by UV embeddings. The former naturally allows for more freedom in the Z ∼ (1, 1) 0 , with masses given by where ω 1 and ω 3 are the Vacuum Expectation Values (VEVs) of the scalars responsible for the spontaneous symmetry breaking (see Section 2.2).
ii) ω 3 > ω 1 as in Ref. [23] which is a preferred phenomenological limit to avoid significant contributions to ∆F = 2 observables. In this case and in the limit g 1,3 g 4 , we obtain While the limit (i) is particularly motivated in composite scalar sector (see for example Ref. [27]), the limit (ii) offers better compatibility between the case of fundamental scalars and phenomenology. For completeness the couplings of the SM gauge bosons, as well as of the heavy bosons U 1 , G and Z are given by [18] g s = g 4 g 3 g 2 4 + g 2 3 , g Y = g 1 g 4 g 2 4 + 2 3 g 2 1 , g U = g 4 , g G = g 2 4 − g 2 s , g Z = 1 2 (2.2) Fermion content. In the non-universal 4321 model, only the 3 rd generation of quarks and leptons are charged under SU (4). The light generation fields, q i L , u i R , d i R , i L , e i R with i = 1, 2, are singlet under SU (4) and have the same quantum numbers as in the SM, with hypercharge as their U (1) X charge and SU (3) c representation as their SU (3) ones. The SU (3) of 4321 can therefore be seen as the color group for light generations. The model also contain one or more vector-like fermions χ j L,R , where j runs over the number of vector-like fermions. After being integrated out, they generate couplings between the 3 rd and the light generations of the SM fermions. The 3 rd generation fermions are unified into quadruplet with quarks (partners) and lepton (partners) as follows The gauged U (1)X is then a combination of U (1) and another U (1), see Ref. [27].
The vector-like fermions (Q j L,R , L j L,R ) and the 3 rd generation fields (q 3 L , 3 L ) in the broken phase of 4321 have the same SM charges as the light generations fields. In the following, we will consider only one vector-like field χ L,R (j = 1).

Scalar content and Lagrangian
The scalar content of the 4321 model used in the analysis is presented in Table 1. The scalar potential is Notice that the ρ 5 term in Eq. (2.5) breaks explicitly a U (1) global symmetry of the scalar potential, thus preventing the appearance of a massless Goldstone boson after spontaneous symmetry breaking.
The Higgs potential is the same as in the SM 6) and the mixing terms between the Higgs sector and the Ω sector are The roles of each scalar are as follows (the explicit Yukawa terms are written in Appendix B.1): • Ω 3 breaks 4321. After getting a vev, since this scalar is charged under both SU (3) light and SU (4) heavy, it will mix the light generation quarks with the vector-like quarks. Therefore, it is responsible for the mass mixing between light quarks and the 3 rd family quarks, and also for the small coupling between the U 1 leptoquark, the colorons and the Z with the light quarks.
• Ω 1 breaks 4321. After getting a vev, since this scalar is charged under SU (4) heavy and has the right U (1) X charge, it will mix the light generation leptons with the vector-like leptons. It has the same role as Ω 3 but for leptons instead of quarks (mass mixing between light and 3 rd family and small coupling of light leptons to heavy gauge boson from 4321 breaking).
• H breaks the electroweak symmetry. This scalar is the same as the SM Higgs. It gives mass to all quarks and charged leptons.
The two scalar fields Ω 1,3 break spontaneously G 4321 → G SM by acquiring the VEVs The β-functions for the couplings of the potential in Eq. (2.4) are given in Appendix A. The scalar content presented in Table 1 is the minimal content needed for the analysis. More scalars might be needed for the phenomenology of the 4321 model [26], however they are not expected to alter significantly the conclusions of this analysis (see Appendix B.3).

Radial and Goldstone modes
The Ω 1,3 scalars can be decomposed as where 1 3 is the 3 × 3 identity matrix and t a are the SU (3) generators. The complex scalars have the following representations under the SM gauge group: The scalar spectrum contains the following mass eigenstates: • two real octets: the massless O a GB = √ 2 ImÕ a 3 and its massive orthogonal combination O a R ; • two complex triplet: the massless T GB = s βT1 − c βT3 , with tan β = ω 1 /ω 3 , and its massive orthogonal combination T R ; • four real singlets, all combination ofS 1 andS 3 : the massless S GB and three massive radials S 0 , S 1 , S 2 .
The Goldstone boson modes O a GB , S GB , and T GB fields are eaten by the coloron G , the Z , and the leptoquark U 1 , respectively and the radial modes acquire the following masses: with u defined as (2.15)

Electroweak symmetry breaking
In this section, we identify the conditions for the spontaneous symmetry breaking of the electroweak sector. Writing the potential Eq. (2.4) in terms of the radial degrees of freedom aligned along the vacua, we obtain where s 1 = ω 1 √ 2 + Re(S 1 ) and s 3 = 3 2 ω 3 + Re(S 3 ).
A necessary requirement for our potential to have a minimum at h = v, s 1 = ω 1 √ 2 and s 3 = 3 2 ω 3 is that the first derivative of the potential evaluated at this point vanishes, i.e.
These three equations give the set of conditions: where for each the vanishing VEV solution has been ignored. From this point forth, since the smallness of ρ 5 is technically natural, we neglect its contribution to the effective mass and quartic coupling of the Higgs. In the limit where v ω 1 , ω 3 , the term proportional to v 2 can be omitted in the second and third equations in Eq. (2.18). The Hessian matrix The squared mass matrix has to be positive definite, which, for a symmetric matrix, is equivalent to all its leading principal minors being positive (2.20) The other minors are guaranteed to be positive from the inequalities above, but it is useful to write them explicitly as they put clear constraints on the mixing parameters η 1 and η 3 , namely In the limit where v ω 1 , ω 3 , the eigenvalue of the square mass matrix Eq. (2.19) corresponding to the physical Higgs mass is .

(2.22)
Given the positive definite conditions Eq. (2.20), it is clear from its definition that λ eff > 0. Moreover, we can define an effective Higgs mass parameter as The tree-level threshold corrections from integrating out s 1 and s 3 in the potential Eq. (2.16) yields an effective potential with the same expressions of λ eff and µ 2 eff as in Eq. (2.22) and Eq. (2.23), respectively. Notice the equivalence between the squared mass matrix diagonalization Eq. (2.19) with a VEV hierarchy approximation and the tree-level threshold corrections by integrating out the heavy states. 2 The EWSB happens when µ 2 eff flips sign, from positive to negative.

Bounded from below conditions
In order for the potential to be Bounded From Below (BFB) and thus for the vacuum to be stable, there must be restrictions on the parameters in the quartic part of the potential V 4 . For the potential V Ω , the BFB conditions have not previously been derived in the literature, and a full stability analysis of the potential is beyond the scope of this paper. A rigorous treatment of vacuum stability for similar models using the copositivity criteria can be found in Refs. [45,46]. Below we list some conditions, obtained by studying the behavior of V 4 along specific directions of the fields Ω 1 and Ω 3 and minimizing the potential against extra parameters such as θ and α as in Ref. [47]. In what follows, we consider strong stability requirements, i.e. V 4 > 0 for arbitrary large values of fields components. Stability in the marginal sense V 4 ≥ 0 (if V 4 = 0 then one must also require the quadratic part V 2 ≥ 0) is obtained by making the following inequalities non-strict. For the Ω potential, we obtained , with r → ∞. For the Higgs potential, the well-known condition is λ > 0. For the mixed part of the potential V ΩH , using the parametrization |H| 2 = r sin θ, we obtained These relations are necessary, but not sufficient, conditions to ensure the vacuum stability.

Analysis of the RG equations system
In this section, we study the RG equations system of the model. In what follows, we explain the procedure for solving this system and list the constraints at the IR boundary.
Methodology. We solve numerically the system of RG equations using the β-functions for the 4321 model presented in Appendix A (neglecting the terms related to the second Higgs doublet) from the scale Λ UV to the matching scale Λ 4 , and the SM β-functions from Λ 4 down to m Z . Starting with the full particle content, we first integrate out the heavy gauge bosons of 4321 together with the radials of the scalars Ω 1,3 and the vector-like fermion(s) at Λ 4 , then the top quark and eventually the Higgs boson at their respective masses. Performing the matching between 4321 and the SM at a common scale is a simplification, because the spectrum of the heavy modes is non-degenerate. However, thanks to the relations imposed by the SU (4) algebra Eq. (2.1), the mass differences between the gauge bosons is still expected to be relatively small and no large logarithms are generated. We choose this scale to be the mass of the leptoquark Λ 4 ≈ M U . Concerning the threshold effects at Λ 4 , we take into account the tree-level corrections due to the decoupling of the heavy scalars Ω 1,3 and neglect any subleading loop-level corrections (see Eqs. (2.22) and (2.23)).
Infrared boundary conditions. We require that all studied benchmarks respect the following phenomenological and formal constraints: • We initialize the running of the SM gauge couplings, the parameters of the Higgs potential and the top Yukawa coupling as follows [48]: • In order to address the B-physics anomalies within this framework, assuming both left-handed and right-handed couplings of the leptoquark, the key condition is [49] Using Eq. (2.1) and Eq. (2.2), the above relation is equivalent to • The lower bounds for the masses of the 4321 gauge bosons are set by LHC searches at high p T [49,50]. Since the masses of the gauge bosons are correlated, the exclusion limits extracted from searches of resonant coloron production with tt and bb final states impose the leading constraints on the whole mass spectrum: M G 4 TeV and M U , M Z 3 TeV .
• The colored radial modes O R and T R can be produced via QCD interactions at colliders. In order to avoid direct detection bounds at the LHC, their masses should be above 2 TeV (see Ref. [17] for the color triplet and Ref. [51] for the octet). For the singlet radial modes S 0,1,2 , the collider bounds are much milder since the only available ones originate from direct Higgs searches [52]. Taking the minimum lower mass bound to be 500 GeV and the mixing couplings η 1,3 to be O(0.1), they can be safely avoided. Note that the above mass bounds translate into constraints for the coupling values ρ i (Λ 4 ) and the VEVs ω 1,3 (see Eqs. (2.10-2.14)).
• Flavor constraints involving the radial T R , e.g. one-loop FCNC via the Yukawa interaction, impose an important lower bound on the ratio tan β T = ω 3 /ω 1 . As discussed 3 The fit of Ref. [49] includes the charged-current anomalies, encoded by the ratios R ( * ) Table  2.1), the neutral-current anomalies R ( * ) Table 2.2) and a plethora of low-energy constraints (see Table 3.1).
in Ref. [23], in case of the universal 4321, the radial T R together with the vector-like leptons enter box diagrams that generate contributions to the amplitudes of B s −B s and D−D mixing. These contributions are proportional to tan β −4 T and for tan β T 2 they become negligible. In the present case of the non-universal 4321 only the D −D mixing is affected, and still requires a similar suppression realized by assuming a hierarchy of VEVs i.e. ω 3 > ω 1 .
• We require that the BFB conditions Eq. (2.25) and Eq. (2.26) are satisfied at all scales.
• The mass parameters m 2 Ω 1,3 at scale Λ 4 are set to the minimum of the potential by satisfying Eq. (2.18).

General dependencies and perturbativity range
In Table 2 we give a benchmark that satisfies the IR boundary conditions and in Figure 1 we show the resulting RG evolution of the couplings. The fit to the B-physics anomalies and the collider bounds on the 4321 gauge bosons prefer a large g 4 coupling i.e. g 4 (Λ 4 ) 2.
On the other hand, the experimental bounds on the scalar masses are marginally surpassed by taking the couplings of the scalar sector within the range O(0.1 − 1). The only exception is the coupling ρ 5 which mainly fixes the mass of S 0 and therefore can be kept one order of magnitude smaller. In fact, as we argued in Section 2.2, the smallness of this parameter is technically natural. In fact, larger values for ρ 5 would imply also larger values for the rest of the ρ i parameters in order to keep the u parameter of Eq. (2.15) real i.e. u 2 > 0, and as we will see next, this severely restricts the extrapolation the 4321 model in the UV.
Ω VEVs values scalar couplings values Even though g 4 is an asymptotically free gauge coupling, the running is slow and its value remains O(1) even at high energies. Consequently, since the β-functions of the various parameters are explicitly or indirectly dependent on g 4 , the whole RG evolution is effectively driven by it. In particular, β(ρ 1 ) contains a term, which solely depends on g 4 4 and has the largest numerical prefactor in comparison to other similar terms appearing in the rest of the β-functions (see Eq. (A.13)). The reason is that the purely gauge correction to the Ω 1 quartic term (box diagram) includes all possible contractions of SU (4), contrary for example to the Ω 3 quartic for which the different contractions are split between the two terms with couplings ρ 2 and ρ 2 . As a result ρ 1 is a monotonically increasing function of the energy scale µ and exhibits the fastest running among the couplings leading to a Landau pole.
Depending on the values of g 4 and ρ 1 at Λ 4 we can estimate the energy scale where the scalar couplings grow larger than ∼ 4π and where the theory becomes non-perturbative. In Figure 2 we plot this scale as a function of the mass of the radial S 1 , which is the scalar mass that is most sensitive to ρ 1 (Λ 4 ). For the benchmark values in Table 2 and the minimum requirement of radials heavier than 500 GeV, we find that the perturbativity range of the theory for g 4 = 2 (M U ≈ 4.5 TeV) reaches almost O(10 8 TeV), while for g 4 = 3 (M U ≈ 6.5 TeV) it only reaches O(10 5 TeV). Furthermore, we see that if we require heavier radials this range drastically drops. Notice that the g 2 = 3 curve falls slower than the g 4 = 2 curve for higher ρ 1 values, because the − 45 2 g 2 4 ρ 1 term in Eq. (A.13) starts to dominate. Characteristically, the Landau poles appear at even lower energies when one further requires all the radial scalar states heavier than 2 TeV. In this case, couplings blow up at O(100 TeV). Considering the above, we argue that the limit of degenerate radials or radials with masses equal or heavier than the leptoquark, which has been used in the literature for simplicity (see for example Refs. [29,43,49,53]), is not favored. Notice also that our results are relevant to all current versions of 4321 models and consequently impose generic upper bounds for the cut-off scale of 4321 model. In Section 3.3 we will see that the indicated scale for new physics is even lower if one wishes to unify in the UV.
In this context of fundamental scalars, the question can be posed whether new degrees of freedom can appear at intermediate energies to tame this effect by generating appropriate (possibly fine-tuned) cancellations in the β-functions. This, however, seems difficult to arrange in our new physics setup (see Appendix B). The reason is that in order to cancel the 99 16 g 4 4 term in Eq. (A.13), besides having the opposite sign, a cancellation term must either have a coupling as large as g 4 or for a coupling less than 1 an unrealistically large prefactor of ∼ 500(100) for g 4 = 3 (2). In the former case, couplings such as Yukawa or scalar couplings, if large, would quickly develop Landau poles themselves. Thus, it is desirable that the large couplings in the IR be asymptotically free. However in the simplest case where they are e.g. gauge couplings, the breaking of the additional gauge symmetry would require the presence of even more scalars and scalar couplings with similar issues in their RG running. As a result, a solution within this context, if at all feasible, would probably be highly fine-tuned.

Radiative electroweak symmetry breaking
Another feature emerges when studying the RG behavior of the 4321 model above Λ 4 : EWSB triggered by radiative effects. Spontaneous EWSB requires µ 2 eff < 0, where µ 2 eff is defined in Eq. (2.23). For the 4321 model benchmark given in Table 2, µ 2 eff flows from positive in the UV to negative in the IR, radiatively triggering EWSB as shown in Figure 3. H , m 2 Ω 1 , m 2 Ω 3 , and µ 2 eff . For each, the quantity plotted is defined as X ≡ sign(X 2 )|X| (X = µ H , m Ω 1 , m Ω 3 , µ eff ). The inset plot focuses on where µ eff passes through 0, radiatively triggering EWSB near Λ 4 = 4.5 TeV. We also show the effect of individually varying each quartic coupling as |δ(λ, η, ρ) i | ≤ 5% at Λ UV . This variation propagates along the RG flow of the mass squared parameters, represented by a band enclosed by dashed lines. We also overlay the effect on µ eff from the smaller variations δ(λ, η, ρ) i in the ±1% range with a dark gray band. 4 To see how radiative EWSB occurs, consider the limit η 1,3 → 0, so that the Higgs is decoupled from Ω 1,3 . In that limit, µ 2 eff = µ 2 H , and the RG flow of µ 2 H is controlled by its β-function given in Eq. (A.4). In the absence of the Ω 1,3 parameters, the one-loop βfunction is proportional to µ 2 H and thus a µ 2 H sign flip from radiative effects alone is not possible. Instead, as µ H approaches 0, it also approaches an approximate fixed point, and β µ 2 H vanishes. Consequently, we see that the coupling of the Higgs to the additional scalars is primarily responsible for triggering EWSB radiatively, as has been pointed out in [41]. In [41], it was argued that the new scalar states must couple to the Higgs with positive mixing couplings in order to ensure µ H flows to negative values in the IR. Here since our Ω 1,3 masses are negative (for 4321 spontaneous symmetry breaking), we choose η 1,3 < 0 to reproduce the same behavior.
Note that once η 1,3 = 0, however, the relevant parameter to track over the RG flow is µ 2 eff 4 We thank Ulrich Haisch for pointing out that the perturbative determination of µ eff in Eq. (2.23) breaks down at energies much higher than the matching scale. Therefore, the black curve in Fig. 3 cannot be trusted in the utmost right part of the figure. Nevertheless, we expect it to be reliable in the region where radiative EWSB happens, i.e. in the window of the inset plot. rather than µ H , and its running is a complicated combination of radiative contributions to µ 2 H , m 2 Ω 1 , and m 2 Ω 3 , as one can infer from Eq. (2.23). In this case, it is not clear that η 1,3 are the dominant contributions to radiative EWSB. Instead, radiative EWSB will be controlled not just by the mixing couplings η 1,3 , but also by the ρ i self-couplings in the Ω 1,3 system, see Eqs. (2.5, A.6-A.7, A.13-A.19, A.21). For the benchmark used in Fig. 3, the η 3 , ρ 1,2,3,4 , ρ 2 all contribute significantly to EWSB breaking in the IR. The coupling η 1 has a smaller impact due to the fact that Ω 1 has fewer degrees of freedom than Ω 3 . Figure 3 shows that for this benchmark, there is a close coincidence of scales between Λ 4 and the scale at which µ 2 eff changes signs. This is due to the scalar mass mixing that defines µ 2 eff , equivalent to accounting for the tree-level threshold effects near the G 4321 → G SM breaking scale [44]. As explained in Section 2.2.3, in the presence of Ω 1,3 scalar fields, EWSB is a collective effect of the full set of scalar fields, and so as the RG evolution scale µ approaches Λ 4 , the onset of Ω 1,3 = 0 is contributing to the emergence of EWSB.
The steep decrease of µ 2 eff from positive to negative values near Λ 4 signals a careful arrangement of parameters required in the UV to reproduce the measured Higgs mass in the IR. The Higgs mass parameter µ eff in the IR is indeed sensitive to small deviations of parameters in the UV, as demonstrated in Figure 4. We present the response of µ 2 eff in the IR to varying the quartic couplings at Λ UV according to where η 0 i , ρ 0 i are the benchmark UV couplings values given in Table 2. One can employ one of the traditional measures [54] of fine-tuning 100 ∆ FT %, defined as where µ 2 eff (y i ) is evaluated at the IR boundary µ = m h phys . The maximum fine-tuning is at the sub-percent level ∆ FT = O(10 3 ). In every coupling except ρ 5 , 5 variations of O(1%) can yield positive µ 2 eff values in the IR, spoiling EWSB altogether, as shown in the shaded regions of Figure 4.
The variations at Λ UV propagate over the RG evolution of the theory, leading to deviations in µ 2 H and m 2 Ω 1,3 at Λ 4 , denoted by the dashed bands in Figure 3. After the spontaneous symmetry breaking at Λ 4 , deviations in the mass parameters lead to a spread in possible µ 2 eff values in the IR, represented by a gray band in Figure 3.
Returning to the limit of η 1,3 → 0 where the scalars Ω 1,3 are decoupled, we find this tuning problem is postponed to the UV. As mentioned above, in this limit we cannot rely on radiative effects to trigger EWSB, and instead µ 2 eff ≈ − (88 GeV) 2 must be arranged directly at Λ UV . A mass hierarchy between µ 2 eff and m 2 Ω 1,3 ∼ TeV 2 will endure throughout Figure 4: Here we show the sensitivity of µ 2 eff (evaluated at the IR scale µ = m h phys ) to small deviations in the Higgs-Ω 1,3 mixing couplings (left) and the Ω 1,3 quartic couplings (right) at Λ UV = 10 3 TeV. For reference, the value of µ 2 eff that reproduces the measured Higgs mass is shown by a dashed, black line. The region where µ 2 eff fails to radiatively trigger EWSB is shaded in gray.
the flow from the UV to IR. This imposed hierachy of scales at Λ UV is the trade-off for alleviating µ eff 's sensitivity to the initial conditions η i (Λ UV ), ρ i (Λ UV ). The tuning apparent in the sensitive dependence of both the Higgs mass and the radiative EWSB mechanism on couplings set at Λ UV is thus rooted in the little hierarchy problem, as discussed in [55].

Unification in the ultraviolet
In this section, we examine the implications of UV constraints on the RG evolution. The constraints we employ are inspired by the UV constructions presented in Refs. [31,34], in which the 4321 model is seen as an effective theory whose gauge group is embedded in multiple copies of the PS group. The scale Λ UV is then the scale where the symmetry group of the UV theory breaks spontaneously to the 4321 gauge group. We also choose this scale as the matching scale between the two theories (for more details see Appendix C). In those UV embeddings the scalar sectors are considerably more complex. However, they feature two common properties: i) The fields Ω 1,3 unify in the same multiplet. As derived in Appendix C.2, this implies that the following relations should hold at the matching scale Λ UV : ii) A second Higgs doublet with 4321 quantum numbers H 2 ∼ (1, 1, 2) −1/2 and mass parameter µ 2 is present. The Higgs potential is that of a Two Higgs Doublet Model (2HDM) as written in Eq. (B.3) and the mixing potential is given by Eq. (B.4). The UV constraints derived in Appendix C.2 give: The 2HDM slightly alters the procedure outlined in the beginning of the section. First of all, we assume a simplified 2HDM potential (B.3) in the Higgs basis i.e. µ 2 2 remains positive at all energies and the H 2 field does not take a VEV. However, in this case, Eq. (3.5) implies that radiative EWSB must happen in order to generate a negative µ 2 eff . Then collider bounds can be simply avoided by requiring that the masses of the heavy Higgs radials are above 500 GeV. 6 Regardless of the choice of λ i , this constraint is readily satisfied for µ 2 2 (500 GeV) (500 GeV) 2 (see Eq. (B.6)). Moreover, the λ i must satisfy the BFB conditions for the 2HDM (see Eq. (B.7)).
Here, we solve the RG equations by using the full β-functions listed in Appendix A. From the scale Λ 4 down to the scale µ 2 , we use the β-functions of the 2HDM (see Appendix B.2) and at µ 2 the degenerate heavy Higgs radials are collectively integrated out. Because µ 2 is much closer to the electroweak scale, the threshold corrections in this case are much smaller and can be omitted. On the other hand, threshold corrections at Λ UV could originate from different scales and induce important deviations to the expectations Eq. (3.4) and Eq. (3.5). We quantify this uncertainty by allowing deviations up to 10 % from the exact equalities. Despite that, due to the fast running of certain parameters, it is not a priori clear whether the RG evolution can be reconciled with even approximate conditions.
Our goal becomes then to search for sets of initial parameter values that satisfies the IR boundary conditions and via their RG evolution yield sets of values that approximately satisfy the UV boundary conditions. To obtain such sets, we perform a numerical scan over the parameter space by solving the RG equations for O(10 5 ) random sets for benchmarks g 4 = 2, 3 and Λ UV = 10, 10 2 , 10 3 , 10 4 , 10 5 TeV. In order to increase the efficiency of the scan in this high-dimensional space, we employ a Markov Chain Monte Carlo algorithm (Hastings-Metropolis) for the generation of trial sets. Additionally, we require that the UV values of the couplings are still perturbative. It is worth mentioning that the non-trivial part of the search concerns only the ρ and m Ω 1,3 parameters. This is because in the absence of concrete IR conditions for η, λ 2,3 , µ 2 and a loose dependence on g 4 , one can always adjust the IR values in order to find appropriate UV ones. 7 In Figure 5 we present the distribution of the obtained initial values of a representative pair of parameters, namely ρ 3 and ρ 4 . The first observation is that we find no acceptable sets for unification scales as low as 10 TeV and very few sets for scales as high as 10 5 TeV (and only in the case of g 4 = 2). This is explained by the fact that the different ρ parameters exhibit different running behaviors and thus convergence to the UV conditions is possible only in the intermediate scales 10 3 − 10 4 TeV. For instance, the relation ρ 1 = ρ 2 + ρ 2 cannot be satisfied at higher energies due to the faster running of ρ 1 and the fact that ρ 2 is monotonically decreasing as opposed to the other two. In total, we found more sets for 6 Notice that this bound is sufficient also for other types of 2HDMs, e.g. where flavor constraints are relevant [56]. 7 We note that in multi-site UV constructions, like the PS 3 model, the Ω and the Higgs fields may unify into their respective multiplets at different scales. Because, as mentioned, the problem of matching to the conditions Eq. (3.5) is much less constrained, one could also achieve this without significant deviations from our discussion.
smaller g 4 because, as already discussed, the overall running in this case is slower and thus more suitable for achieving unification.
On the other hand, for lower unification scales the RG evolution spans a shorter energy range and the running has to be accelerated in order to meet the UV conditions. This is the reason why in Figure 5 we see that the initial values of the parameters increase as Λ UV decreases so that the terms in the β-functions proportional to the couplings themselves are enhanced. Altogether, the attempt to unify the couplings requires Λ UV approximately three orders of magnitude lower than the cut-off scales discussed in Section 3.1 for the generic case. We conclude the discussion by illustrating in Figure 6 the unification of the masses and couplings using as an example the RG evolution generated by one of the obtained sets.  The y-axis, scalar masses, is defined as in Figure 3.

Conclusions
In this work, we performed an analysis of the RG evolution of the 4321 model, which revealed several interesting features inherent to the scalar sector of the model. We observed that the large value of g 4 plays an important role in the overall behavior of the running and it can lead to the appearance of Landau poles even though the gauge coupling itself is asymptotically free. Phenomenological constraints on the Ω 1,3 radial modes prevent their quartic couplings from assuming arbitrarily small IR values, leading these quartic couplings to blow up, and compromising the validity of the 4321 model, at energies far below the GUT scale. Additionally, we found that it is possible to satisfy relations from motivated UV embeddings within the energy range O(10 2 − 10 4 ) TeV.
Intriguingly, the emergence of Landau poles seems to be intrinsic to a large class of the models proposed as combined explanations for both charged-and neutral-current B-physics anomalies. Depending on the model, the problem may appear in different sectors of the theory but typically involves the presence of fundamental scalars (see Refs. [23,57,58] for examples of Landau poles in the Yukawa couplings). As discussed, it is highly non-trivial to address this problem at the level of β-functions by introducing additional fields. This might indicate that the 4321 gauge symmetry is broken not by fundamental scalars, but rather à la technicolor by the condensate of a strongly coupled sector as in Ref. [27,33]. Another interesting possibility is that the leptoquark arises directly as a composite resonance of a new strongly coupled sector as in Refs. [59,60]. For example, the model in Ref. [60] with composite scalar leptoquarks and Higgs is free of low-energy Landau poles.
Upon studying the RG flow of the mass squared parameters in the scalar sector, we found that the 4321 model can accommodate radiative EWSB due to the additional scalar fields coupled to the Higgs. While radiative EWSB requires positive contributions to the βfunction generated by the H-Ω 1,3 mixing, the IR value of the parameter that controls EWSB, µ 2 eff , relies on the whole system of scalar quartic couplings, mixing couplings, and Ω 1,3 mass parameters. We quantified the sensitivity of the Higgs mass to the UV values of the parameters. The resulting fine-tuning reflected the little hierarchy problem associated with the new TeV scalars.
Ultimately, if the B-physics anomalies are confirmed as clear signals of new physics, leptoquarks arise as the most plausible explanation. Given that all leptoquark models contain additional scalar particles, we stress the importance of carefully examining their sectors, and in particular, the RG evolution of the scalar couplings. As we have seen in the case of 4321, new features can emerge upon analysis of radiative effects, constraining the possibilities for UV constructions. A β-functions of the 4321 model The β-functions for the gauge couplings at one-loop level are given by the master formula: The first term contains the quadratic Casimir of the gauge boson representation C 2 (G), the second and the third the sums over the Dynkin indices for the representations of all the fermions S 2 (R f ) (κ f = 1/2 or κ f = 1 for 2-or 4-component fermions) and the scalars S 2 (R s ), respectively. Taking into account the full particle content of the 4321 model the resulting values for b i are where n VL is the number of vector-like fermions.
The β-functions for the couplings of the potential in Eq. (2.4) are given at one-loop order:

B.1 Yukawa sector
Neglecting the Yukawa couplings of the light generations, the Yukawa sector of the 4321 model reads withH = iσ 2 H * and where we used the freedom of rotating the multiplet Ψ L , χ L to remove the mixing termΨ L χ R . The biggest couplings are y t ∼ 1 and y + ∼ 1/4, while the other couplings are expected to be smaller, i.e. λ q, ∼ 10 −1 , y b ∼ 10 −2 and y − ∼ 10 −3 from phenomenological considerations [49]. The influence of these terms on the scalar quartics β-functions were computed and did not yield any relevant contribution. For example, the contribution should be compared to the g 4 4 term in Eq. (A.13).

B.2 Two higgs doublet model
and the mixing potential in Eq. (2.7) by For the Yukawa couplings, the same terms with H in Eq. (B.1) can be written with H 2 .
Radial modes. In the Higgs basis where only H takes a VEV, the two doublets can be parametrized as where η ± W , η Z become the W ± and Z radial modes, respectively, and h is the physical Higgs boson. The scalar radial degrees of freedom h, h ± , h R and h I aquire the following masses [62,63] BFB conditions. For the potential (B.3), necessary and sufficient conditions to ensure boundedness from below were found to be [64] λ > 0 , For the mixed part of the potential V ΩH , similarly to Eq. (2.26), the derived conditions are β-functions. The β-functions above the 4321 breaking scale are presented in Appendix A with the gray terms to be included. Another small difference is that the coefficient in front of y 3 t in Eq. (A.23) becomes 6 instead of 11/2.
Below 4321 breaking scale, we can run in the 2HDM until the mass threshold of the heavier Higgs. The most general Yukawa sector is This field could also be used to split the leptons and quarks masses directly, as with H 15 , by introducing higher dimensional operators, e.g. 1 ΛΨ L Ω 15H Ψ + R . In the broken phase of 4321, this field decomposes into Both H 15 and Ω 15 have a more involved scalar potential which could give rise to the appearance of new Landau poles similarly to the Ω 1,3 potential. Moreover their mixing coupling with the Ω 1,3 sector will give contribution to the quartic couplings studied in the main analysis with the same sign as the g 4 4 contribution in Eqs. (A.13-A.17). Their effect will therefore only strengthen the conclusions about the appearance of Landau poles.
Another important effect of extending the scalar content of the theory is its influence on the running of the gauge couplings as can be seen in Eq. (A.2). At one loop, the condition for asymptotic freedom is b i > 0. The β-functions coefficients b i for i = 1, 2, 3, 4 with the minimal scalar content are presented in Eq. (A.3). Since the IR value of g 4 is large and drives most scalar quartic beta-functions, the faster it decreases in the UV, the later the Landau poles arise. Adding more scalars is equivalent to decreasing b 4 and slowing down its approach to asymptotic freedom. The influence of each scalar on the b i coefficients is presented in Table 3. We can conclude that adding all these scalars to the theory is not enough to directly ruin the one-loop asymptotic freedom condition for gauge couplings, but they would accelerate the appearance of the Landau poles in the scalar sector as discussed in Section 3.1.  Table 3: Dynkin representation multiplicity factor of each scalar.

C Constraints from UV embedding
Assuming a PS 3 origin of the 4321 model as in [31] yields specific relations between the coefficient of the scalar potential and the Yukawa couplings. By considering the representations of the fields in this UV setup, their decompositions under the subgroup L ∼ (4, 2, 1) 3 → Ψ L ∼ (4, 2, 1) 0 . (C.2)

C.1 Yukawa sector
We consider a Higgs bi-doublet field that couples mostly to the third generation of fermions. Neglecting light generations coupling, two Yukawas terms can be written with Φ c = σ 2 Φ * σ 2 . The complex bi-doublet Φ can be written as where H and H 2 are SU (2) L doublets andH (2) = iσ 2 H * (2) . Since Ψ Generally, the Higgses will get the aligned VEVs (C. 5) and give rise to the masses where we defined v 2 = v 2 1 + v 2 2 and tan α =

C.2 Scalar sector
Higgs potential (2HDM). The most general potential for a bi-doublet Higgs is [47] V UV  Ω potential. The UV relations of the Ω potential coefficients would require considering all degrees of freedom from the decomposition in Eq. (C.2). The next step would be to decouple triplets of SU (2) L , so that Ω only takes a VEV along the SU (2) L -preserving direction. Neglecting the contribution from integrating out Ω 1,3 , we write a multiplet where these degrees of freedom are omitted from the start, as well as the SU (2) L and SU (2) 3L indices, where Ω 3 and Ω 1 are both anti-fundamental of SU (4) 3 and the vector is an SU (4) multiplet.