Perturbative unitarity constraints on generic Yukawa interactions

We study perturbative unitarity constraints on generic Yukawa interactions where the involved fields have arbitrary quantum numbers under an $\prod_i SU(N_i) \otimes U(1)$ group. We derive compact expressions for the bounds on the Yukawa couplings for the cases where the fields transform under the trivial, fundamental or adjoint representation of the various $SU(N)$ factors. We apply our results to specific models formulated to explain the anomalous measurements of $(g-2)_\mu$ and of the charged- and neutral-current decays of the $B$ mesons. We show that, while these models can generally still explain the observed experimental values, the required Yukawa couplings are pushed at the edge of the perturbative regime.


Introduction
Yukawa interactions are, together with gauge and scalar self-interactions, the building blocks of renormalizable theories. In the Standard Model (SM) they are a crucial ingredient as they describe the interactions between quarks and leptons and the Higgs boson, with the latter being ultimately responsible for the generation of fermion masses after electroweak symmetry breaking (EWSB).
They are also ubiquitous in new physics (NP) theories that try to address the shortcomings of the SM. For example, they appear in theories that generate neutrino masses both at tree-level, as for the case of the well know seesaw mechanism [1][2][3][4][5][6][7][8][9], as well as at higher orders [10][11][12][13] and in models where a fermion Dark Matter (DM) candidate is connected to the SM through a scalar portal, see e.g. [14] for a review. Interestingly, the existence of new scalar bosons or some beyond the SM (BSM) fermions that possess Yukawa interactions with the SM could solve some anomalies reported in the recent years in low energy data. This is, for example, the case of the measurement of the anomalous magnetic moment of the muon (g − 2) µ , for which the recent measurement by the E989 experiment at Fermilab [15], which is in agreement with the previous BNL E821 result [16], implies a ∼ 4.2σ discrepancy with respect to the SM prediction [17], although a recent lattice calculation seems in agreement with it [18]. It is also the case of other long-standing anomalies in semileptonic decays of B mesons both in charged- [19][20][21][22][23][24][25][26][27] and neutral-current [28][29][30] decays, usually dubbed R D ( * ) and R K ( * ) , that can be accounted for by various models involving additional Yukawa sectors. Some or all of these anomalies can be solved by postulating the existence of leptoquarks (LQ) , i.e. new colored states which connect quark and leptons, or by extending the SM with new heavier scalars and vector-like fermions [79][80][81][82][83][84][85][86][87][88][89][90][91][92]. Obviously, new Yukawa interactions are constrained by a large variety of experimental searches, ranging from direct production of new on-shell degrees of freedom at high energy colliders to low energy precision measurements. The results of these analyses are generally expressed as limits on combinations of couplings and masses, and the resulting bounds strongly depend on the specific structure of the NP realization and on the experimental settings.
There exists, however, and old tool of theoretical physics, namely perturbative unitarity (PU), that can be used to set an upper limit on the magnitude of the couplings, above which the perturbative expansion is expected to break down. Most famously this tool, that we review in Sec. 2.1, has been applied to set an upper bound on the Higgs boson mass [93][94][95][96] and on the masses of quarks and leptons participating in weak interactions [97,98] 1 if weak interactions were to remain weak at all energies. It has then been widely used in the literature to assess the range of validity of both renormalizable and effective operators [101][102][103][104][105][106][107][108][109][110][111][112][113][114]. In this work we consider the problem in more generality and we answer the following question: Given a Yukawa interaction between a scalar and two fermions with generic quantum number under a group G = i SU (N i ) ⊗ U (1), what is the maximum allowed value for the coupling with the requirement of PU?
To answer this question we consider the most general form of Yukawa-type interactions and all possible 2 → 2 tree-level scatterings in the high-energy limit. We obtain compact expressions for the upper limit on the value of the Yukawa coupling up to which perturbation theory could be trusted, and highlight their dependence on the various fields' quantum numbers under G.
More specifically, we start by computing all the necessary ingredients for building the partial wave scattering matrix, namely the Lorentz parts of the scattering amplitudes and the group structure factors entering the amplitudes themselves, in a set of phenomenologically relevant toy models where the various fields are only charged under a single SU (N ) factor. We firstly show how the SU (N ) group structure of the interaction can lead to en enhancement of the scattering amplitudes and thus to a tightening of the partial wave unitarity bounds, while the role of the U (1) charge is to enforce a selection rule that makes some amplitudes vanish. We then use these toy models as building blocks for more complicated theories, where the various fields are charged under multiple SU (N i ) factors, giving as a working example the case of the SM quark Yukawa sector. We then apply our results to different NP models which solve, by introducing a new Yukawa sector, the aforementioned anomalies in (g − 2) µ and/or in the semileptonic decay of B meson, showing that while the proposed theories can generally still provide an explanation to these measurements, their model parameters are stretched close to the limit of our-tree level unitarity bound criteria.
Altogether the results presented in this work are of practical use and can be used to analyze scenarios beyond the examples presented in the text. While we restrict only to a limited number of irreducible SU (N ) representations under which the various fields can transform (singlet, fundamental and adjoint), we believe that our computations furnish the necessary ingredients to study a large set of NP theories with additional Yukawa interactions.
The paper in organized as follows. In Sec. 2.1 we review the tool of PU and clarify the physical interpretation of the inferred bounds, while in Sec. 2.2 we discuss the general properties of Yukawa interactions relevant for the study of PU. In Sec. 3 we introduce the toy models and discuss how they can be used to construct the most general partial wave scattering matrix. Then in Sec. 4 we study the partial wave unitarity bounds for the first type of toy models which have a Dirac type structure for the Yukawa interaction. In Sec. 5 we apply our formalism to the case of the SM quark Yukawa sector, also highlighting the role that multiplicity due to a flavor structure can have in the determination of the bound. In Sec. 6 we then discuss a second class of toy models, which present Majorana type Yukawa interactions. Then in Sec. 7 we show some phenomenological applications, finally concluding in Sec. 8. We also add few relevant appendices. In App. A we list our conventions for the calculation of the partial wave matrix, while in App. B we report the results for other Dirac type theories not included for brevity in the main text.

The tool of perturbative unitarity
Perturbation theory is a powerful tool to provide approximate solutions to physical problems. Within this approach the relevant result is expressed in terms of a power series of some small parameter . In general, however, it is not easy to determine if some specific numerical value of allows for a good approximate solution of the problem under consideration. The aspect we want to face in this Section is to identify a reasonable criterium to state whether a parameter is too large to be treated in perturbation theory.
A level zero criterium for considering a parameter as perturbative is to ask the expansion parameter entering the β functions = g 2 (4π) 2 , to be small. The requirement < 1 implies g < 4π, which is the maximum value frequently adopted in the literature. This condition can however be improved by analyzing the problem in more detail. Consider an abelian gauge theory with coupling g and N f copies of matter fields charged under the local symmetry. The scaling of the contribution to the one-loop two point function for the gauge field is The requirement of small expansion parameter then suggest a more refined version of the naive criterium, since a small value of g is in fact not enough if N f becomes very large. One can then require g < 4π √ N f , where it is clear that the multiplicity of the matter fields has to be taken into account for a more refined version of the perturbative criterium. One can however do even better than the former proposal. It is possible to use results that hold beyond perturbation theory to motivate a more stringent criterium based on partial wave PU.
The key point of our analysis are the so-called partial waves, i.e. the scattering amplitudes with fixed total angular momentum J. In the case of 2 → 2 scatterings in the high-energy massless limit they are defined as [115] Here θ is the polar scattering angle in the center of mass frame and √ s the center of mass energy, d µ i µ f (θ) are the small Wigner d−functions where µ i = λ i 1 − λ i 2 and µ f = λ f 1 − λ f 2 are defined in terms of the helicities of the initial and final states, and (2π) 4 δ (4) , with S the S-matrix, that defines the scattering amplitude. We report in App. A.1 the definition and the explicit expressions of the small Wigner d−functions used throughout our analysis. The unitarity condition on the S−matrix, where the sum runs over all the intermediate states h. By focusing on elastic channels i = f and restricting the sum over h only to 2-particle states one obtains the condition This last equation defines a circle in the complex plane, the Argand circle, inside which the amplitude must lie at all orders in perturbation theory Since for 2 → 2 high-energy scatterings the tree-level elastic amplitudes are real 2 , this suggests the following unitarity bound While the factor 1/2 is somewhat arbitrary, it gives a reasonable indication of the range of validity of the perturbative expansion, since a tree-level value which saturates Eq. (2.5) needs at least a higher-order correction of ∼ 40% in order to re-enter the Argand circle, thus signalling the breakdown of the perturbative expansion itself. In order to extract the best PU bound, one then needs to identify the optimal elastic channel and this corresponds to diagonalizing the partial wave scattering matrix a J f i . The most stringent limit will be then set by the largest, in absolute value, eigenvalue.

The structure of the Yukawa interaction
In this section we specify the class of models we are interested in and the assumptions we make in our analysis.
First of all we assume that mass terms are negligible and all the computation are performed in the high energy regime. In this limit we can consider massless fields as the most natural degrees of freedom. Regardless of the symmetry structure of the interactions, we can always (re-)write a generic Yukawa interaction between a set of N φ real scalar fields φ α and of N ψ Weyl fermion fields ψ i L in the following way: where ψ c L,i = Cψ T L,i , C = iγ 2 γ 0 is the charge conjugation matrix, Y αij = Y αji and the index α runs from 1 to N φ , while i, j = 1, . . . N ψ . Notice that this form is the most general one. For example, a complex scalar field can be always expressed in terms of two real fields with specific restrictions on the phases of the coupling Y αij , in a similar way symmetry properties of the Yukawa interactions are manifest through the presence of null elements or by specific relations among them. We are going to clarify these aspects in what follows with explicit toy models.
Our task is now to compute the partial wave matrices a J f i of Eq. (2.1) and extract their eigenvalues. We start discussing the helicity and Lorentz structure of the scattering amplitudes. To this end we firstly need to compute the T f i ( √ s, cos θ) amplitudes between the initial and final states. It is useful to write the scattering matrix T in the following form where we have indicated the total helicities of the initial and final states, µ i,f , as well as the helicities of the single particles involved in the scattering. Each of the amplitudes entering the various blocks of the T f i matrix are themselves matrices, whose dimensionalities depend on the number of fermions and scalars of a given theory. Notice that the matrix of Eq. (2.7) has many zero entries, which correspond to the empty blocks. In particular the T ±0±± , T ±0+− and T ±000 amplitudes vanish because of total angular momentum conservation, while T 0000 is present if one adds also a potential for the scalar fields. Finally, the other null amplitudes are strictly zero only in the massless limit we are considering. Since we are working in the high-energy massless limit it is useful to compute the non vanishing amplitudes by working with helicity eigenstates, following the conventions of Jacob and Wick [115]. We refer the reader to App. A.2 for the details on the choice of the spinor helicity basis. The remaining non vanishing amplitudes can be computed explicitly starting from the interaction in Eq. (2.6). All the non vanishing amplitudes are reported in Eq. (A.10) The usefulness of Eq. (2.7) is that, when projecting the amplitudes onto the J−th partial wave via Eq. (2.1), only a subset of the non-zero blocks survives. In particular the channels with µ i = µ f = 0 have a non-zero projection only on J = 0, since the relevant amplitudes do not depend on θ, see Eq. (A.10). Again because of total angular momentum conservation channels with µ i,f = ± 1 2 project only on half-integers values of J. Hence, for integer J > 0 only the channels T +−+− and T 00+− contribute and the scattering matrix is effectively separated in three different blocks, allowing us to consider each of them independently when studying partial waves with different values of J. Interestingly, we will show that the stronger bound might arise from a partial wave different from J = 0. In particular in Sec. 4 we will see that the tighter limit can come from the analysis of J = 1 2 or J = 1, while higher partial waves give a weaker bound. Table 1. Dirac type and Majorana type theories described by the Lagrangians of Eq. (3.1) and Eq. (3.2), respectively. If S is either a singlet or transforms in the adjoint representation of SU (N ) and also has a vanishing U (1) charge, then it is a real scalar field. Model 5 of the first class is only present in the case of SU (3), while model 2 of the second class only in the case of SU (2).

Toy models and scattering amplitudes
The Lagrangian of Eq. (2.6) describes the Yukawa interaction between any set of scalar and fermion fields, where the entries of the Yukawa matrix Y αij are at this level completely generic. By assigning definite quantum numbers under a group G = i SU (N i ) ⊗ U (1) to the scalar and fermion fields involved in the interaction, the Yukawa matrix, and consequently the scattering matrix, acquires a well definite structure. In this Section we consider two types of toy models which we use as building blocks for the study of more general Yukawa theories. More specifically we consider two theories, described by the following Lagrangians: − L Dirac = ySχη + h.c. , (3.1) and where χ and η are left-handed and right-handed fermion fields respectively, and S a scalar field that can be either complex or real. We dub these theories as Dirac type and Majorana type respectively. To begin our study, we start by assuming that all fields are charged under a single SU (N ) factor and show later in Sec. 5 how the case of multiple SU (N i ) charges can be dealt with. For both theories the parameter y can be chosen to be real without loss of generality by a proper field redefinition. For concreteness we restrict our study to the case where all the fields transform in the trivial, fundamental or adjoint representation of SU (N ) and we allow them to have arbitrary U (1) charges. For simplicity we also consider theories where at most one field transforms in the adjoint SU (N ) representation. Under these assumptions the various models that can be written are reported in Tab. 1. Note that in the Dirac type class, model number 5 can only be written in the case of SU (3) while in the Majorana type class, model number 2 can only be written in the case of SU (2). A comment regarding other possible models in the Majorana type class is in order: • For the case of SU (2) one can write a gauge invariant interaction with χ ∼ q and S ∼ 1 2q which however identically vanishes sinceχεχ c = 0, where ε is the SU (2) totally antisymmetric tensor. One can restore this interaction by charging the field χ under a second SU (N ) factor, or by considering different flavors for χ, since in this case one can antisymmetrize in the additional gauge and/or flavor index.
• For the same reason, in the case of SU (3) also the model where χ ∼ q and S ∼ 2q vanishes if the states are not charged under another SU (N ) group or flavor multiplicity is not added.
We will comment on these possibilities when considering in more detail the Majorana type class of models, presenting in Sec. 6.3 a phenomenologically relevant example in which the involved fields are charged under more than one SU (N i ) factor.
For any given toy model, the task is to to build the a J f i partial wave matrices and compute their largest eigenvalues. This can be done mechanically by brute force by building the Yukawa matrix Y αij entering Eq. (2.6), computing all the amplitudes 3 and building the a J f i matrices explicitly through Eq. (2.7). Although straightforward, this process turns out to be highly inefficient, due to the rapid increase in the a J f i matrix dimension when considering SU (N ) factors with large N . As an example for the third theory of the Dirac type with a complex scalar field, the transition matrix has dimension 406 for N = 3 and 2346 for N = 5. When considering the possibility of fields charged under more than one SU (N i ) factor the dimensionality of the transition matrix dramatically increases, making also the numerical calculation inefficient.
The situation drastically simplifies if one realizes that when considering any 2 → 2 scattering, each amplitude can be decomposed into a Lorentz part which depends only on the spin and helicity of the involved fields, and a group-theoretical part that depends on their SU (N ) quantum numbers, while the role of the U (1) charge is to enforce a selection rule that will make some amplitudes vanish. More concretely any 2 → 2 scattering amplitude among particles i 1,2 and f 1,2 with helicities λ i 1,2 and λ f 1,2 can be written schematically as T where T is the Lorentz part of the scattering amplitude and F m,r f 1 f 2 i 1 i 2 (N ) is a function that contains the group part coefficient for the scattering through the Mandelstam m−channel in the SU (N ) r irreducible representation that can be built from the initial and final state particles, while d r stands for the dimensionality of r. The direct sum runs over all the irreducible representation through which a scattering can proceed. One of the necessary ingredients are thus the T functions for the two theories of Eq. (3.1) and Eq. (3.2), which can be computed from the Lagrangian of Eq. (2.6). We report them, normalized by the common y 2 factor, in Tab. 2 and Tab. 3 for both the real and complex scalar S case 4 . Here we clearly see the role played by the U (1) factor. As 3 Here one has to consider the relevant factors of 1/ √ 2 for identical particles, which can occur only for two-fermion states when µi and/or µ f = 0 and for two-scalar states. 4 As mentioned in Sec. 2.2 we find now convenient, in the case of a complex scalar field, to directly work with its complex components instead of the real ones, since as we will see the presence of the U (1) symmetry allows us to simplify the scattering structure. Table 2. Lorentz part of the amplitudes for the models of the first class divided by y 2 . In the complex scalar case also the 0 * 0 * + − amplitudes are zero. In the helicity amplitudes the notation 0 * indicates the scattering involving the conjugate of the scalar S.

Dirac type models
an example, when S is a complex scalar field the amplitude χη →χη in T ++−− is zero, because this process violates the conservation of the U (1) charge. Analogous selection rules appear in other scattering channels. In order to fix the idea let us make an explicit example and consider the first theory of the Dirac type class, where both χ and η transform under the fundamental representation of SU (N ) and S is a scalar singlet. The scatteringχη →χη which proceeds through the + + ++ helicity channel has only an s−channel contribution, with an amplitude which is proportional to −y 2 , see Tab. 2. Since ⊗ = 1 + Adj, this scattering can only proceed through the singlet and the adjoint channels. Thus the scattering amplitude, by applying Eq. (3.3), reads 5 Table 3. Lorentz part of the amplitudes for the models of the second class divided by y 2 . In the complex scalar case also the 0 * 0 * + − amplitudes are zero. In the helicity amplitudes the notation 0 * indicates the scattering involving the conjugate of the scalar S.
With the same procedure one can build the amplitudes among irreducible representations for all the possible scatterings of the theory. In each of the separated subsectors that we have identified (J = 0, half-integer J and integer J > 0, see Eq. (2.7)), the matrix can be decomposed into scattering blocks among the various irreducible SU (N ) representations. Barring the convolution with the Wigner d−functions and the integration over the angular variable θ, finding the eigenvalues of the partial wave matrix, and thus extracting the partial wave unitarity bound, is then a trivial task. One only needs to compute the F m,r f 1 f 2 i 1 i 2 (N ) factors. The advantage of this procedure with respect to the mechanical brute force one previously described is clear: being the group factor proportional to the identity in group space, one needs in practice to consider for each representation only one scattering among the d r ones, since all of them will give the same result 6 .

Dirac type theories
In the previous Section we have described the general strategy for computing the a J f i partial wave matrices and presented the Lorentz part of the amplitudes that are needed to compute them. In this Section we compute, for the various models of the Dirac type class presented in Tab. 1, the F m,r f 1 f 2 i 1 i 2 (N ) group factors. For brevity of presentation we report in the main text only the results for the first two type of models, while we defer to App. B for the remaining ones.

First model
In this model S is an SU (N ) singlet, real if q = q , while η and χ transform under the fundamental representation of the SU (N ) group. By choosing as basis 7 where a runs from 1 to N . The Yukawa matrix of Eq. (2.6) thus reads where α = 1 in the real scalar case and α = 1, 2 in the complex scalar one. As previously stressed we can separately consider the sectors with J = 0, half-integer and integer J > 0, since no amplitude has a non zero projection on more than one sector. Let us start by considering J = 0 with a real scalar S. The two particle states with ++ helicities decompose as while the states with −− helicities are their conjugates. Here above we indicate with S and AS the totally symmetric and antisymmetric irreducible representations that arise from the tensor decomposition ⊗ = S ⊕ AS. Note, however, that the antisymmetric combination of two identical fermions identically vanishes for even J, see e.g. [115]. In order to compute the scattering amplitudes we need to explicitly write the two-particle states. We define them as where ψ = χ and/or η, T A are the SU (N ) generators, and T A S are N (N + 1)/2 symmetric matrices that we choose to be the symmetric SU (N ) generators, with the addition of 1 N √ 2N , which is normalized to preserve the canonical trace normalization Note that in this case the symmetric combination is always built through two identical states: this is the reason of the extra factor 1/ √ 2 in the symmetric two particle state with respect to the adjoint one. By direct computation one obtains the following non zero group factor amplitudes relevant for the scattering in J = 0, see again Eq. (2.7), where the − − −− and − − ++ amplitudes are equivalent since they are obtained by conjugation. Note that among the vanishing + + ++ group factors we have that, for example, the one in the adjoint channel is zero because S is a scalar singlet and the amplitude thus turns out to be proportional to the trace of the SU (N ) generators, while the one in the symmetric channel is zero since the theory does not mediate processes such as ηη → ηη. From the explicit form of the group factors it is clear that it's the singlet channel that manifests an enhancement of the scattering amplitude due to the SU (N ) group structure. By using the general expression of Eq. (3.3), the J = 0 partial wave in the singlet channel written in the (χη, χη) basis reads whose largest eigenvalue in absolute value is y 2 16π (2N + 1). The eigenvalues relative to the scatterings in the other SU (N ) irreducible representations are all smaller and thus the PU condition of Eq. (2.5) leads to the bound If the scalar is complex the conservation of the U (1) charge forbids scattering in the ±±∓∓ channel. The matrix of Eq. (4.6) becomes thus diagonal with eigenvalues y 2 16π N and the bound now reads which is weaker than in the real scalar case. For half-integer J the two-particle states with +0 helicities decompose for real scalar S as while the states with −0 helicities are their conjugates. The group factors read which are equivalent to their conjugates. Since there is no scattering between the +0 and −0 states in the massless limit, we can consider only the +0+0 scattering channel, while the −0 − 0 will be its conjugate. For example for the former the scattering in the fundamental channel for J = 1/2 reads thus leading to the bound One can understand that there is no multiplicity factor due to the SU (N ) group structure since for all the diagrams, both in the s− and u−channel, the SU (N ) index is never contracted between initial or final states, but is instead conserved between them. If S is a complex scalar the U (1) charge conservation enforces to treat separately the scattering of two distinct Mandelstam channels. As an example both the ηS → ηS and the ηS * → ηS * scattering proceed through the fundamental representation, but the former via an s−channel diagram, while the latter via u−channel one. The different angular function of the two amplitudes makes the one in the u−channel dominate. For the scattering ηS * → ηS * one obtains which leads to the bound For both the real and complex scalar the bounds in the J = 1/2 sector are weaker than in the J = 0 one. Finally we can consider the relevant scatterings for integer J > 0, starting again with the case of a real scalar field. Here the two particle states among fermions decompose as while, the scalar field being an SU (N ) singlet, only the trivial representation exists for 00. Bose symmetry however forbids the scattering in the 00 + − channel for odd J [115], while the relevant group factors in the + − +− channel are 8 (4.14) For the lowest partial wave J = 1 all the eigenvalues of the partial wave matrix are ± y 2 32π , therefore the bound is simply In the case of a complex scalar field instead one has additional non vanishing processes as SS * →χχ and SS * → ηη that can however proceed only via the singlet channel, since the scalar belongs to the trivial SU (N ) representation. The scattering among singlet is thus modified with respect to the real scalar case. The group structure for these scatterings is and explicitly one has in the (χχ, ηη, SS * ) basis where • represents the Hadamard product 9 . Here we have split the Lorentz and group part of the amplitudes to highlight that the different helicity channels are associated to different group coefficients. The largest eigenvalue in absolute value is y 2 64π (1 + √ 1 + 16N ) and the bound thus reads We then report in Fig. 1 the bounds on y obtained in the J = 0, 1/2 and 1 partial waves.
In particular we see that in the case of a real scalar field the strongest bound is obtained in the J = 0 channel. On the other hand for a complex scalar field the strongest bound is obtained through the analysis of the scattering in the J = 1 channel, while the one in the J = 0 channel dominates only in the case of an abelian theory. This is a non trivial result, highlighting the role that higher partial waves can have in deriving a PU bound.
In this model the scalar transforms in the fundamental representation of SU (N ), and is thus always a complex field. By fixing the basis where a runs from 1 to N , one has for the Yukawa matrix of Eq. (2.6) with i, j = 1, . . . , N + 1. We start again by considering the J = 0 partial wave. Here the two particle states decompose as and the −− states are the conjugates. Again, the antisymmetric combination identically vanishes for even J. The two-particle states can be written in analogy with Eq. (4.4), where clearly the singlet combination is now the trivial state |ηη , to which we add the state in the antifundamental as |χη a = |χ a η . In this theory all the scatterings in the ± ± ∓∓ channels vanish together with the ones that proceed through the singlet and symmetric channels in ± ± ±±. The only non vanishing amplitudes are the ones in the (anti)fundamental channel with a group factor that reads In this case the partial wave matrix is already diagonal and after the trivial integration over the angular variable one obtains the bound Moving onto the J = 1/2 channel we can consider the +0 + 0 scattering. Here we have and the group factors for the non zero amplitudes are with again the −0 − 0 being the conjugates. We can consider the singlet channel, which exhibits a N multiplicity factor. In the +0 + 0 sector one has This is the largest eigenvalue for N > 1 while for N = 1 it's the scattering in the (anti)symmetric channels that dominates due to the u−channel amplitude which scales as T +0+0 u ∼ 1 cos θ 2 and has eigenvalue ±2y 2 . Altogether we obtain (4.26) Finally for J = 1 the relevant two-particle states decompose as where however now it's the symmetric combinations of the two identical scalars that vanishes identically for J = 1 [115]. The group factors for the non-zero amplitudes are Note that the 00 * + − scattering in the singlet channel has two contributions, coming from the two possible ways of making an SU (N ) singlet from the two fermions. The strongest bound turns out again to be the one arising from the scattering among singlets. Explicitly one has in the (χχ, ηη, SS * ) basis where we again have split for convenience the Lorentz and group structure of the amplitude. The eigenvalues of this matrix have a complicated form for generic N and we thus show the numerical results in Fig. 2. There we see that for N < 4 it's the J = 1 partial wave that enforces the strongest bound while for N ≥ 5 is the J = 1 2 one. As for the case of the previous toy model we see that the stronger limit can arise from partial waves different from J = 0. We also note that when N = 1 we recover, with no ambiguity, the same bounds obtained for the first toy model in the case of a complex scalar field.

The case of the Standard Model Yukawa sector
The two toy models discussed in the previous Section can be used as building blocks through which it is possible to study more involved theories where, e.g., the fields are charged under multiple SU (N i ) factors and/or where more than one generation of fields with the same quantum numbers is present. We highlight this by discussing in detail the case of the SM Yukawa couplings, focusing on the down type quark sector, for which the Yukawa Lagrangian is

Multiple SU (N i ) factors
We start by discussing the role played by multiple SU (N i ) factors under which the various fields can be charged. To this end we consider a single generation of SM fermions, y 11 d = y d . The rule of Eq. (3.3) is readily generalized, by considering that now each group factor coefficient F is the product of the various group coefficient factors for the different SU (N ) groups and the dimension of the identity matrix is the product of the dimensions of the considered irreducible representations for each SU (N ) factor. It is again instructive to work out the most important scattering amplitudes for the case of J = 0, half-integer J and integer J > 0. We start with J = 0. Since the scalar is complex the only non-zero amplitudes are in the ±±±± scattering channel. Working in the (qd, qd) basis one has that the scattering proceeds through the (anti)fundamental channel for what concerns SU (2) L with a group factor proportional to the identity, see Eq. (4.22). As regarding SU (3) c one can again consider the scattering in the singlet channel, which exhibits a group factor enhancement, see Eq. (4.5). The J = 0 partial wave explicitly reads Here the presence of two SU (N ) groups simply increases the eigenvalue multiplicity, given that SU (2) L group factor is proportional to the identity and the scattering matrix is already diagonal. The bound in this case is simply the one of Eq. (4.8) For J = 1/2 we can consider the +0+0 sector and it's convenient to consider theqH ↔qH in the antifundamental channel for SU (3) and singlet channel for SU (2). The group factors can be read from Eq. (4.9) and Eq. (4.24) and are the identity for SU (3) and N 2 = 2 for SU (2). Also in this case the presence of two SU (N ) factors simply increases the eigenvalue multiplicity. The partial wave reads The situation is more involved for J = 1. Focusing on the singlet channel scattering for both the SU (3) c and SU (2) L groups one has the group coefficients of Eq. (4.14), Eq. (4.16) and Eq. (4.27). Explicitly then in the (q L q L , d RdR , HH * ) basis one obtains Here we see the non trivial interplay between the two SU (N ) factors, which for this partial wave gives a bound which is stronger than the one obtained considering only one of the two SU (N ) factors, see Eq. (4.18) and Eq. (4.28). Overall in the case of the SM Yukawa sector the most stringent bound turns out then to arise from J = 0 partial wave. The limit of Eq. (5.3) implies that if there were additional quarks acquiring mass from EWSB, their mass should have been 500 GeV in order to preserve PU. Analogously, additional leptons should have a mass 700 GeV.

Multiple generations
We want now to highlight what is the role played by the presence of mutiple states with the same quantum numbers, as in the case of multiple generations of SM fermions. We then go back the the general case of Eq. (5.1) and, for simplicity, work with only two generations of fermions. Through a biunitary rotation acting on the fermion fields q L → U L q L and d R = U R d R it is possible to go to a basis where the Yukawa matrix becomes diagonal with real and non negative entries, namelyỹ whose largest eigenvalue is N 3 16π (ỹ 2 d,1 +ỹ 2 d,2 ) and the bound is thus on the geometric mean of the two Yukawa couplings This is not the case for the scattering in J = 1/2, where different generations do not communicate since the scatterings proceed through the exchange of an s− or u−channel fermion. In this case the eigenvalues give independent bounds on the two couplings separately, which readỹ 11) in total analogy with the single family case. In J = 1 the situation is again more involved due to the presence of the 00 two particle state which is common between all generations. The eigenvalues of the scattering matrix have a complicated analytical form, but numerically one can see that the strongest bound is always given by the J = 0 partial wave.

Majorana type theories
In this Section we study the Majorana type theories, described by Eq. (3.2). As discussed in Sec. 3 and indicated in Tab. 1, when only one SU (N ) factor is present, or flavor multiplicity is not added, only one model, other than the one where all fields are SU (N ) singlets, can be written. We firstly study these two theories in turn and then in Sec. 6.3 we present an explicit, phenomenologically relevant, example in which the involved fields are charged under multiple SU (N ) factors and the above caveat can thus be evaded. Given that the two models of Tab. 1 cannot be written for arbitrary SU (N ), in this Section we directly illustrate our findings without explicitly presenting the group factors F, as opposed to the thorough derivation of Sec. 4 for the Dirac type theories. We also do the same for the explicit example of Sec. 6.3. In this case the various amplitudes can be derived analogously to the examples of Sec. 4.
6.1 First model: In this model both the fermion and the scalar are SU (N ) singlets, where the latter is real if q = 0. Clearly, no group factor is present in this theory and the partial waves can be easily built directly from the Lorentz amplitudes of Tab. 3. In the complex scalar basis the Yukawa matrix of Eq. (2.6) is simply Y α = y, where α = 1 in the real scalar case and α = 1, 2 in the complex scalar one. Let's start again by discussing J = 0 with real S. In this case both the ± ± ±± and ± ± ∓∓ helicity channels contribute. The partial wave matrix is readily computed and, after integration in the (χ, χ) basis, reads which gives the bound Again, if the scalar is complex there is no scattering in the ± ± ∓∓ helicity channels. The partial wave matrix of Eq. (6.1) becomes diagonal and the bound relaxes to If S is a complex scalar there is now a contribution to the partial wave matrix from the 00 * + − scatterings. In the basis (SS * ,χχ) and after the angular integration the partial wave matrix is which gives the bound y 2 < 8π . (6.7) 6.2 Second model: In this model the fermion transforms in the fundamental of SU (2) while S in the Adjoint representation, and is then a real scalar if q = 0 and complex otherwise. In the first case we choose as basis 8) and the Yukawa matrix is then where α = 1, 2, 3 and i, j = 1, 2. When q = 0 then S is a complex field. In this case we can choose φ = (S A , S A * ) T , A = 1, 2, 3 , (6.10) and the Yukawa is now where now α = 1, . . . , 6. Proceeding in a similar manner as for the Dirac type models, we find that in J = 0 the bound is the same for both real and complex scalar, due to a cancellation between the s-, t-and u-channels in the ± ± ∓∓ amplitudes. Moreover, since for the ± ± ±± transition there is only the s-channel exchange of S, the only non-vanishing scattering has the fermions in the SU (2) triplet configuration, giving Moving to the J = 1/2 partial wave, again the strongest bound is the same for real and complex scalars, coming from the scattering in the 4 of SU (2): Finally, in J = 1, the best bound for real S is obtained in the adjoint channel, where the partial wave matrix reads with eigenvalues − y 2 256π (1 ± √ 65), thus giving the bound For complex S, instead, the singlet channel gives the strongest constraint. The partial wave matrix is and it has eigenvalues − y 2 256π (3 ± √ 57). The bound therefore is 6.3 The case of the S 1 leptoquark As already mentioned above, there are more possibilities for the Majorana type models once one allows for the fields to be charged under more than one SU (N ) group. Of particular phenomenological interest is the case of the leptoquark S 1 , that will be discussed in more detail also in Sec. 7. This field transforms under the SM gauge group as S 1 ∼ (3, 1, 1 3 ). One can thus write the following interaction term with the SM quark doublet where the colour indices are contracted with the totally antisymmetric tensor ε abc of SU (3), compensating the SU (2) contractionqεq c . The bounds on the coupling y in this case can be obtained along the same lines as the ones in the previous sections, and we therefore quote only the results for the three considered partial waves

Phenomenological applications
In this Section we apply our results to some illustrative models which present additional Yukawa interactions formulated to solve several anomalies reported in low energy measurements, such as the muon anomalous magnetic moment (g − 2) µ and the anomalies in the charged-and neutral-current decays of B−mesons, commonly dubbed as R D ( * ) and R K ( * ) anomalies respectively. The former is an anomaly in the B(B → D ( * ) τ ν)/B(B → D ( * ) ν) observable in b → cτ ν charged-current transitions, with = e, µ, while the latter is an anomaly in the B(B → K ( * ) µ + µ − )/B(B → K ( * ) e + e − ) observable in b → sµ + µ − neutralcurrent transitions. In order to explain the R D ( * ) anomaly, a ∼ 15% modification with respect to the theory prediction is required. However in the SM the partonic process b → cτ ν occurs at tree-level, hence when one tries to explain the experimental measured value through some additional NP contribution one might encounter several problems. Since the NP contribution to this observable scales, in case of a tree-level effect, as where g NP and m NP are the coupling and the mass of the relevant NP state, a large effect can be obtained either with a small NP mass or with a large NP coupling. However given that the suppression scale for the SM effective operator 1 Λ 2 (q 2 γ µ σ A q 3 )( 3 γ µ σ a 3 ) that can address this anomaly is Λ 3 TeV, in the former case one has to face stringent limits from direct searches from, e.g., the LHC, while in the latter case the coupling might be pushed at the edge of perturbativity. On the other hand the partonic process b → sµ + µ − entering the R K ( * ) anomaly occurs in the SM at one-loop level, with a V tb V ts CKM suppression. When considering NP models that try to explain this measurement also at one-loop level, again one can obtain couplings which might be in conflict with the requirement of perturbative unitarity. The purpose of this Section is to apply our results to phenomenologically relevant models and show that the requirement of PU can enforce significant bounds that might deserve further investigation. Since typically in the models that we will consider more than two couplings at the same time can enter the expression of the PU limit, our strategy will be to trade some of them for other measurements and/or constraints and then to depict the PU bound in the region of the two remaining independent couplings.

Scalars and fermions for R K ( * ) and (g − 2) µ anomalies
The first model that we study extends the SM by adding new scalars and fermions in order to generate contributions to b → sµ + µ − and (g − 2) µ , both at loop-level and it is based on [79,80]. We first consider the simplest extension which contains only left-handed (LH) couplings and then we evaluate the consequences of adding right-handed (RH) couplings.

Left-handed scenario
In the LH scenario, the NP states couple only to LH SM quarks and leptons. We can consider two models with the following schematic interactions 10 • Model a) with one additional scalar Φ and two additional fermions Ψ q and Ψ • Model b) with two additional scalars Φ q and Φ and one additional fermion Ψ where q and are the SM quark and lepton doublet respectively and where the NP fields quantum numbers under the SM gauge group are at this level unspecified. Here however we wish to assess how constraining the PU requirement could be and since, as shown in Sec. 4, bounds are generally stringent when the theory features a real scalar field, we wish to consider models that feature a real scalar. In model a), however, by making Φ a real scalar one obtains an exact cancellation of the various contributions to b → sµ + µ − [88], an option disfavored if one is willing to explain the R K ( * ) anomaly. This is not the case for model b), where one can choose Φ to be a real scalar. Altogether we consider the following quantum number assignments under the SM gauge group for the two models where by fixing X = 1 2 one has that Φ is a real scalar in model b). Regarding the flavor structure of the theory, since the goal is to generate a contribution to b → sµ + µ − , we only need couplings to the second and third quark families 11 , λ q 3 ≡ λ b , λ q 2 ≡ λ s and to the second generation of leptons λ 2 ≡ λ µ . The loop-level diagrams responsible for generating the NP contribution to b → sµ + µ − are shown in Fig. 3 By fixing for simplicity all the masses of the NP states at a common value m NP ∼ O(TeV), the most stringent bound for the couplings to quarks comes from B s −B s oscillations, where using the result in [117] we get 10 Note that one can also construct a model where Φ and Ψ couple to SM quarks while the conjugate fields Φ c and Ψ c couple to leptons. This however leads to very similar phenomenological results. 11 We work in down-quark aligned basis. Figure 3. Loop-level diagrams responsible for generating the NP contribution to b → sµ + µ − for the models of Eq. (7.2) and Eq. (7.3).
This relation can be inserted in the expression for the ∆C µ 9 = −∆C µ 10 coefficients 12 for reproducing the neutral-current anomaly R K ( * ) from where one has for both model a) and model b). By plugging Eq. (7.5) into Eq. (7.6) one can set a lower bound on the |λ µ | coupling [80] |λ where we use the updated 1-dimensional fit ∆C µ 9 = −∆C µ 10 = −0.41 ± 0.07 in [118]. Hence, by saturating the bound in Eq. (7.5), i.e., by imposing |λ s λ b | ∼ 0.15 m NP 1 TeV , we can compute the bound set by PU to see if one can explain at the same time the observed value for ∆C µ 9 relevant for the R K ( * ) anomaly. In order to do so we fix m NP = 1 TeV, which for SU (3) c charged NP states is at the edge of exclusion from direct searches at the LHC, and plot the allowed regions from PU in the (λ µ , λ b ) parameter space, accounting for the ∆C µ 9 value from [118]. We illustrate this for model a) in Fig. 4 where, as explained before, the scalar is a complex field. There in green (yellow) we illustrate the regions compatible with the measured value of ∆C µ 9 at 1σ and 2σ while in gray we show the one compatible with PU. We see that, in this case, there is an overlap between the two regions and the R K ( * ) anomaly can be explained with couplings whose magnitude is compatible with perturbative unitarity.  Since in order to solve the neutral current anomaly we need a coupling to the muons, it is natural to ask if one can reproduce the (g − 2) µ anomaly and how large the relevant coupling has to be to achieve the correct NP contribution. For the observable a µ ≡ (g − 2) µ /2 we consider the recent value of the Fermilab Muon g − 2 experiment [15] for which one has a ∼ 4.2σ discrepancy with respect to the SM prediction [17] ∆a µ = a exp µ − a SM µ = (251 ± 59) × 10 −11 . (7.9) We want to see if this anomaly can be explained with the λ µ coupling in a perturbative regime. To illustrate this we consider the case of model a), since it's the one for which one has a less stringent PU bound. Using the results in [80], it turns out that to explain the muon anomaly at the 1σ level by saturating λ µ = 4.0 as per Eq. (7.8) one needs to have For a common NP mass of 1 TeV one needs a quite exotic and large value for the hypercharge X = 10.6. In this case one has to assess the validity of the perturbative regime studying scattering processes that involve gauge bosons. Even more extreme hypercharge values are needed for smaller values of λ µ 13 . The need for a large muon coupling in order to explain the (g − 2) µ anomaly arises because the process needs a chirality flip, which can be obtained in this LH model only via the muon mass term with a contribution proportional to m µ , which forces the couplings to be too large to account for the anomaly. This fact could be solved by introducing NP that couples to RH SM muons together with a mixing term among the NP states, since in this case one can generate a chirality flip proportional to m NP m µ , that allows for a smaller NP coupling. We present some models which include RH couplings in the next Section.

The inclusion of right-handed couplings
The possibility of adding RH couplings to scalar-fermion models has been largely discussed before in the literature, also in the context of DM physics, since for some choices of the field representations one can have suitable DM candidates [88,92,[119][120][121]. Introducing a coupling to RH leptons requires at least one new scalar or fermion field. A mixing term among the NP fermion fields can be generated through the interaction with the Higgs boson, while dangerous mixing terms between the Higgs, a SM and a NP field can be forbidden by introducing an extra symmetry like a Z 2 symmetry or a U (1) charge. As explained before, the motivation for introducing RH couplings in these kind of models is to be able 13 In [80] other representations for the fields are discussed, where for a value of hypercharge |X| = 1 the muon anomaly can be accounted for with |λµ| ≥ 3.7, although in that case the PU limit also tightens to |λµ| ≤ 2.5 to account for the (g − 2) µ anomaly, while keeping the NP couplings in a perturbative regime. This can be achieved provided that we have a chirality flip contribution bigger to the one proportional to the muon mass. In a recent work [92], the Authors investigate two models containing a good DM candidate while explaining at the same time the b → sµ + µ − and the (g − 2) µ anomalies. In particular one of the two scenarios is the extension of model b) of Sec. 7.1.1 with a real scalar Φ , whose Lagrangian reads (7.11) where we labeled explicitly the chirality indices L, R in the new fermions Ψ, Ψ . The field quantum numbers that we consider in this case are Again, we restrict our analysis to the case where the flavor structure enforces only couplings to b and s quarks and to muons. Hence we are left with 5 parameters that will allow us to explain the muon anomalous magnetic moment: λ b , λ s , λ µ , λ e µ and λ H . By fixing, e.g., λ H = 0.1 it turns out that one can explain the (g−2) µ anomaly with perturbative couplings. In order to assess whether the neutral-current anomaly can be explained in this scenario while remaining in the perturbative regime we proceed similarly to the case of the LH scenario and start by saturating the B s −B s bound in Eq. (7.5), which fixes the quark coupling combination |λ b λ s |. For what concerns the b → sµ + µ − observable, we now have no longer the ∆C µ 9 = −∆C µ 10 pattern, so that we take the 2D fit result from [118] ∆C µ 9 = −0.68 ± 0.16 ∆C µ 10 = 0.24 ± 0.13 , (7.13) and the ∆C µ 9 and ∆C µ 10 expressions from [92]. We then focus on two benchmark points presented in [92], that can account for both the (g − 2) µ and the R K ( * ) anomalies and the DM relic density, while being compatible with the bounds from direct searches, namely which uniquely fix the values of (λ µ , λ e µ ). We then show in gray in Fig. 6 the region compatible with the requirement of PU in the (λ µ , λ e µ ) plane, as well as the allowed region for reproducing (g − 2) µ at 1σ and b → sµµ at 2σ, which are depicted in brown and yellow respectively. In the figures we also show the (λ µ , λ e µ ) for both benchmark points. Since this model features a real scalar field and the λ e µ needs to be small to satisfy flavor observables, (g − 2) µ and R K ( * ) , the PU bound is again dominated by the results of the

Scalar leptoquarks
Scalar LQs are a natural candidate to explain the charged-and neutral-current R D ( * ) and R K ( * ) anomalies since they couple quarks to leptons, and are thus an ideal scenario to be tested with the tool of perturbative unitarity. Among all the scalar LQs the SU (2) L triplet S 3 and singlet S 1 are the most robust candidates to explain the anomalies [43, 46, 53, 64, 67, 69-73, 76, 77]. Under the SM gauge group they transform respectively as S 3 ∼ (3, 3, 1 3 ) and S 1 ∼ (3, 1, 1 3 ). When both LQs are combined so as to explain both the B−meson anomalies and the (g − 2) µ , and their mass is set to O(1) TeV, the SM discrepancies can be explained without suffering from PU constraints. However, for higher masses, the couplings are required to be tuned to higher values and perturbativity might be lost. In a recent work [76] the Authors have considered the following SM extension where φ + is an SU (3) c and SU (2) L singlet scalar with Y = 1. This model aims at explaining the B−anomalies, the anomalous magnetic moment of the muon and the so called Cabibbo Figure 7. Left: In gray we show the allowed region from PU and in purple the region compatible with R D ( * ) at 2σ. The best fit point in [76] is shown in blue. Right: In gray we show the allowed region from PU while in green and yellow the regions compatible with R K ( * ) at 1σ and 2σ respectively. Between the dashed red lines constraints from τ → µγ are satisfied. The best fit point in [76] is shown in blue.
Angle Anomaly [122][123][124] with the following flavor structure: (7.17) Setting m NP ≡ m S 1 = m φ + = 5.5 TeV, the best fit point of this model appears in Eq. (12) of [76]. Taking the best fit values, each coupling turns out to be in the perturbative regime when considering one of them at the time. However, when considering the contribution from all the couplings simultaneously, perturbative unitarity is lost. This is mainly due to the fact that λ u cτ has to be very large in order to account for the R D ( * ) anomaly. Explicitly one has [59,125]  where ∆R D ( * ) =R exp D ( * ) − R SM D ( * ) . We then try to see whether it is possible to explain all the anomalies considered in [76] while remaining in a perturbative regime, without changing the values of m NP = 5.5 TeV chosen by the Authors. In the left panel of Fig. 7 we show in the (λ q bτ , λ u cτ ) plane the region compatible with PU, depicted in gray, and the 2σ region where R D ( * ) can be reproduced, depicted in purple, where the other couplings are set to their best fit value so that the other relevant anomalies, (g − 2) µ and the Cabibbo Angle, can be reproduced. We see that there is a small region where R D ( * ) can be satisfied while being compatible with PU, provided that we lower the value of |λ u cτ | from the best fit value to |λ u cτ | = 2.5 and keep the best fit value for λ q bτ . Thus, by fixing λ u cτ = −2.5 and keeping the other couplings at the best fit as indicated in [76], we show in the right panel of of Fig. 7 in the (λ u cµ , λ µτ ) plane the region where the R K ( * ) anomaly can be reproduced at 1σ (green) and 1σ (yellow). There we see that there is compatibility between this requirement and the one of PU, although with a slightly different benchmark point than the one of [76]. It is important to mention that the coupling λ u cµ is introduced in the model in order to cancel undesired effects in τ → µγ due to the large value of λ u cτ . For this region in the right panel of Fig. 7 we include the region allowed by τ → µγ [126], to show the compatibility with this latter measurement.

Yukawa sector in vector leptoquark models
Other than scalar LQ, a compelling possibility to simultaneously solve the R D ( * ) and R K ( * ) anomalies is through a vector LQ. The most remarkable candidate is the U µ 1 vector with SM quantum numbers U µ 1 ∼ (3, 1, 2/3) which has triggered a large theoretical activity aiming at providing an UV completion [34, 35, 37, 46, 48-52, 56, 59, 60, 68, 78]. Generally, in order to address the flavor anomalies, the models including U µ 1 also require the presence of new vector-like fermions and new scalars that couple to the SM via Yukawa couplings which can be constrained by PU considerations. Here we focus as an example on the model presented in [60], usually dubbed in the literature as 4321 model , since it possesses a gauge symmetry G = SU (4) × SU (3) × SU (2) × U (1). The Yukawa part of the theory can be divided in a SM-like part and a part which includes the NP fields L = L SM−like + L mix . Explicitly where we refer to [60] for the field definitions and their quantum numbers under G. Here we focus on the last term of L mix which contains the mixing between the new vector-like fermions Ψ and the scalar Ω 15 , whose quantum numbers under G are By computing the PU unitarity bound one obtains that the strongest limit is obtained from the J = 1 channel and reads λ 15 2.1 . Regarding the viability of perturbative couplings of the 4321 model , while in the original work [60] the Authors set λ 15 2.5 in order to introduce a mass-spliting between new heavy vector-like quarks and leptons, which would then be in contrast with the perturbative unitarity limit that we have derived, with the new experimental world averages for the R D ( * ) and R K ( * ) anomalies, one can easily lower the Yukawa coupling to, e.g., λ 15 ∼ 2, while remaining compatible with ∆F = 2 observables. Thus the model is still viable, although the parameters are stretched to the edge of perturbativity according to our criteria.

Right-handed neutrinos for B anomalies
There are also models that can account for the anomalies with the addition of a RH neutrino, thus connecting the flavor tensions with the one of the neutrino mass generation. In [59] the authors have proposed a model that can address the R D ( * ) anomaly by adding a new decay channel B → D ( * ) τ N R into a right-handed sterile neutrino N R while simultaneously solving the R K ( * ) anomaly at one-loop level, through the exchange of a scalar leptoquark S 1 . The Lagrangian of the theory is In order to reproduce the neutral current anomaly one has to tune (7.25) in order to avoid violation of lepton flavor universality in b → c ν processes, see again [59], where V is the CKM matrix. With this tuning one has that the neutral-current anomaly is reproduced for |λ q bµ | 2 0.87 + 3.15 where we have normalized the expression to the latest best fit for the ∆C µ 9 coefficient [118]. Barring the mass of the RH neutrino, there are four couplings and one mass in this model.  Figure 8. In black we show the region compatible with PU of the Yukawa couplings λ d bN and λ q bµ for a LQ mass of 1 TeV (left) and 2 TeV (right). The brazilian band represents the region compatible with the R K ( * ) anomaly at 1σ and 2σ, while in the purple region the measured value R D is reproduced, having fixed λ u cτ so as to reproduce R D ( * ) . A LQ mass of 1 TeV is excluded by experimental bounds on B s mixing.
One parameter is eliminated by the tuning of Eq. (7.25), while we can eliminate, e.g., the value of λ u cτ by asking to reproduce the R D * anomaly, which is the one with the smaller experimental error. This leaves two independent couplings, λ d bN and λ q bµ , on which we can check the constraints imposed by PU. We show the results in Fig. 8 for two representative values of the LQ mass. In those figures the region compatible with the PU of the Yukawa couplings is shown in gray, while the brazilian band plot illustrates the region of parameter space that can explain the R K ( * ) anomaly at 1σ and 2σ. Finally in purple we show the 2σ region compatibility for R D , having fixed λ u cτ so as to reproduce R D ( * ) . Altogether we see that for a LQ mass of 1 TeV (left panel) we can simultaneously explain both anomalies while remaining in the perturbative regime. However for this value of the LQ mass, the solutions to the R K ( * ) anomaly is excluded by the experimental bounds on B s mixing, see again [59]. We can restore the compatibility with this measurement by raising the LQ mass up to 2 TeV (right panel), where however now the λ d bN and λ q bµ couplings are pushed at the edge of the perturbativity.

Conclusions
Yukawa interactions are ubiquitous in NP theories that try to address the shortcomings of the SM and are largely employed in models that try to solve experimental anomalies reported in the recent years in low energy data, as for the case of the muon (g − 2) µ and semileptonic decays of B−mesons. In this paper we have studied the constraints imposed by PU on generic Yukawa interactions where the fields involved have arbitrary quantum numbers under an i SU (N i ) ⊗ U (1) group.
By considering all 2 → 2 tree-level scatterings in the high-energy limit we have constructed the general form of the partial-wave matrices a J f i and derived compact expressions for the upper limit on the value of the Yukawa interaction up to which perturbation theory can be trusted. This has been achieved by computing all the necessary ingredients for building the partial-wave matrix, namely the Lorentz parts of the scattering amplitudes and the group structure factors entering the amplitude themselves. We have started by considering a set of phenomenologically relevant toy models with Dirac type and Majorana type interactions, where the various fields are only charged under a single SU (N ) factor, working for concreteness in the case where all the fields transform in the trivial, fundamental or adjoint representation of SU (N ) and allowing them to have arbitrary U (1) charges. We have shown how the SU (N ) group structure of the interaction can lead to an enhancement of the scattering amplitudes and thus to a tightening of the partial wave unitarity bound, while on the other hand the presence of the U (1) symmetry enforces a selection rule that makes some amplitudes vanish. Interestingly, we obtained that the stronger bound might arise from a partial wave different from J = 0.
The results obtained for these toy models can then be used as building blocks for more complicated theories, where the the various fields are charged under multiple SU (N i ) factors. To highlight the strategy we have provided a guided working example, by focusing on the case of the SM quark Yukawa sector. For this case we have also stressed the role that a non trivial flavor structure has in determining the PU bound. We have then applied our results to various more complicated NP models which solve the aforementioned anomalies in (g − 2) µ and/or semileptonic B−meson decays by postulating the existence of new Yukawa interactions. We have highlighted that, while the proposed theories can generally still provide an explanation to these measurements, their models parameters are stretched close to the limit where perturbation theory cannot be trusted and care must be taken in deriving any conclusion.
Finally, the results presented in this paper and illustrated in Figs 1, 2, 9 and 10, are of practical use, and their applicability lies beyond the simple examples presented in the text. While we have restricted only to a limited number of irreducible SU (N ) representations under which the various field can transform, the expressions that we have derived furnish the necessary ingredients to study the limits imposed by the requirement of PU in a large set of phenomenologically relevant NP theories that present additional Yukawa interactions.

A.1 Wigner d−functions
The small Wigner d−functions are defined in the angular momentum basis as whereĴ y is the generator of the rotations around the y−axis. The explicit expression of these functions used throughout our analysis are

A.2 Helicity spinor formalism
The fields entering Eq. (2.6) can explicitly be expanded in terms of creation and annihilation operators as where we choose the spinor basis to be where σ µ = (1 2 , σ i ),σ µ = (1 2 , −σ i ) and σ i are the Pauli matrices. We can choose ξ to be an eigenstate of σ 3 , i.e. ξ + = (1, 0) and ξ − = (0, 1) corresponding to spin up and down along the z-direction and we fix η + = (0, 1) and η − = (−1, 0) with the same convention. By building the helicity operatorλ where s, r = ±1 indicate helicity ± 1 2 for both particle and antiparticle and where for the antiparticle the helicity is defined with the opposite sign according to standard definitions of helicity spinors, see e.g [97]. Then the field ψ L of Eq. (A.4) annihilates negative helicity states ψ − and creates positive helicity states ψ + while the conjugate field annihilates positive helicity states ψ + and creates negative helicity states ψ − . To compute the relevant amplitudes we also need rotated spinors that can be built as is the rotation matrix in the x − z plane by an angle θ with respect to the y−axes. An analogous expression holds for v s (p ). where the s, t and u supscripts indicate the Mandelstam channel through which the relative amplitude proceeds.

B Other Dirac type theories
Here we present the results for model 3,4 and 5 in the Dirac type class.
B.1 Third model: In this model S transforms under the adjoint SU (N ) representation and is thus a real field if q = q , complex otherwise. We choose as basis where the index a and A run from 1 to N and N 2 − 1 respectively. With this choice the Yukawa matrix in the real scalar case reads with i, j = 1, . . . , 2N and α = 1, . . . , N 2 − 1 while in the complex scalar case it is instead where now α = 1, . . . , 2(N 2 − 1). In the case of J = 0 the relevant two particle states decompose again as in Eq. (4.3). The relevant group factors for the non zero amplitudes in the real scalar case are The partial wave matrix for J = 0 in the (χη, χη) basis 14 after integration on the angular variable is where × denotes the Kronecker product. For N > 1 the largest eigenvalue of this matrix comes from the singlet channel and is equal to y 2 32π N 2 −1 N , which thus gives the bound If S is a complex scalar again the amplitudes in the ± ± ∓∓ channels are zero. The only non vanishing scatterings when the matrices in Eq. (B.5) are diagonal are the ones in the adjoint channel. The largest eigenvalue is y 2 32π and the perturbative bound becomes Moving now to the scattering in the J = 1/2 partial wave the two particle states now decompose as where r 1 and r 2 are the two irreducible representations arising from the tensor decomposition × Adj = + r 1 + r 2 . In tensor component this reads having indicated with A i and B ·j k two tensors transforming in the antifundamental and adjoint SU (N ) representation. The first line of Eq. (B.8) indicates the fundamental representation, while the second and third are symmetric and antisymmetric tensors in i, j with null traces with respect to k. As an example, in the case of SU (3) this reads 3 ⊗ 8 = 3 ⊕ 6⊕ 15, and in SU (4) it is 4 ⊗ 15 = 4 ⊕ 20 ⊕ 36. The two-particle state in the fundamental representation is easily built as For the other two irreducible representations r 1,2 one needs to build by hand the basis for the vector space. Let's start with the representation with the higher dimension r 2 . Here one can split the vector space of the last line of Eq. (B.8) in three categories. Tensors where i = j = k, which are trivially traceless, tensors with i = j, i = k, which again are trivially traceless, and tensors which are traceless but where the null trace arise because of the sum of non zero elements 15 . One can count the dimensionality of these three categories to be , N (N − 1) and N (N − 1) respectively, whose sum is (N +2)N (N −1)

2
, matching the dimensionality of r 2 . For r 1 , instead, one has that only the tensor with i = j are non vanishing due to the antisymmetry in those indices. One can build then two categories for the tensor basis with dimensions N (N −1)(N −2) 2 and N (N − 2), whose sum is N 2 (N 2 − N − 2) which matches the dimension of r 1 . Note that this representation vanishes for the case of SU (2). In order to compute the group factor entering the scattering amplitude it's enough to explicitly build only one of this states, since all of them will give the same result. For example we construct the state with unit norm belonging to the first category for r 2 as This works for N ≥ 3, since for SU (2) it is not possible to have i = j = k. In order to be able to compute the scattering in the representation r2 also for N = 2 (r2 = 4), one can construct e.g. the states with i = j = k as |Sψ IK where I, J, K label the irreducible representation and range from 1 to d r 2 , i, j, k range from 1 to N and A from 1 to N 2 − 1. Analogously one can build the state in r 1 of the same category as The group factors entering the amplitudes for J = 1/2 turn out to be +0 + 0 In the case of a real scalar S and considering the scattering in the fundamental channel and for the +0 + 0 helicity amplitude, which is then ηS → ηS, one has explicitly For what concerns the scattering in the r 1,2 channels, they all have ± y 2 32π eigenvalues. The bound is thus If the scalar is complex one has that in the +0 + 0 helicity channel the ηS → ηS and χS →χS scatterings proceed through s−channel, while the ηS * → ηS * andχS * →χS * through u−channel. Considering the s−channel diagrams in the fundamental channel the eigenvalue is y 2 64π N 2 −1 N while the u−channel diagrams in the r 1,2 channel have eigenvalue ± y 2 32π . The bound is thus Finally in the J = 1 channel, while the +− two-particle states decompose again as in Eq. (4.3), the same is not true for 00, since the scalar field now belongs to the adjoint representation. Here one has the decomposition Adj ⊗ Adj= 1 ⊕ Adj ⊕ Adj + . . . . Given the irreducible representations that can be built out from two fermions, only the singlet and adjoints channels are relevant. For real scalar fields the two particle states can be built as where f ABC and d ABC are the antisymmetric and symmetric SU (N ) structure constant respectively 16 . However, both the singlet and the symmetric adjoint state do not contribute in J = 1 [115]. The relevant group factors read The channels with the highest eigenvalues are the singlet and adjoint ones, for which the partial wave matrix explicitly reads, in the (χχ, ηη, SS) basis, where for convenience we have written explicitly the expressions of the T amplitudes and where 2 csc θ = 1 tan θ 2 + tan θ 2 is the angular part arising from the sum of the t− and u− channel contributions which have a different sign because of the antisymmetry of f ABC entering the amplitude computation. We further note that, as expected, the matrix in the adjoint channel is complex but hermitian, yielding thus real eigenvalues. The strongest bound comes from the adjoint channel for N < 5 and from the singlet channel for N > 5 and reads When the scalar is a complex field one has a non zero contribution also in the singlet channel of the 00 * + − scattering, since one can build a non symmetric state, and the scatterings proceed now via the u−channel. Also in the adjoint channel the contribution from the symmetric state built with the d ABC structure constant, see Eq. (B.17), no longer vanishes. The correct normalization for the state is now 16 Note that d ABC identically vanishes in SU (2). and the group factors become where we have indicated with Adj−f and Adj−d the two contributions from the 00 twoparticle states built with the antisymmetric and symmetric SU (N ) structure constants respectively. For simplicity we only write the amplitude in the singlet channel, which is the one yielding the stronger limits: From the largest eigenvalue of this matrix one obtains the bound Altogether the limits arising from the various partial waves are reported in Fig. 9.
B.2 Fourth model: χ ∼ q , η ∼ Adj q , S ∼ q−q In this model S is always a complex scalar field and by fixing the basis ψ L = (χ a , η c,A ) T , φ = (S a , S * a ) T , a = 1, . . . , N, , A = 1, . . . , N 2 − 1 , (B.25) where a and A are indices running from 1 to N and N 2 − 1 respectively, the Yukawa matrix reads As in the previous cases, in the J = 0 partial wave there is no scattering in the + + −− sector because the scalar is complex, so we focus on the + + ++ channel where only theχη →χη process is non vanishing. We then use the tensor decomposition for × Adj = + r 1 + r 2 of Eq. (B.8) and built the two particle states analogously to Eq. (B.9), Eq. (B.11) and Eq. (B.12). One obtains non vanishing amplitudes only in the fundamental channels + + + + F s, χηχη = N 2 −1 2N , (B.27) In J = 1/2 we can again focus only on the +0 + 0 helicity amplitude. Here one decomposes the two particle states as + 0 ηS ∼ + r 1 + r 2 χS ∼ 1 + Adj ηS * ∼ + r 1 + r 2 χS * ∼ S + AS , and the group factors for the non vanishing scatterings are Also in this case the channel yielding the stronger limits depends on the value of N . One obtains where the first comes from the scattering in the antisymmetric channel while the second from the one in the fundamental. Finally, for J = 1 we can decompose the two-particle states as For brevity we report only the partial wave in the singlet channel, which is the one giving the most stringent bound, which reads, in the basis (χχ, ηη, SS * ), a J=1 = y 2 32π (B.32) The eigenvalues of this matrix have a complicated form and we report the numerical results in Fig. 10, together with the limits from the other partial waves. In the J = 0 partial wave only the ±±±± scatterings proceeding through the antisymmetric channel are non zero, due to the antisymmetry of ε. The group factor is simply 2 and the bound turns out to be y 2 < 4π . (B.35) 17 We fix ε123 = 1.
For J = 1/2 the two particle states decompose as Finally for J = 1 all the two-fermion states decompose as ∼ 1 + Adj and the same is true for the SS * scalar state, which is the only one which leads to non zero amplitudes. The group factors for + − +− are all ±(N − 1) for scatterings in the singlet channel and ±1 for scatterings in the adjoint one. On the other side in the 00 + − helicity channel for the group factors one obtains ±2 in the singlet channel and ±1 for the adjoint one. The most stringent bound is then obtained from the scattering among singlets where one has, in the (ηη, ηχ,χη,χχ, SS * ) basis and after the angular integration, which gives the bound which is the most stringent among the various partial waves.