Spontaneous CP Violation and Horizontal Symmetry in the MSSM: Toward Lepton Flavor Naturalness

We study the contributions of supersymmetric models with a $U(1)$ horizontal symmetry and only spontaneous CP breaking to various lepton flavor observables, such as $\mu \to e\gamma$ and the electron electric dipole moment. We show that both a horizontal symmetry and a lack of explicit CP violation can alleviate the existing bounds from such observables. The undetermined $\mathcal{O}(1)$ coefficients in such mass matrix models muddle the interpretation of the bounds from various flavor observables. To overcome this, we define a new fine-tuning measure for different observables in such setups. This allows us to study how naturally the observed IR flavor observables can emerge from a given mass matrix model. We use our flavor-naturalness measure in study of our supersymmetric models and quantify the degree of fine tuning required by the bounds from various lepton flavor observables at each mass scale of sleptons, neutralinos, and charginos.


Introduction
The origin of fundamental particle masses is a central issue in modern particle physics. The Standard Model (SM) partly addresses this issue by the higgs mechanism. However, it leaves two major questions unanswered. The first question is why the three generations of quarks and leptons couple to the higgs field so differently. The couplings range from O(10 −5 ) to O(1), and it seems unlikely that they are determined completely randomly. The second question is on the stability of the scale of the electroweak symmetry breaking (EWSB) under quantum corrections, commonly known as the naturalness problem or the hierarchy problem [1][2][3][4]. A significant fine tuning is required for the higgs mass-squared parameter to realize the electroweak scale much smaller than UV mass scales such as the Planck scale. These two questions have driven numerous theoretical and experimental studies of physics beyond the SM (see Refs. [5,6] for reviews).
The stability of the electroweak scale is elegantly ensured by supersymmetry (SUSY). In SUSY models, a large quantum correction from a SM particle to the higgs mass-squared parameter is canceled with that of its corresponding superpartner. If SUSY is realized in nature, it must be broken; however, SUSY breaking can radiatively drive the higgs masssquared parameter negative so that the electroweak symmetry is dynamically broken. While the discovery of the higgs boson with the mass of around 125 GeV at the Large Hadron Collider (LHC) indicates that the superpartner of the top quark is significantly heavier than the weak scale (see Ref. [7] for a recent review, and Refs. [8][9][10] for early references), to avoid a severe tuning, the mass scale of supersymmetric particles should be within just a few orders of magnitude above the electroweak scale (for a recent study of tuning of the minimal supersymmetric SM (MSSM) or other SUSY extensions of the SM, see Ref. [11]).
Nonetheless, a random breaking of SUSY generates significant flavor-and CP-violating processes, in contradiction with precise experimental data. One possible solution to this problem is to engineer a mechanism to realize flavor independent masses of superparticles. Gauge mediation is the most popular example (see, e.g., Refs. [12,13] for reviews), but it predicts a light gravitino which generally leads to cosmological problems such as overclosure of the universe [14][15][16][17]. Gravity mediation with a cosmologically harmless gravitino of O(10) TeV appears to be more natural, but it is not equipped with a mechanism to suppress flavor and CP violating processes.
An intriguing approach to both the SM Yukawa hierarchy and the SUSY flavor problem is provided by a U (1) horizontal symmetry and its spontaneous breaking [18]. Depending on U (1) charge assignments of quark and lepton supermultiplets, their Yukawa couplings are suppressed by some power of a small parameter ≡ S /Λ where Λ is some UV mass scale and S is a vacuum expectation value (vev) of a flavon field S charged under the U (1). Appropriate charge assignments can lead to the observed SM fermion mass and mixing pattern. Since the same symmetry is applied to SUSY breaking masses of superparticles, the mass matrices of quarks and leptons and their superpartners can be approximately diagonal in the same basis. This alignment can suppress dangerous flavor violating processes induced by supersymmetric particles [19][20][21][22][23]. The U (1) symmetry can arise as an accidental symmetry enforced by discrete gauge symmetries, or as a low-energy remnant of an anomalous gauged U (1) obtaining a mass from the 4d Green-Schwarz mechanism [24]. These options allow very good approximate global symmetries, while being compatible with the expectation that exact global symmetries do not exist in gravitational theories [25][26][27].
The quark sector with three generations contains a CP violating complex phase called the Cabibbo-Kobayashi-Maskawa (CKM) phase. Recently, a nonzero phase has been hinted at in the lepton sector as well [28][29][30]. The existence of CP violation (CPV) in nature and the SUSY CP problem are reconciled by the idea of that CP is a gauge symme-try [31,32], which is broken spontaneously [33]. Multiple flavon fields can realize this idea: their vevs are complex and break CP spontaneously. In Ref. [33] this mechanism gives rise to the correct CKM phase, while suppressing the neutron electric dipole moment (EDM) through CP violating squark mass matrices. Taking account of a hinted phase in the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix and the stringent upper bound on the electron EDM (eEDM) [34], 1 it is natural to apply the same idea to the lepton sector.
The µ → eγ process is a well-known observable which gives a strong constraint on new physics with lepton flavor violation (LFV). The current most stringent limit on this process comes from the MEG experiment [37]. Future experiments such as the upgraded MEG II experiment [38] will increase the sensitivity by an order of magnitude. Moreover, the µ → e conversion process has been currently searched for and its future limit is expected to be dramatically improved. Thus, it is timely to assess up-to-date constraints from CPV and LFV observables on supersymmetric models with a U (1) horizontal symmetry and discuss future prospects.
In this paper, we consider a U (1) horizontal symmetry and its spontaneous breaking to realize the correct SM fermion masses and mixings, including CP violating phases, in the SUSY framework. We focus on the lepton sector. We do not study the squarks and gluinos here, since they will not contribute to the observables of interest at leading order. A detailed study of the quark sector observables in our setup is left for future works.
To accommodate a nonzero phase in the PMNS matrix without inducing sizable EDMs, the idea of spontaneous CP violation [33] is extended to the lepton sector. We explore current constraints from CPV and LFV observables on models with two flavons and discuss future prospects of the constraints. The allowed mass scale of supersymmetric particles is lowered in the models with two flavons compared to the case with anarchic soft masses containing uncontrolled complex phases. 2 The assessment of constraints on models with horizontal symmetries is generally obscured by uncertainties originating from undetermined O(1) numbers in mass matrices. Horizontal symmetries only tell us the scaling of mass matrix entries with the parameter , while undetermined O(1) numbers multiplied in each entry can have profound effects on SUSY contributions to CP and flavor violating processes as well as the obtained SM fermion masses and mixings. We improve this situation by proposing a new fine-tuning measure to assess the CP and flavor constraints that captures the effects of undetermined O(1) numbers. The new measure makes it possible to estimate constraints on models with horizontal symmetries in a more numerically precise way. This technique is generic and not limited to the present 1 See also Refs. [35,36] for interpretations of these bounds on SUSY models' parameter space. 2 By "anarchic," we mean completely uncontrolled by any flavor symmetry, as opposed to the well-known neutrino anarchy scenario [39] in which the left-handed leptons have a random mixing matrix but the righthanded leptons' mass matrix is structured. SUSY models. Similar tuning measures can be defined for other SUSY or non-SUSY models that explain the SM flavor structure.
The rest of the paper is organized as follows. In Sec. 2, we summarize the SM flavor structure, as well as the current and projected future bounds from measurements of eEDM and LFV processes. The severity of the SUSY CP and flavor problems is illustrated with the case of anarchic sfermion masses. In Sec. 3 we introduce a supersymmetric model with a horizontal symmetry and spontaneous CPV and discuss the effect on suppressing disastrous eEDM and LFV processes. We also introduce a new notion of naturalness for these observables and use it to quantify their bounds on a handful of sample charge assignments. Section 4 is devoted to conclusions and discussions. We also include further details in three appendices. We review the contribution of generic UV models to the relevant dipole operator in the IR in App. A and list the diagrams contributing to this operator in a generic MSSM model in App. B. We also elaborate on how we generate random numbers in the mass matrices in our numerical scans in App. C

Flavor & SUSY
In this section, we review the existing literature on various ingredients in our study. We start by reviewing the lepton flavor structure in the SM and the current and projected future bounds on the leptonic CPV and LFV observables that we study in this paper. The SUSY flavor problem and contributions to CPV and LFV observables in the MSSM are then summarized. We will see that the current bounds on these observables severely constrain the SUSY parameter space without further structure or severe tuning; we will systematically study new structures that may open up the parameter space in the next section.

The SM flavor structure
There are clear hierarchies in the flavor parameters of the SM. Within the strong sector, the quark masses span five orders of magnitude. Moreover, the diagonal elements of the CKM matrix are of order one, while off-diagonal elements are suppressed. Those hierarchies suggest the existence of a structure with a fundamental small parameter λ ∼ 0.2, originally introduced by Wolfenstein [40] for the CKM matrix elements, and later applied to the quark masses hierarchies [21].
Within the lepton sector, the charged lepton masses are hierarchical. The mass ratios can be parameterized by the Wolfenstein parameter λ as  The absolute scale of neutrino masses is still unknown, but a hierarchy is observed in the ratio of mass differences for the normal (inverted) hierarchy [41]. In contrast to the CKM matrix, the PMNS matrix of the lepton sector does not show any clear structure or hierarchy [41]. The PMNS matrix has one (three) CPV phase in the case of Dirac (Majorana) neutrinos. Only the Dirac phase δ CP can be measured by oscillation experiments. There is some tension between measurements of this phase. The NOνA [29] experiment favors δ CP < π, still allowing a CP-conserving phase, while T2K [30] favors δ CP > π and disfavors CP-violation at the level of 2σ.

Lepton flavor observables
There are many IR observables that can probe deviations from the SM flavor structure outlined above. In Table 1 we show a handful of such observables, including both current constraints and those that are expected in the next decade. The current strongest constraint on the eEDM comes from the ACME experiment, using the thorium monoxide (ThO) molecule [34]. Proposed future experiments are expected to improve the experimental sensitivity by a few orders of magnitude in the near future; see, for instance, Refs. [42][43][44][45][46][47]. In terms of a simple parametrization of the scale probed by an operator, the eEDM is the most sensitive lepton-sector observable, as indicated by the scale Λ quoted in Table 1. However, we will see that model-dependent details can significantly modify this conclusion.
A number of stringent constraints arise from charged lepton flavor violation observables probing µ → e transitions. A useful, recent summary prepared as an input to the European Strategy for particle physics may be found in Ref. [48]; see, especially, its Fig. 1 for the sensitivity of several planned experiments. The current most stringent bound on the LFV observable BR (µ → eγ) comes from the MEG experiment [37]. Approximately an order of magnitude improvement is projected for this observable after three years of data-taking at the upgraded MEG II experiment [38]. The conversion of a muon to an electron in the Coulomb field of a nucleus has been constrained by the SINDRUM II experiment [49]. Orders of magnitude increase in sensitivity are expected in the near future as well from experiments such as Mu2e [50] or COMET [51,52]. The SINDRUM experiment also provides the current most stringent limit on the BR (µ → 3e) process [53]. A new search for this process at the Mu3e experiment will improve the limit by a few orders of magnitude [54,55]. All three charged LFV processes (µ → eγ, µ → 3e, and µN → eN ) are sensitive to similar dipole operators (see Eqs. (6) -(8) below). Taking into account only these operators, the current most stringent constraints arise from µ → eγ, but will be surpassed by µN → eN and µ → 3e in the future. Because µN → eN can receive additional important contributions from box diagrams in portions of the MSSM parameter space, we will focus on the dipole operator contributions to µ → 3e as our benchmark for sensitive future charged LFV constraints.
There are other similar observables, e.g., BR (τ → µγ), that will be studied in the near future. In this paper we neglect these observables since the bounds on them are far less constraining than the ones studied in Table 1 in the foreseeable future, for the class of flavor models that we consider.
The ongoing experimental efforts in measuring eEDM and LFV observables will provide us with a trove of new data that can be used to probe various new physics models. The leading new physics contribution to all the observables in Table 1 arises from the dipole operators in the IR, where i, j indices refer to SM fermions, P R = (1+γ 5 )/2, F µν is the U (1) em field strength, and with c L ij being the Wilson coefficient of the operatorf i σ µν P L f j F µν . Any generic new physics model with new fermions or bosons can contribute to this operator through the diagram of Fig. 1. In App. A we review the contribution of generic new heavy scalars and fermions that can run in the loop to this IR operator. We explicitly check that our calculation is in agreement with the previous calculations in the literature [56][57][58]. Discussions of the renormalization group evolution of such operators may be found in Refs. [59,60], but we do not require this level of precision and so do not include them in our calculations.
In terms of these dipole operators, the eEDM is given by In addition to the eEDM observable, Eq. (3) contributes to the LFV observables such as BR (µ → eγ): where m µ is the muon mass and G F is the Fermi constant. The dipole operator in Eq. (3) can also give rise to conversion of µ to e in a nucleus. When this dipole operator is the only source of this conversion, we have [61] where the O(1) prefactor depends on the nucleus [61,62]. However, this observable can get contributions from other penguin and box operators as well. The dipole operator of Eq. (3) can also give rise to an enhanced BR (µ → 3e). This contribution dominates the four-fermion diagrams contribution and gives rise to [61] BR The contribution from all sub-dominant four-fermion operators to this observable can be found in Ref. [61] as well. Below, we will consider only the dipole-operator contributions to µ → eγ and µ → 3e.

SUSY breaking and R-symmetry
In this paper, we will assume a rather generic supersymmetric extension of the Standard Model. The breaking of SUSY can be parametrized through a chiral superfield X whose auxiliary field F X acquires a SUSY-breaking vev, F X /Λ ∼ m soft . For simplicity, we will assume that X couples in the most generic way possible through Λ-suppressed operators, unless they are forbidden by additional symmetries. If the cutoff is taken to be the Planck scale, this is the scenario generally known as gravity mediation. Any supersymmetric extension of the SM must confront two dramatic problems: the stability of the proton, and the mass scale of the higgs doublets (i.e., the µ-problem). We invoke a standard solution to both of these problems, namely a Z 4R R-symmetry under which the MSSM matter superfields Q,ū,d, L, andē carry charge +1 and the higgs superfields H u , H d carry charge 0 [63] (see also Refs. [64,65] for similar scenarios). We assume that the SUSY-breaking field X and the flavon-sector fields that we will introduce below also carry zero charge under the R-symmetry. As a result, Yukawa terms like d 2 θ y u H u Qū are allowed, but the dangerous µ-term d 2 θ µH u H d is forbidden. Forbidding a large superpotential µ-term allows for the Giudice-Masiero solution of generating the µ and b µ terms at the scale of soft SUSY breaking masses via the Kähler potential operators d 4 [66]. The Z 4R symmetry forbids not only the renormalizable baryon-and lepton-number violating R-parity odd operators H u L, LLē,ūdd, and QLd from the superpotential, but also the dangerous nonrenormalizable superpotential operators QQQL andūūdē that can induce proton decay. It allows unsuppressed, holomorphic contributions to the A-terms via the operators d 2 θ X Λ H d Lē, and to gaugino masses via the operators d 2 θ X Λ W α W α . In the future, it may be interesting to explore models with additional structure in the SUSY-breaking sector, especially those with mildly split SUSY that offer one of the simplest explanations for the observed higgs boson mass. For example, if X carries a charge under some symmetry, then gaugino masses and A-terms will be suppressed because holomorphic operators with a single X insertion will be forbidden. In this case, anomaly mediation could supply small gaugino mass terms [67,68], leading to a mildly split spectrum [69][70][71][72][73][74]. Similar phenomenology arises if the imaginary part of X is a periodic axion field, such that X comes equipped with a discrete shift symmetry. Kähler potential soft terms involving X + X † are still allowed, as is the gaugino mass term d 2 θXW α W α , provided that the coupling is adjusted to match the periodicity of Im X to that of the gauge theory's θ-angle. This symmetry suppresses A-terms, but also tends to suppress gaugino masses by a loop factor relative to scalar masses (see, e.g., Refs. [75]). Such structure is relatively common in scenarios with SUSY breaking mediated by moduli arising from extra dimensions, which sometimes have further structure modifying the form of the EFT [76][77][78][79][80][81]. For simplicity, in this paper we do not consider this range of possibilities. Instead, we will consider mediation only through a neutral SUSY-breaking field X, with couplings that are unsuppressed save for the additional spurions needed to account for flavor and CP symmetries, as we discuss below.

The SUSY CP and flavor problem
Supersymmetric extensions of the SM as described in Sec. 2.3 are generically severely constrained by flavor and CP violating observables. We intend to study the bounds on such models from the observables of Table 1. In such models sleptons, neutralinos, and charginos can run in the loop of Fig. 1 and contribute to the observables from Table 1. This can introduce stringent bounds on these particles' masses. 4 In this section, we review the basics of lepton flavor physics in the MSSM and argue that, without further structure, the bounds from the observables in Table 1 severely constrain its parameter space.
In the MSSM, the scalar (the heavy fermion) in the loop in Fig. 1 is replaced by various sleptons (charginos or neutralinos) to generate the dipole operators of Eq. (3) in the IR. Once the couplings between the various sleptons, charginos/neutralinos, and the SM fermions are determined, we can simply use the results of App. A to find the MSSM contribution to the dipole operator.
The chargino and neutralino mass eigenvalues are dominated by the gaugino masses and µ; see, for instance, Ref. [56] for their mass matrices. We work in a basis where these mass matrices are diagonalized. We assume the mixing between different mass eigenstates, which is proportional to the ratio of SM gauge boson masses to these mass eigenvalues, is small. This suggests that the gauge eigenstates and the mass eigenstates have substantial overlap, allowing us to refer to the mass eigenstates with their dominant component in the gauge basis, i.e., gaugino or higgsino.
While the MSSM model has many degrees of freedom, our qualitative conclusions are independent of the numerical relations between most of them. As a result, we use a few working assumptions about the MSSM spectrum to simplify our analysis. We list some of these assumptions below.
The gaugino masses are dominated by their soft masses M 1 and M 2 . It is usually expected that these soft masses are comparable. For simplicity, we assume the bino and wino soft masses are related as M 2 /M 1 = 2, as is often the case in models with gaugino mass unification; changing this ratio does not affect our final qualitative conclusions. In the SM fermion mass basis, the charged slepton mass matrices are parameterized as while the sneutrino mass matrices are simplyM 2 ν =m 2 δ LL . Herem is a shared mass scale for all the sleptons. The dominant part of the diagonal blocks δ LL/RR come from the SUSY breaking slepton masses, all of which are assumed to be of the orderm, while where M f is the SM leptons' 3 × 3 mass matrix, t β = H u / H d is the ratio of the two higgs vevs, µ is the usual mass scale in the MSSM higgs sector, and A is the leptons' Aterms matrix. 5 We also assume |µ| =m and A =mδ A , which suggests the higgsino masses and the scale of A-terms are comparable to the sleptons mass scale. Had we neglected the A-terms, in the SM fermion mass basis, δ LR would have been flavor-diagonal. All the δ matrices (δ LL , δ RR , δ LR , δ A ) have dimensionless entries whose magnitude is determined by the details of the UV model. Notice that unlike the charginos and neutralinos, in our working basis the sfermion mass matrix is not diagonal. The off-diagonal entries can give rise to LFV effects. The contributions of all sfermions and charginos/neutralinos to the dipole operators using the mass insertion approximation are included in App. B. The diagrams in the appendix are only meant to provide us with an intuition about the contribution of different states. In our numerical calculations in the next section, we properly diagonalize these mass matrices and calculate the dipole operators exactly.
Without any suppression of flavor-violating couplings, i.e., with arbitrary mixing pattern among sfermions in Eq. (9), the observables of Table 1 put very strong bounds on the mass scale of SUSY particles. The eEDM and BR (µ → eγ) in this scenario depend on the UV parameters as (see, e.g., [82,83]) Here, we assume all slepton masses are the same order asm and thatm ∼ |µ| M 2 ∼ M 1 . Off-diagonal elements of the δ A matrices are taken to be δ A ∼ 1 so that the A-term 5 Generically the mass scale of the LL block differs from that of the RR, at least due to different running from the scale of SUSY breaking. We neglect this O(1) effect for the sake of simplicity of the numerical study. Taking this effect into account typically eases the constraints from experiments as it decreases the L − R slepton mixing.
contributions are dominant for both the eEDM and BR (µ → eγ). The CP phases in δs are also assumed to be maximal, contributing to the eEDM.

Horizontal Symmetry and Spontaneous CPV
It is interesting to consider simple mechanisms that can suppress LFV and CPV effects. In the present section, we will see that introducing a new horizontal symmetry and forbidding explicit CPV can suppress these observables in the MSSM (see Eq. (12) and Fig. 2), thus opening up the parameter space for lighter supersymmetric particles.

A U (1) horizontal symmetry
A simple mechanism for suppressing the LFV effects in SUSY models is alignment [20,23]. In this mechanism, the slepton mass matrices are diagonal in the SM fermion mass basis, eliminating any sources of LFV. The most straightforward way to achieve this is by augmenting the model with a U (1) horizontal symmetry [18], i.e., a U (1)-augmented MSSM, under which superfields of different generation have different charges. As shown in Ref. [19], such a horizontal symmetry can not be exact. Generically, it is assumed that this symmetry is broken spontaneously via the vev of some "flavon". Although we do not expect UVcomplete theories to have exact global symmetries (spontaneously broken or not), a very good approximate U (1) symmetry could arise, protected by discrete gauge symmetries or perhaps as a remnant of an anomalous gauged U (1) obtaining a Green-Schwarz mass.
By charging the MSSM fields under this symmetry, we not only can generate the SM pattern of leptons' masses and mixings, but we can also suppress the disastrous LFV contributions from slepton loops. (The contribution of different superpartners are still captured by the diagrams in App. B.) In such a setup, the eEDM and BR (µ → eγ) observables depend on the UV parameters as [83] |d e | e ∼ 10 −29 10 TeṼ m where we use a shared slepton mass scalem withm ∼ |µ| elements of the soft mass squared matrices are suppressed by a horizontal symmetry as well, and the A-terms are neglected. We assume the fermion charges under this symmetry have the right values implied by the right SM fermion masses and lepton mixings (see Eq. (19) below). The off-diagonal elements δ LL,RR used here in the Yukawa diagonal basis are given in Eq. (23) below, and the CP phases are assumed to be maximal.
Comparing these equations to Eq. (11) indicates that, generically, with the right horizontal symmetry charges the disastrous contributions to the eEDM and BR (µ → eγ) can be suppressed compared to a structureless setup with random CPV entries in the mass matrices. In light of this, in what follows we introduce a general setup in which we augment the MSSM with a horizontal U (1) symmetry. We will study the bounds on this setup from the observables listed in Table 1.
Let us start by assuming we have a single flavon field S 1 . We can normalize all the charges such that this flavon's charge under the horizontal symmetry is [S 1 ] = −1 and that its vev's phase is rotated away using an appropriate field redefinition. For simplicity, we assume all the lepton superfields have integer charges under this symmetry too. Let us also assume the small expansion parameter S 1 /Λ ≡ 1 , where Λ is a UV cut-off where the details of the S 1 potential become relevant, and should not be confused with the SUSY breaking scale Λ SUSY . We also assume that S 1 is generated through supersymmetric dynamics, so that the flavons do not themselves mediate SUSY-breaking (and conversely, we will assume that SUSYbreaking spurions do not provide additional breaking of flavor symmetries). Conventionally, we assume 1 ∼ λ. 6 In addition to the horizontal symmetry that can suppress disastrous flavor-violating effects, it was argued in Ref. [33] that by prohibiting explicit CPV effects in the UV theory, light squarks could evade the stringent bounds from the neutron EDM. In such a theory, CP is assumed to be a gauge symmetry that is only broken spontaneously and via the interference of two flavon vevs, one of which has a physical phase. We now extend the idea of spontaneous CPV to the leptonic sector of the MSSM. In our setup we introduce another flavon field, S 2 , with [S 2 ] = −c 2 and S 2 /Λ ≡ 2 ∼ λ n 2 , with c 2 and n 2 integer numbers greater than 1. We assume S 2 gets a complex vev and spontaneously breaks CP.
Let us study in more details how spontaneous CPV can suppress different CPV effects in the MSSM. In such a setup, the only source of CPV is the interference between terms with different powers of the two flavons. For instance, if a holomorphic operator made of different superfields has the horizontal charge [O] = c O 0, the coefficients multiplying it can have the following flavon spurion insertions: An effective holomorphic superpotential term can be generated by inserting vevs in nonholomorphic Kähler potential terms, as in the Giudice-Masiero mechanism. In this case, the resulting term is suppressed by both SUSY-breaking spurions and flavor spurions. For non-holomorphic insertion of 1 or 2 , we have In writing these equations we used the fact that 1 = * 1 . It can be shown that the least λ suppression happens for terms where both flavon spurions 1,2 have non-negative powers, i.e., the holomorphic condition in Eq. (13) is satisfied. Moreover, from Eq. (13) we find that when c 2 = n 2 , different interfering terms have the same λ scaling, thus their interference can give rise to O(1) phases. This is the maximally violating CP scenario, as in cases with c 2 = n 2 the interfering phase is suppresed by powers of λ. Since we still need to generate the observed CPV in the CKM and PMNS matrices, we focus on the case with c 2 = n 2 and assume Arg( 2 ) = π/4, which maximizes the amount of CPV in our setup.
Nonetheless, even when c 2 = n 2 , spontaneous CPV can suppress various disastrous contributions to the eEDM. In particular, we consider cases where the higgses are uncharged under the horizontal symmetry. Thus, the µ term multiplies a neutral operator (c O = 0) and, according to Eq. (13), the holomorphic contribution to µ is real. The leading CPV contribution to µ from non-holomorphic Kähler potential terms is suppressed by c 2 which suppresses its contribution to the eEDM. It should be noted that this suppression crucially depends on the two higgses being neutral under the horizontal symmetry. We use this assumption in our model too.
With neutral higgses, the SM charged lepton mass matrix is then given by while the neutrinos' Majorana mass matrix is Here H u is the up-type higgs vev, i, j are SM fermion generation indices, C f /ν,ij a 1 ,a 2 are O(1) numbers, and M ν is an undetermined UV scale where new heavy neutrinos appear [84][85][86][87]. Notice that all C ij a 1 ,a 2 s are real; the only source of CPV is 2 . The hermiticity of the Lagrangian forces the undetermined numbers in the neutrino mass matrix to follow C ν,ij a 1 ,a 2 = C ν,ji a 1 ,a 2 . 7 7 In writing these equations we have neglected possible kinetic mixing terms among different superfields and start from a basis where the kinetic terms are canonical. As discussed in appendix B of Ref. [21], including these terms will merely give rise to some corrections to the undetermined O(1) coefficients in the mass matrices, i.e., slightly changes their distribution.
We then perform a rotation to the mass basis of the SM fermions. Similar to Ref. [18,23], we can show that the PMNS matrix entries scale like The large mixing angles among the LH leptons [88] suggests that various [L i ] should be close to each other. Comparing the fermion mass and mixing matrices with the observed values suggests the following relations between different fermions' charges: The first line is a result of the observed neutrino mixing pattern (V PMNS ), while the next lines are due to the observed charged lepton masses. Notice that thanks to undetermined O(1) numbers in the original mass matrices, we can still obtain the observed masses and mixing patterns even if we slightly deviate from these relations. The holomorphy of the SUSY models does not allow terms with negative powers of 1 or 2 to appear in the fermion mass matrices; this is enforced by the a i 0 condition above. The SUSY-breaking contributions to the sfermion mass matrices, on the other hand, can include non-holomorphic terms. The charged sfermion mass matrices are given by where the complex conjugate superscript ( * ) appears when a 2 0. Meanwhile, δ LR is given by where the last term captures the effect of A-terms from Eq. (10) and Recall that we assume all the A-terms have a shared scale, which is the same as the shared slepton mass scalem. We also assume |µ| =m. As mentioned before, we also assume M 1 = 0.5M 2 . Thus, apart from the horizontal charges of the superfields, the only free parameters in our setup are M 1,2 ,m, and tan β.
With a hierarchical choice of charges for superfields of different generations, the horizontal symmetry pushes us near the alignment limit for the sfermion mass matrix. Similarly, spontaneous CPV suppresses the eEDM contributions. This allows us to lower the scale of new SUSY partners in our setup, without violating the stringent LFV or eEDM bounds included in Table 1.
The C ij parameters in all the mass matrices originate from undetermined UV dynamics. Depending on their values, the model prediction for different IR observables can change significantly. As a result, in our study, we treat these numbers as random O(1) coefficients. For any fixed set of charges, we generate all the mass matrices many times and calculate the contribution of each trial to different IR observables. Further details on how we generate these mass matrices and the distribution of these random coefficients are included in App. C, where we also show the resulting distribution of the slepton masses and in particular show thatm is indeed the average slepton mass scale.
To illustrate the effect of these undetermined coefficients, in Fig. 2, for fixed values of the MSSM mass parameters, we show the contribution to the eEDM and BR (µ → eγ) in four different extensions of the MSSM. In each panel of the figure, we start in the SM fermion mass basis and generate the UV mass matrices for the sleptons 1000 times with different random O(1) coefficients reflecting our agnosticism about the UV model parameters. More details on the distribution of these random numbers are included in App. C. In the figure, we also fixm = 10M 1 = 10M 2 = 10 TeV 8 and t β = 5 in all the scenarios. The four scenarios we consider are summarized below.
• In the top left plot of Fig. 2 we consider a case with random complex numbers in all the slepton mass matrix entries and with no flavon spurion suppression. This is essentially a structureless (anarchic) MSSM model with explicit CPV. We assume Arg(µ) = π/4. Such a structureless mass matrix, where each entry is an O(1) number times a shared mass scale, does not typically give rise to mass eigenvalues that follow the SM hierarchical pattern of masses. 9 By starting from the SM fermion mass basis, we are neglecting this shortcoming of these structureless models in favor of a better comparison of their contribution to the eEDM and BR (µ → eγ) observables to the scenarios with horizontal symmetries. Eqs. (11)-(12); we go back to M 1 = 0.5M 2 assumption for the analysis of our models in the rest of the paper. 9 For that to happen, the random numbers in the charged leptons mass matrix had to roughly be spread in the range 10 −5 − 10 −2 . Assuming these random numbers are drawn from similar distributions (say a uniform distribution between −1 and 1), the probability of this outcome is exceedingly small. See the text for further details on each scenario. We observe that both the horizontal symmetry and spontaneous CPV can suppress the stringently constrained observables BR (µ → eγ) and eEDM; this opens up parameter space for lighter superpartners. Similar conclusions can be drawn using other LFV observables from Table 1 instead of BR (µ → eγ).
Eqs. (20)-(22) (without any 2 and 1 ∼ λ), for these charges we find times random complex numbers in each entry. Again we assume Arg(µ) = π/4. • Finally in the bottom right plot of Fig. 2 we consider a scenario with the same charges as the second scenario, but this time we introduce a second flavon with [S 2 ] = −2 (similar to the third scenario) and assume the only source of CPV is through its interference with the first flavon. This is a typical scenario with an horizontal symmetry and spontaneous CPV that we will be studying in more details in this work.
Comparing the first scenario with the second scenario, we see the substantial effect of the horizontal symmetry in suppressing both LFV and eEDM. The order of magnitude of the contribution of these two scenarios to the eEDM and BR (µ → eγ) are given by Eqs. (11)-(12), respectively. The contribution in the eEDM of the second scenario is dominated by the µ term phase, which is kept constant in different trials; hence, the distribution in the eEDM has a much smaller spread. The third scenario with only spontaneous CPV, and no further flavor structure, also shows a great suppression in the eEDM thanks to the suppressed µ argument. Finally, the type of models we are studying in this work are represented by the fourth scenario. We find a significant suppression in at least one of the two observables compared to all the previous scenarios.
The spread in the distribution of the observables shows the importance of undetermined O(1) numbers in the UV theory. As mentioned earlier, we treat these undetermined coefficients as random numbers that change in different trials. The sampling of the random numbers is further explained in App. C. Furthermore, this wide spread undermines the accuracy of naive order of magnitude estimations such as Eqs. (11)- (12). In the upcoming section, we develop a framework for interpreting the bounds on such models that turns this substantial effect of the undetermined UV parameters into a naturalness argument.

A new measure of naturalness for different observables
In the previous section, we highlighted the effects of a horizontal symmetry and spontaneous CPV in the MSSM that allow us to obtain the SM in the IR, while relaxing various LFV and eEDM bounds. However, a proper calculation of the contribution of these models to the IR operators is muddled owing to numerous undetermined UV theory coefficients in the mass matrices. As is evident from Fig. 2, these coefficients can give rise to orders of magnitude spread in the predictions from various flavor observables.
In particular, we argued that there is always some configuration of these random numbers that can completely suppress these observables. This suppression can happen for any scale of superpartner masses. As an extreme example, when all the off-diagonal random numbers in the slepton mass matrix of Eq. (9) are zero, the LFV observables do not get any contribution from this model, regardless of the scale of the slepton masses. Nonetheless, as we lower the superpartner scale, the fraction of space of the undetermined UV numbers that can give rise to suppressed values of these observables decreases.
In this section, we develop a framework for interpreting the bounds from different observables on models with horizontal symmetries while remaining agnostic about the undetermined coefficients in the UV theory. We can model our ignorance of these coefficients by randomly sampling the space of all these coefficients and see how natural it is for a given UV model to give rise to a viable set of IR observable values. One can view this from a Bayesian viewpoint: absent any strong prior reason to prefer one UV model to another, the models that are more likely (as we vary their parameters) to give rise to IR physics resembling what we see in our universe are more likely to be correct.
While, due to the effect of undetermined UV coefficients, we can not ever claim a U (1)augmented MSSM is definitely inconsistent with the results of a given set of precision experiments, we can calculate the fraction of the space of all the random numbers in which these observables are suppressed enough such that the model is not constrained by the aforementioned experimental results. This motivates us to propose the following measure of fine tuning in these models.
• Let's assume we identify the space of all the undetermined UV O(1) numbers that give rise to the right pattern of lepton masses and mixings as observed by the experiments and call its volume V SM . We then find the fraction of this space that further predicts the value of some observable O is in a certain range I and call its volume V O∈I . We take the ratio of these two volumes to define a measure of fine tuning associated with observable O. As an example, condition I could refer to the observable O being smaller than some fixed valueÔ, i.e., Let us clarify this definition further through an example. Let's assume we fix all the charges and mass parameters of a U (1)-augmented MSSM model and generate 10,000 different mass matrices for fermions and sfermions, each with their own unique random numbers multiplying different matrix entries. Let's assume we find that 100 of these trials give rise to the observed pattern of masses and mixings for leptons and that of these 100 "good" trials, 25 of them give rise to O Ô (O being some observable of interest). Here a trial refers to a set of fermion and sfermion mass matrices with its unique set of random numbers; see App. C for further details about distribution of these random numbers that we use in our study. Then t O Ô = 25/100 = 25%, i.e., 25% of the "good" trials gives rise to O Ô . The study of mass matrix models always suffers from the large number of free parameters. This complicates a rigorous study of the bounds on such models from IR observables, since the bounds can change significantly depending on what values are used for the undetermined free parameters. Our notion of naturalness enables us to sift through the entire parameter space of a model, a hitherto overlooked task that allows us to systematically update our priors on what models are preferred in light of experimental results.
Notice that this definition relies on a definition of a "good" trial. In our work we assume a trial is "good" if the following conditions are satisfied.
• The charged lepton masses it predicts are within a factor of λ of the observed lepton masses, • After we rescale all the neutrino masses to get the observed total neutrino masses [89] ν m ν = 0.12 eV, we find the neutrino mass differences to be within a λ of the observed values: In the last equation, the acceptable range of PMNS entries are comparable to the uncertainties in the measured value of each entry [90].
It should be noted that the ranges in Eqs. (25)- (27) are simply one possible choice; we could use similar criteria with slightly different acceptable ranges for these parameters. The exact numerical value we find for the tuning depends on these criteria for a good trial (and on the range we draw our random numbers from). However, by taking the ratio in Eq. (24), we suppress this sensitivity in the tuning calculation. For instance, one might have defined a "good" trial as one reproducing the measured masses and mixings within experimental uncertainties, but this is much less computationally tractable (and more of a time-dependent goal) and we expect that it would not substantially change our final results. 10 We have explicitly checked that changing the criteria for good experiments or the range from which we draw the random numbers at most gives rise to an O(1) change in the final value of the tuning.
In studying the bounds on our setup, we use the tuning measure t O O lim , where now O is (the absolute value of) any of the observables from Table 1 and O lim is the experimental limit on it. Notice that t O Ô is a non-decreasing function ofÔ. As a result of this, we can always go to lower values ofÔ by introducing more fine tuning. This suggests that we can always push down the sleptons to lighter masses without violating the existing experimental bounds from observables such as the eEDM or BR (µ → eγ) but the required tuning is worse.
In practice, we approximate the volume of these subspaces by sampling them multiple times and use the ratio of total number of trials in each sample to calculate t O Ô . To probe very small tuning values, we need to generate a large enough sample size of trials. Larger sample sizes not only allow us to probe lower-tuning subspaces, but also reduce the uncertainty in estimation of the tuning values.
A similar definition of fine tuning can be used in the context of any other U (1)-augmented model and for other flavor observables, e.g., the quark sector observables. Similarly to how we measure the naturalness of a given model from the perspective of the higgs mass hierarchy problem using some fine-tuning measure, we can use the fine-tuning measure of Eq. (24) to see how natural different models are from the perspective of various flavor physics observables.
It is also interesting to compare our naturalness measure to commonly-used naturalness measures of the higgs mass, such as the Barbieri-Giudice measure [91]. While in the latter the tuning measure quantifies the sensitivity of an IR observable (e.g., the higgs mass) to a certain UV model parameter at a given point in the parameter space of the model, our tuning measure from Eq. (24) quantifies how naturally the experimental results, here the experimental results on a given observable O, emerge from the undetermined UV parameters of the model as a whole. An analogue of the Barbieri-Giudice measure might be used when we fix all the undetermined UV parameters to some value and study the sensitivity of the observable at that specific point in the parameter space. We, however, want to study the entire parameter space of a model and see how typically the SM and other IR observables emerge from that setup. We aim to compare different models, rather than individual points in the parameter space of a single model.
It will be interesting to apply our analysis to any BSM model to explain the SM flavor structure and see how much tuning should be introduced in order to satisfy various flavor physics constraints at a given mass scale of that model.

Fine tuning and flavor bounds on the MSSM
In this section, we apply our tuning measure to the U (1)-augmented setup introduced in Sec. 3.1 to see how natural the light superpartners are given the current and future bounds on the observables of Table 1. The specific question we want to answer is: for a given U (1)-augmented MSSM model (with fixed charges under the horizontal symmetry and mass scales), what tuning values should be tolerated in order to avoid the existing constraints on the LFV observables and the eEDM?
We now consider a handful of different U (1)-augmented MSSM models and, for various amounts of tuning, calculate the bounds on their parameter space from LFV observables and the eEDM. We scanned over different charge assignments (assuming integer charges, and |[f ]| 4 for any superfield f ) and tan β values. For any reasonable tan β, we found on the order of ∼ 1000 different charge assignments that can give rise to the right pattern of lepton masses and mixings (with the criteria in Eqs. (25)- (27)) at least once from among 1000 trials. We use the same assumptions on the mass spectrum as outlined in Sec. 3.1. To review, we assume that all the A-terms and the sleptons masses and |µ| are proportional to the same scalem, M 2 = 2M 1 , and the µ phase is given by Eq. (15).
In this work we go a step beyond the common practice in the literature. Instead of studying an arbitrary single charge assignment, we study five (out of many) charge assignments which span a wide range of t β . We show these charge assignments in Table 2. It should be noted that these are merely sample charges and it is possible that by scanning over the entire set of charges, we could find charge assignments that are less constrained. Table 2: A few sample charge assignments that give rise to SM in the IR quite often. A total of 5 × 10 5 trials are generated for each charge assignment. In the last column we indicate the fraction of these trials that give rise to the observed pattern of masses and mixings among the leptons, i.e., "good" trials as defined in Sec. 3.2. We do find some flexibility in the conditions from Eq. (19) on the charge assignments; yet, by and large, these conditions give the right ballpark of the charges needed to get a non-negligible number of good trials.
Our study is an intermediate step towards an exhaustive study of large number of charge assignments, which requires the development of novel scanning methods. 11 The study of a large number of charge assignments emphasizes which features are a universal outcome of horizontal symmetries and which features are unique to a specific charge assignment.
For the charges we study in Table 2, the charges in the first and fourth rows of the table obey the exact relations predicted by Eq. (19). 12 In the second, third, and fifth rows, the relation between the LH fermion charges slightly deviate from Eq. (19).
For each setup, we scan over the common slepton massm as well as gaugino masses M 2 = 2M 1 . 13 For each point in the mass plane, we generate N tot = 5 × 10 5 different trials.
In the last column of Table 2 we report the fraction of trials from each setup that gives rise to the right pattern of lepton masses and mixings as reviewed in Sec. 2 and in Sec. 3.2. We find that the first and fourth row's charge assignments, which follow the relations in Eq. (19), give rise to a slightly better "good" trial ratio than others. Yet, other charge assignments with some deviation from these relations still give rise to SM-like pattern of masses and mixings in the IR with comparable frequency. It is worth noting that getting the SM masses and mixings in the IR is itself quite rare; we only find a sub-percent fraction of the trials that give rise to the SM masses and mixings in the IR. Nonetheless, the U (1)-augmented models give rise to SM-like masses and mixings far more often than a completely anarchic model would, where merely getting a sufficiently small electron mass (let alone all the other masses and mixings) would require a ∼ 10 −5 tuning cost. 11 See Refs. [92,93] for recent proposals for efficient methods of finding all potentially suitable charge assignments. 12 Notice the third row is an anarchy scenario in the neutrino mass matrix [39], favored by the observations from various neutrino experiments. 13 Our calculation can straightforwardly be repeated for other MSSM spectra with different relationships between various mass scales.
We then calculate the eEDM and the LFV observables of Table 1 for each of these trials. With these numbers, we can calculate t O Ô for different values ofÔ for different observables. This, combined with the experimental bounds from Table 1, allows us to find the bounds on the slepton-gaugino mass plane for different fine tunings.
We repeat this calculation for the five charge assignments in Table 2 and for different SUSY parametersm and M 2 = 2M 1 . The results are shown in Figs. 3-7. 14 In the top (bottom) row of each figure we show the current (future) bounds on the model; see Table 1. We report the bounds on the models with different values of tuning allowed. The conclusions derived about the model from Figs. 3-7 are more or less similar, so we discuss the main takeways only from Fig. 3; the discussion for other models would be closely analogous.
In the upper panels of Fig. 3, we find that if we want to avoid any significant tuning in the model, i.e., t ≈ 50%, and if all the particles are comparable in mass, the current bounds on sleptons and charginos/neutralinos from different observables are in the 10 -50 TeV ballpark. For TeV scale superpartners to avoid the current bound from each observable of Table 1, a sub-percent-level tuning is required in this setup, and an even higher level of tunning is needed in order to push the superpartners masses all the way down to the current bounds from the LHC [94,95]. Note however that we present the average slepton mass scale, from which the mass eigenstates might differ by an O(1) factor. Without any correlation, the values of tuning are multiplied over all the observables of Table 1 to find the total degree of tuning required in the model with a given mass scale.
As discussed in Sec. 3.1, the assumption of c 2 = n 2 for the S 2 flavon charge and its spurion's λ scaling guarantees that we have the maximal CPV possible in our setup. Even so, we still find that, with the current experimental results summarized in Table 1, BR (µ → eγ) gives rise to slightly more stringent bounds on the parameter space than the eEDM, while the current BR (µ → 3e) bounds are comparable to the eEDM bounds for this charge assignment. With the projected bounds of Table 1, we find that the bounds from all these observables are comparable.
Similar deductions are made for the other charge assignments studied in Figs. 4-7. While the bounds on superpartners' masses can change by O(1) numbers across different charge assignments, the final conclusions are similar for different charge assignments studied: that different models' superpartners are currently ruled out all the way up to O(10) TeV masses if the model is completely natural, and to realize O(1) TeV scale (or even lighter) superpartners requires fine-tuning. 15 The similarity of the results emphasizes that the observed patterns   Table 1) from eEDM (top), BR (µ → eγ) (middle), and BR (µ → 3e) (bottom), with different amount of tunings. A total of N = 5 × 10 5 trials are generated and only the good trials are used in calculating the bounds. The left (right) column shows the current (future) bounds from each observable. We find that for this model to be natural, i.e., tuning of ∼ 50% for each observable, the superpartners should be heavier than ∼ 10 − 50 TeV; around percent-level tuning for each observable is required to realize superpartners of mass O(1) TeV.   Table 1). are a general feature of horizontal symmetries and not of a unique charge assignment. It raises the question of how one can study a more complete set of good charges. Our tuning measure allows one to study different charge assignments and find explicit bounds on its parameter space from these flavor observables (as well as other observables neglected in this study). An exhaustive study of charge assignments might uncover a unique preferred set of charge assignments, as well as exposing universal features of horizontal symmetries.
We should also reemphasize that a similar calculation can be done for any model with undetermined UV parameters and for other observables. In particular, it will be interesting to repeat this analysis for U (1)-augmented models (either supersymmetric or not) in the quark sector. We leave such a study for a future work.
In calculating the tuning in this section we studied one observable at a time and considered the fraction of the UV parameter space that suppresses that observable. Instead, we could look at the subspace that simultaneously gives rise to a small enough value for different observables, which enables us to take account of correlations between different observables. We will leave a study of such correlations for a future work.

Conclusion
The MSSM remains an extremely appealing candidate for a UV completion of the SM. The minimal version of the theory predicts disastrous contributions of superpartners to various LFV observables and the eEDM. In this work we focused on a U (1)-augmented MSSM, and carried out a complete calculation of the contribution to these observables from different sleptons and charginos or neutralinos. We showed that the combination of a horizontal symmetry and spontaneous CPV can suppress the contributions to the full set of leptonic observables.
The existence of undetermined O(1) numbers in mass matrix models always muddles the interpretation of the various experimental bounds on these models. To properly interpret the bounds from various observables on such models, we developed an intuitive notion of naturalness. We treat the undetermined UV coefficients as random parameters in the mass matrices and define the tuning measure as the fractional volume of the space of all these possible coefficients that satisfy a given condition on an observable, while giving rise to the SM patterns of masses and mixings in the IR. In particular, for any flavor observable we can use the condition that the model gives rise to a suppressed contribution to that observable. This provides us with a measure of how natural various mass matrix extensions of the SM are, considering the existing experimental results on any specific observable. We applied this naturalness notion to our study of the U (1)-augmented MSSM.
In all the charge assignments we studied, for a fixed value of tuning, the current bounds       Table 1). from LFV observables were stronger than those from the eEDM. We found that if we want to avoid any significant tuning in the model, the current bounds on sleptons and charginos/neutralinos from different observables are in the 10 -50 TeV ballpark. It is conceivable that by introducing more structures in the UV model, the TeV scale superpartners can become viable with mild tuning. Furthermore, the supersymmetric setups studied in this work also contribute to the muon g − 2 via the c R 22 Wilson coefficient; this motivates a future repeat of our tuning analysis with muon g − 2 included in the list of observables.
We studied a handful of sample charge assignments, in contrast to the common practice in the literature of studying a single charge assignment; nevertheless, there are many more charge assignments that can give rise to SM-like theories in IR with comparable frequency. We emphasize that the development of novel methods is called for, in order to efficiently study all such charge assignments.
We can repeat our analysis for different charge assignments or, more generally, for any mass matrix model. Our naturalness notion enables a new rigorous way for interpreting the bounds from various flavor physics observables on such models in general. It also puts various flavor physics observables on the same footing as the higgs mass and allows us to compare the naturalness of different flavorful BSM models. While we only focused on a handful of lepton flavor observables, our analysis can be repeated for the quark sector observables in any mass matrix model as well. After integrating out the new heavy particles and matching the UV diagram onto the dipole operator, we find where m φ (m χ ) is the heavy scalar (fermion) mass while Q φ (Q χ ) is its electric charge, and r = m 2 χ /m 2 φ . The loop functions are given by A(r) = 2r 3 + 3r 2 − 6r + 1 − 6r 2 log(r) 6(r − 1) 4 , B(r) = r 3 − 6r 2 + 3r + 2 + 6r log(r) 6(r − 1) 4 .
We have explicitly checked that our formulas are in agreement with the corresponding equations from Refs. [56][57][58]. As evident from Eq. (29), the A(r) and B(r) loop functions correspond to diagrams with a heavy fermion mass insertion in the loop, whileĀ(r) andB(r) correspond to diagrams with a mass insertion on the external fermion legs. The terms with A(r) andĀ(r) (B(r) andB(r)) come from diagrams with the photon emitted from the scalar (fermion) in the loop. Evidently, only the term proportional to y L y R can contribute to the CP-odd observables.
B Diagrams Generating the Dipole Operator in the

Mass Insertion Approximation
As mentioned in Sec. 2, there are many diagrams in the MSSM that can contribute to the dipole operator of Fig. 1. In the absence of RPV, we only have diagrams with sleptons and higgsinos or gauginos running in the loop. To properly calculate the contribution of all these diagrams, one should go to the mass eigenbasis of these particles and use the result of the previous appendix with the new physical states and their couplings (which can be found in Ref. [56]). We use this basis to carry out all the calculations in deriving our numerical results. In this appendix, we provide a better analytic understanding of these contributions to c R ij using the mass insertion approximation. We have checked explicitly that these approximations are in great agreement with our full numerical calculation as well as those obtained from publicly available code susy flavor v2.5 [96][97][98].
We choose a common scalem for all the sleptons, such that their mass eigenvalues are comparable. We keep the original non-diagonal form of the slepton mass matrices and treat the δs from Eq. (9) as perturbations on the propagators. Since the δs all depend on the horizontal charges of the superfields, this calculation allows us to study the effect of varying different charges.
The gauginos and the higgsinos, on the other hand, can have vastly different masses. To keep track of their contributions, here we go to their mass eigenbasis. This gives rise to small modifications of their couplings to leptons and sleptons, proportional to the mixing parameter between the gauginos and the higgsinos. This mixing parameter can be treated as a perturbative parameter in our parameter space. We denote this parameter by ε n (ε c ) in case of neutralinos (charginos). Its smallness suggests that the mass eigenstates and the gauge eigenstates overlap substantially.
We show the leading order contributions to c R ij from various particles in the loop in this basis in Tables 3-6. We consider all possible combinations of fermions and sleptons in the loop. In each row we show the fermions in the loops and the SM coupling that appears in each vertex. (Other factors of O(1) at the vertices are neglected; we also denote both g 1 and g 2 gauge couplings of the SM by g for simplicity.) The charge and the chirality of the slepton running in the loop can be determined from this information.
There are three different perturbative parameters that help us identify the most dominant diagrams. These are (i) the SM fermion yukawa couplings y µ and y e , (ii) the neutralinos (chargino) mixing parameters ε n (ε n ), and (iii) the inverse of the common slepton soft mass 1/m 2 , which can be introduced through the loop functions, as well as through δ LR s (see Eq. (10)). Thus, we interchange this last expansion parameter with δ LR . For simplicity, we denote the mixing between all neutralinos (either the two gauginos with each other or one gaugino and one higgsino) by ε n . In each row we consider the diagram with the least number of mass insertion δ and at most one δ LR entry insertion. These δs are the source of CPV (in Table 3) and LFV (in Tables 4-6) 16 . We emphasize that, depending on the flavor 16 Notice that the hermitian property of δ LL and δ RR forces us to have at least three δ insertions in contributions to the eEDM in Table 3 when the slepton chirality is preserved in the loop. model, there are cases in which more δ insertions are actually the leading term. As our tables capture all terms with the correct flavor and chirality structure, it easy to obtain all higher order terms by using a replacement of the sort: In Table 3 we show the diagrams contributing to the imaginary part of c R 11 , which is relevant for calculating the eEDM. As a result, only the couplings of the form y R y L from Eq. (29) are relevant. In Tables 4-6 we show the contributions to c R 21 that respectively come from the terms proportional to y R y L , y L y L , and y R y R from Eq. (29). Results similar to those shown for c R 21 in Tables 4-6 can be obtained for any c R f f with trivial replacements. 17 In Tables 4-6, we include two different δ insertions in the rows with δ LR . The first (second) insertion denotes the contribution of the A-terms (the µ term) from Eq. (22). 18 For these terms, our organizing principle (including the diagrams with the fewest number of δ insertions) could actually miss the dominant contributions to the LFV observables (for a given χ, c 1 , and c 2 ). For instance, in our setup, in the contribution from the third row of Table 4, the term with δ LR µµ × δ RR µe can dominate the δ LR µe contribution thanks to the factor of tan β in Eq. (21). 17 It should be noted that then the c L ij coefficients can be calculated using Eq. (4). 18 We should keep in mind that the A-term contributions to δ LR can be non-diagonal in the SM fermions mass basis, while the contribution from the first term in Eq. (21) will be proportional to M f , i.e. diagonal in the SM fermions mass basis.

The Diagram
W 0 ε n g ε n y e δ RR µẽ W 0 ε n y µ ε n y e δ LR µe , δ LR µµ × δ RR µẽ W 0 ε n y µ g δ LL µe Table 4: Diagrams with sleptons and gauginos/higgsinos in the loop contributing to c R 21 that come from the y L y R term in Eq. (29). This Wilson coefficient contributes to LFV observables such as BR (µ → eγ). The general topology of the diagram is shown on left. The second column shows the gaugino or the higgsino in the loop. The SM couplings appearing in each vertex, as well as εs that characterize neutralino or chargino mixings, are shown in the next two columns. We consider all possible combinations of couplings with all possible fermions and sleptons in the loop. The last column includes the δ insertions on the slepton propagator. We are only including diagrams with the fewest number of δ insertions and at most one δ LR entry. These δ insertions are the source of LFV in our setup. As evident from the table, whenm is not too large, i.e. δ LR is not too small, the diagram with a Bino in the loop and a chirality flip on the propagator dominates, while other diagrams have more suppression thanks to fermion mixings ε n/c or the SM yukawas.

C Generating Random Mass Matrices
As explained in Sec. 3, we study the effect of undetermined coefficients in the mass matrices by generating the mass matrices multiple times, with random O(1) numbers replacing these coefficients. The distribution of these random numbers gives rise to a spread in the predictions for various observables across different trials. In this appendix, we elaborate more on generating these random numbers and how we introduce them into the δ LL,RR,A matrices. We emphasize that different prior distributions give rise to different results. Therefore one has to make a choice of some reasonable distribution. In what follows we explain and motivate our choices. Our philosophy is to choose the simplest distribution that is close to a uniform distribution but also satisfies key physical constraints like hermiticity or positivity.
We start by considering the four scenarios considered toward the end of Sec. 3.1. This is the procedure used to prepare Fig. 2. In each scenario, we assumed that we start from the SM leptons' mass basis. It should be noted that in the notation of Sec. 3.1, each entry of the δ matrices can get contributions from different interfering spurion insertions. Each of these interfering terms can be multiplied by their own random coefficients (hence the subscripts a 1 and a 2 in C ij a 1 ,a 2 in Eqs. (16)- (17) and Eqs. (20)- (21)). As a result, each δ matrix can be broken into two pieces: (i) a 3 × 3 matrix with random entries denoted by C L,R,A , and (ii) a piece containing all interfering spurion insertions. (The latter only appears in the scenarios when we do have horizontal symmetries.) We then multiply these two matrices entry-by-entry. This is repeated for each of the interfering terms separately. We finally add up these interfering matrices to make the final δ LL , δ RR , and δ A matrices.
It should be noted that while δ A can be a generic 3 × 3 matrix, δ LL and δ RR must be hermitian and positive definite. (The requirement that they be positive definite follows from the phenomenological need to avoid tachyons that would break electromagnetism, not from theoretical consistency.) In the horizontal symmetry scenarios the spurion part of δ LL and δ RR are hermitian matrices that are very close to the identity matrix, thus we can guarantee the hermitian and positive definite properties of the overall δ LL and δ RR by making sure that C L and C R are hermitian and positive definite. 19 To do that, in scenarios with spontaneous CPV and for each interfering term, we draw 9 random real numbers each from a uniform distribution between (−1, 1) and call the resulting matrices X L,R . Then we define C L,R = X L,R · X L,R † . (Note that while X L,R are uniformly distributed, C L,R are not. Nevertheless, this construction guarantees the physical constraints are satisfied by 19 The δ LR matrix could potentially spoil this property; however, its entries are suppressed compared to δ LL,RR matrices by roughly H d /m, which can be a small number. As a result of this, we can guarantee the slepton mass matrices are positive definite and hermitian by only focusing on their diagonal blocks δ LL,RR . ����������� �  = �� � ��� Figure 9: The distribution of different slepton masses across different good trials for the charge assignment in the first row of Table 2 withm = 100 TeV. The orange histograms show the distribution of different mass eigenvalues, while the black histogram denotes the average slepton mass in different trials. Note that the average masses are sharply concentrated aroundm. While the heaviest eigenvalues are within an O(1) factor ofm, the lightest eigenvalue for this charge assignment can be, in rare occasions, more than an order of magnitude smaller than the average slepton masses.
C L,R .) We also generate 9 more random numbers uniformly distributed between (−1, 1) and directly use them for each interfering term in C A . In Fig. 9 we show the distributions of the slepton mass eigenstates for an arbitrary point of the parameter space. We find that the average slepton mass agrees withm at the percent level, while the individual mass eigenstates typically differ from the average by an O(1) number.
For the scenarios with explicit CPV, the only difference is that the random numbers will have a magnitude and phase that are drawn from a uniform distribution in the range (0, 1) and (−π, π), respectively. In each scenario of Fig. 2, we repeat this process 1000 times. In each trial, we calculate the mass eigenstates and their contribution to the eEDM and BR (µ → eγ) observables.
A similar process is repeated for our calculations with the charge assignments of Table 2 in Sec. 3.3, which were used to prepare Fig. 3, 4, 5, 6, and 7. Unlike the scenarios in Sec. 3.1, in these models we start from an arbitrary basis with off-diagonal entries in the SM fermion mass matrices. We generate the random matrices X l,ν for the charged leptons and neutrinos mass matrices, as well as X L,R,A that will eventually appear in the corresponding δ matrices.
Physical consistency requires the neutrino mass matrix to be symmetric. Thus, we use C ν = (X ν +(X ν ) T ) and C R are also made positive definite and hermitian as prescribed above. All the random numbers are still drawn from a uniform distribution between (−1, 1). We then rotate to the SM fermion mass basis. Then, similar to the scenarios in Sec. 3.1, we go to the mass basis of the new SUSY particles and calculate the dipole operator and the observables of Table 1 for each trial. Note that depending on the distribution from which we draw the random numbers, the fraction of "good" trials in Table 2 can change. Using the uniform distribution reflects our complete lack of knowledge about the UV structure giving rise to the undetermined coefficients. While the good trial fraction can be affected by this change of distribution, it is natural to think the tuning ratio of Eq. (24) is far less sensitive to this distribution. This is why we choose to use this ratio as our measure of tuning. It not only captures the constraints from a specific observable (rather than getting the SM mass and mixing structure), but also partially cancels the sensitivity of the final tuning value to the distribution of these numbers in the numerator and the denominator.