Observable flavor violation from spontaneous lepton number breaking

We propose a simple model of spontaneous lepton number violation with potentially large flavor violating decays, including the possibility that majoron emitting decays, such as μ → eJ, saturate the experimental bounds. In this model the majoron is a singlet-doublet admixture. It generates a type-I seesaw for neutrino masses and contains also a vector-like lepton. As a by-product, the model can explain the anomalous (g − 2)μ in parts of its parameter space, where one expects that the branching ratio of the Higgs to muons is changed with respect to Standard Model expectations. However, the explanation of the muon g − 2 anomaly would lead to tension with recent astrophysical bounds on the majoron coupling to muons.


Introduction
Neutrino masses may be non-zero due to the violation of lepton number, L, in which case neutrinos are Majorana particles. The literature is abound with Majorana neutrino mass models which simply add explicit lepton number violating interactions or mass terms to the Lagrangian, but the violation of L could well be spontaneous in origin. The spontaneous breaking of a global continuous quantum number generates a massless Goldstone boson, in the case of lepton number usually called the majoron [1,2], J. The original majoron of [1,2] is a complete gauge singlet. On the other hand, it is also possible to use larger multiplets to break lepton number spontaneously, resulting in the doublet and triplet majorons [3][4][5]. In general, the majoron can be an admixture of these three representations, or even more exotic possibilities.
A massless boson in the particle spectrum certainly will affect phenomenology. However, whether these changes with respect to the explicit lepton number violating models are quantitatively important depends very strongly on the model and the nature of the majoron. The pure singlet majoron interacts so weakly with all Standard Model (SM) particles, that it is highly unlikely it will ever be observed experimentally. In the other extreme, pure doublet and triplet majorons have been ruled out by LEP [6], but majorons with sufficiently large singlet admixture can escape this constraint.

JHEP01(2022)098
The interaction of the majoron with charged leptons can be described in a model independent way as [7], L J = J¯ β S βα L P L + S βα R P R α + h.c. , (1.1) where α,β are the standard light charged leptons and P L,R are the usual chiral projectors. The S L,R couplings are dimensionless coefficients, and we consider all flavor combinations: βα = {ee, µµ, τ τ, eµ, eτ, µτ }. Due to the pseudoscalar nature of majorons, the diagonal S ββ = S ββ L + S ββ * R couplings are purely imaginary. Apart from the invisible Z-boson width, measured at LEP, there are a number of laboratory and astrophysical constraints on majorons. First, neutrinoless double beta decay will occur with the emission of a majoron. The best current constraints come from the EXO-200 experiment [8], limiting the effective coupling of the majoron to electron-type neutrinos, to values roughly below O(10 −5 ). Astrophysics also constrains the majoron parameter space. The production of majorons inside stars and their posterior emission constitutes a very efficient stellar cooling mechanism. This allows one to set very stringent constraints on the majoron couplings to charged leptons. For instance, white dwarfs allowed the authors of [9] to set the bound Im S ee < 2.1 × 10 −13 , whereas ref. [10] considered the supernova SN1987A and found Im S µµ < 2.1 × 10 −10 . Last but not least, one can also derive indirect bounds on the majoron couplings to charged leptons from the bounds on the majoron couplings to photons, since the former induce the latter at the 1-loop level. Using results from the OSQAR experiment [11], a light-shining-through-a-wall experiment, ref. [7] found the approximate bounds S ee 10 −7 and S µµ 10 −5 . However, these are less stringent than the stellar cooling bounds mentioned above.
In this paper we study charged lepton flavor violation (LFV), connected to the spontaneous violation of lepton number. While our numerical results are obtained for one particular model realization, it is possible to discuss some general features qualitatively. The S L,R couplings in eq. (1.1) can be generically written as Here, Y is a general matrix in flavor space and M is the diagonal charged lepton mass matrix, M = diag(m e , m µ , m τ ). The coefficients a and b depend on the model under consideration. In the case of a type-I seesaw and a pure singlet (or doublet) majoron, both a and b are zero at tree-level, but generated at 1-loop [2,12,13]. The 1-loop diagram is further suppressed by the small mixing between the right-handed neutrinos and the active states, thus one expects that decays such as µ → e J are unobservable. For the triplet majoron, on the other hand, a is non-zero at tree-level. It is generated via the mixing of the triplet with the SM Higgs, due to a coupling of the form λ σH∆H, where σ is the scalar singlet, H the SM Higgs and ∆ the triplet scalar. b again can be generated only radiatively and LFV interactions with majorons are expected to be tiny. How about type-III seesaw? In the type-III seesaw, the charged leptons of the SM are mixed with the charged components of the fermionic triplet Σ. In the spontaneous version of this setup JHEP01(2022)098 one would thus expect some non-diagonal coupling of the majoron to charged leptons to appear in the mass eigenstate basis [14]. However, the corresponding mixing is related to the small neutrino masses and thus, in the end, the rates for α → β J decays are typically very small.
This discussion can serve as a basis to establish the criteria a model has to fulfill in order to have sizeable off-diagonal couplings between the majoron and charged leptons. First of all, if the majoron couplings to charged leptons are induced via majoron mixing with the SM Higgs doublet, they will be exactly diagonal in the charged lepton mass basis. Therefore, in order to obtain sizable off-diagonal couplings, the majoron must couple directly to the charged lepton sector. This coupling can be either to the light charged leptons themselves or to some additional heavy charged leptons which, after symmetry breaking, mix with them. In the first case, sizeable off-diagonal couplings can be obtained with a non-universal lepton number assignment, whereas in the second case it is crucial that the light-heavy mixing is not suppressed by neutrino masses. We refer to [15] for a discussion along similar lines.
In this paper we propose a relatively simple model that induces large off-diagonal majoron couplings to charged leptons at tree-level. Our model adds a singlet and a doublet scalar to the SM, both with lepton number. It also extends the SM symmetry by imposing lepton number conservation. Finally, we introduce three right-handed neutrinos and a vectorlike lepton, which we choose to be an SU(2) L singlet for simplicity. After the electroweak and lepton number symmetries are spontaneously broken, neutrinos acquire non-zero masses via a TeV-scale type-I seesaw mechanism and the vector-like lepton mixes with the SM charged leptons, inducing in this way large LFV majoron couplings. Furthermore, extending the SM lepton sector with vector-like fermions will affect a number of observables, most notably the anomalous magnetic moment of the charged leptons, as well as their coupling to gauge bosons. As a by-product of our construction, the model can explain the observed anomaly in the muon anomalous magnetic moment [16,17] in parts of its parameter space. It can also lead to observable effects in Higgs boson decays, most notably in h → µµ.
We mention in passing that branching ratios for µ → e J decays as large at the experimental limit have been found in [18]. The underlying model is supersymmetric with spontaneous violation of R-parity. This provides one, albeit quite complicated, model example, where off-diagonal majoron couplings to charged leptons are induced at tree-level and can be large. Nevertheless, we note that this decay can in principle saturate the experimental bound also in models that generate the majoron couplings to charged leptons at the 1-loop level [13].
The rest of this paper is organized as follows. In the next section we introduce the model, discuss the majoron profile and the mass matrices for leptons and scalars. Section 3 is devoted to a discussion of the majoron couplings, while section 4 discusses the possible phenomenological signatures. In section 5 we present our numerical results, before closing the paper with a short summary in section 6. Some more technical aspects of our calculations are given in the appendices. 3  3  3  3  3  3  1 1 1 1

The model
We consider a type-I seesaw with spontaneous lepton number violation. The quark sector remains as in the SM, whereas the lepton sector is extended with the addition of 3 generations of singlet right-handed neutrinos, a pair of singlet vector-like leptons, F L and F R , the scalars σ and S and a U(1) L global symmetry, where L refers to lepton number. The full particle content of the model and the representations of all fields under the gauge and global groups are shown in table 1.
As usual, the SM SU(2) L doublets can be decomposed as whereas the new S doublet can be decomposed as Under the above working assumptions, the most general Yukawa Lagrangian allowed by all symmetries can be written as is the usual SM Lagrangian, with H = i τ 2 H * and Y u,d,e 3 × 3 Yukawa matrices in flavor space. We have omitted SU(2) L contractions and flavor indices to simplify the notation.
The new terms are given by Here Y ν and κ are 3 × 3 matrices, ρ is a 3 × 1 matrix and Y S is 1 × 3 matrix. The parameter M F has dimensions of mass and, again, SU(2) L contractions and flavor indices have been omitted for the sake of clarity. The guiding principle when writing eq. (2.5) is

JHEP01(2022)098
the conservation of lepton number. 1 Finally, the scalar potential of the model also includes new terms involving the σ and S fields. It can be written as with φ = H, σ, S, and (2.8) Here all m 2 φ parameters have dimensions of mass 2 , whereas µ is a parameter with dimensions of mass.

Scalar sector
The scalars of the model take the vacuum expectation values (VEVs) These relations define the VEVs v H , v σ and v S , which break the electroweak and U(1) L symmetries. As a result of this, the W and Z gauge bosons acquire non-zero masses, given by where v 2 = v 2 H + v 2 S and g and g are the SU(2) L and U(1) Y gauge couplings, respectively. Here v 246 GeV is the usual electroweak VEV, which receives contributions from both scalar doublets H and S. The tadpole equations obtained by minimizing the scalar potential read 14) The trilinear µ term in eq. (2.8) plays an important role. In the limit µ → 0, the scalar potential has an accidental U(1) symmetry, under which all scalar fields can have arbitrary JHEP01(2022)098 charges. Therefore, its presence is crucial to explicitly break this global symmetry and avoid the appearance of an unwanted Goldstone boson. The µ term also induces a tadpole for each of the scalar fields if the other two have non-vanishing VEVs. Assuming that CP is not violated in the scalar sector, namely that all scalar potential parameters and VEVs are real, we can split the neutral scalar fields into their real and imaginary components as The scalar potential contains the piece V mass = V N mass +V C mass , with mass terms for the neutral (V N mass ) and charged (V C mass ) scalars in the model. The neutral scalars mass terms read where z = {H 0 , σ, S 0 } and M 2 R and M 2 I are the 3×3 squared mass matrices for the CP-even and CP-odd neutral states, respectively. The prefactors of 1 2 are due to the fact that Re(z i ) and Im(z i ) are real scalar fields. It is straightforward to get the analytical expressions of the mass matrices, which can be computed as (2.20) Then, using the expressions above, one finds 2

21) and
It proves useful to compute the matrix M 2 I in a general R ξ gauge, since this allows for the proper identification of the Goldstone boson that becomes the longitudinal component of the Z boson. However, we present here the results in Landau gauge (ξ = 0).

JHEP01(2022)098
One can now use the tadpole equations in eqs. (2.12)-(2.14) to evaluate these matrices at the minimum of the scalar potential. We obtain and (2.24) The physical CP-even mass eigenstates The model has thus 3 CP-even neutral scalars. One of them, presumably the lighest, is to be identified with the Higgs boson discovered at the LHC, H 1 , with m H 1 ≈ 125 GeV. Similarly, diagonalizing the mass matrix M 2 I , we can obtain the profile of the three CP-odd mass eigenstates. We end up with a massive state, that we denote by A, and two massless states. One of the massless states is the Goldstone boson, z, eaten by the Z boson, while the other is the majoron, J, the Goldstone boson associated to the spontaneous breaking of lepton number. In the basis Im {H 0 , σ, S 0 } = {P H , P σ , P S }, the mass eigenstates are given in terms of the original gauge eigenstates as with their masses given by As already discussed, the µ parameter breaks an accidental U(1) symmetry that would lead in its absence to the appearance of an additional massless Goldstone boson. This can be observed in m 2 A , that

JHEP01(2022)098
would vanish if µ = 0. We have also found that the majoron has a non-vanishing component in the doublet directions. Therefore, in order to avoid phenomenological problems with a doublet majoron, such as a sizable invisible Z-boson width, we are forced to impose the hierarchy of which guarantees that the majoron is mostly singlet. We turn now to the charged scalar mass matrix. In this case, the scalar potential contains the term with

HS +λ
(2) (2.33) Again, the application of the tadpole equations leads to One of the eigenvalues of this matrix vanishes. This corresponds to the Goldstone boson eaten by the W boson, w. The other state is the massive charged scalar C ± . In the basis {H ± , S ± }, they are given in terms of the gauge eigenstates as and their masses are

Lepton masses
The light neutrinos get their masses by means of a standard type-I seesaw mechanism. Defining the 3 × 3 matrices in generation space m D and M R the neutral leptons mass term is given by

JHEP01(2022)098
with the 6 × 6 matrix M N defined as The resulting mass matrix corresponds to the standard type-I seesaw matrix. If m D M R , the light neutrinos mass matrix is given by the well-known formula m ν = −m T D M −1 R m D . We note that the hierarchy m D M R follows naturally from the hierarchy in eq. (2.31). On the other hand, the mass term of the charged leptons reads where the 4 × 4 matrix M C is given by and we have defined The matrices m e , m ρ and m S are 3 × 3, 3 × 1 and 1 × 3, respectively. The mass matrices M N and M C can be brought to diagonal form as where V ν and V L,R are unitary matrices. 3 It proves convenient to derive approximate expressions for the matrices involved in eq. (2.45). We now follow [19] to obtain approximate expressions for the matrices V L,R by performing a perturbative expansion in powers of the inverse of the largest scale in M C . 4 We assume m e , m S m ρ M F , (2.46) consistent with eq. (2.31). Then, we can write the matrices V L,R as where U L,R and D L,R are unitary matrices. The matrices U L,R will be responsible for the block-diagonalization of M C while D L,R will diagonalize the light and heavy sub-blocks. The matrices D L,R can be written in the form

JHEP01(2022)098
where D L,R e are 3 × 3 matrices and D L,R F are just complex phases. While the non-vanishing elements of the D matrices are expected to be of O(1), the elements of U will instead be sensitive to the hierarchy of the scales involved in M C . We can decompose the 4 × 4 unitary matrices U L,R in the block form where U L,R ee , U L,R eF and U L,R F e are 3 × 3, 3 × 1 and 1 × 3 matrices, respectively, and U L,R F F are complex numbers. We impose that the following unitary transformation brings M C into a block-diagonal matrix, namely that Eq. (2.50) imposes constraints on the U L,R matrices and hence reduces their numbers of independent parameters. In particular, it requires U L,R to lead to two vanishing 3 × 1 and 1 × 3 submatrices. Therefore, each of them must have three degrees of freedom only. We then formulate the ansätze for U L and U R where I 3 is the 3 × 3 identity matrix and L and R are 3 × 1 matrices. The L and R matrices must be determined perturbatively as a function of the parameters in M C by expanding in powers of 1/M F . One must also Taylor-expand the square root. In case of L, the expansion is given by L = L 1 + L 2 + L 3 + · · · (2.53) and where the L i matrices are proportional to M −i F . Analogous expansions can be given for R. The coefficients of the expansion are computed recursively, imposing that the off-diagonal sub-blocks of M C vanish at each order in M F . Using the aforementioned hierarchy among scales, this procedure leads to

JHEP01(2022)098
where we kept only the leading order contribution in m ρ to the O(M −3 F ) coefficients L 3 and R 3 . The block-diagonal masses for the light and heavy charged leptons are finally given by In summary, m light ≈ m e and m heavy ≈ M F , with corrections to these zeroth order results entering at different orders in 1/M F .

Majoron couplings
In the gauge basis, the interaction terms of the majoron with neutrinos and charged leptons read The majoron profile in eq. (2.28) has been used in the derivation of eqs. (3.1) and (3.2). We can now focus on the interaction Lagrangian involving charged leptons and write it in the fermion mass basis. This results in where α, β are flavor indices, specified here for the sake of clarity, we have defined the four component array in flavor space X = (e, F ) and is the matrix in eq. (3.2). By comparing to eq. (1.1), one finds a dictionary between the S L,R effective couplings and the parameters of the model under consideration. In the case of the flavor violating couplings, with β = α, one finds This matching only holds for the light charged leptons, hence α, β = 1, 2, 3 here. We note that the matching is completely specified by the charged lepton mass matrix M C . In the case of the flavor conserving couplings, with β = α, one gets We note that there is a mismatch of a factor of 2 with the matching holding for the off-diagonal couplings, which prevents us from writing a single matching relation. The coupling in eq. (3.7) is purely imaginary, as expected for a pure pseudoscalar. The analytic proof of this result is given in appendix A. Finally, one can find approximate expressions for the S L,R couplings by using the expressions derived for the V L,R matrices in the previous section. One finds

Phenomenology of the model
The model can be probed thanks to its signatures in low-energy flavor experiments and high-energy colliders.
Lepton flavor violating signatures. The first and most evident consequence of a massless majoron with tree-level LFV couplings is the existence of large LFV rates in wide regions of the parameter space. This includes the usual LFV processes, such as α → β γ or α → 3 β , which constitute important constraints for our model. For the radiative processes α → β γ we use the general formulas in [20], whereas for the 3-body LFV decays γ , we use the general expressions in [21]. The effective coefficients for the 3-body decays are generated in our model at tree-level and are listed in appendix B. In addition, we must consider processes involving the majoron in the final state. In particular, the LFV decay α → β J is induced by the off-diagonal S βα A scalar couplings, with A = L, R, defined in eq. (1.1). At leading order in m β /m α , the α → β J decay width is given by [7] Searches for α → β J have been performed by several experiments. The nonobservation of these processes has been used to set stringent limits on the corresponding charged lepton LFV branching ratios. In the following we will focus on the muon decay µ → e J. In this case, the strongest limit was obtained at TRIUMF, finding BR (µ → e J) < 2.6 × 10 −6 at 90% C.L. [22]. However, the high polarization of the muon beam used in this experiment implies that the limit is only strictly valid for a purely JHEP01(2022)098 right-handed µ − e − J interaction, with S eµ L = 0, as discussed in [18]. This reference estimates a more general limit of the order of BR (µ → e J) 10 −5 . A very similar bound was recently obtained by the TWIST collaboration [23].
Anomalous magnetic moment of the muon. The enlarged lepton sector in our model induces contributions to many leptonic observables. We have already discussed flavor violating observables, which vanish in the SM. In addition, flavor conserving observables also receive new contributions, and these may potentially induce deviations from the SM predictions. For instance, the new states contribute to the anomalous magnetic moment of the muon, an observable that has received a lot of attention recently.
The anomalous magnetic moments of charged leptons, with α = e, µ, τ , are described by the effective Hamiltonian [24] where F µν is the electromagnetic field strength tensor. One can obtain the anomalous magnetic moment of the charged lepton α as A discrepancy between the SM prediction for the muon g − 2 and its experimentally determined value has existed for a long time. The interest in this deviation has increased notably after the Muon g − 2 experiment announced its first results [17]. The combination of their measurement with that obtained by the E821 experiment at Brookhaven [25] leads to a 4.2σ discrepancy with the SM prediction [16], which can be quantified as One should note, however, that a recent calculation of the hadronic vacuum polarization contribution does not favor such a large deviation [26]. We therefore need additional experimental data, soon to be provided by the Muon g − 2 collaboration, as well as more theoretical cross-checks, to firmly establish the presence of new physics in the muon g − 2.
In our model, new contributions to the muon g − 2 are induced at the 1-loop level, as shown in figure 1. This figure includes diagrams with the massive CP-even bosons H k , with the massive CP-odd state A and with the majoron J. For the massive states we use the analytical expressions given in [24], whereas the majoron contribution was recently computed in [7].
Higgs boson decays. The lightest CP-even scalar mass eigenstate, H 1 ≡ h, can be identified with the 125 GeV state discovered at the LHC. Therefore, it is crucial that its properties and decay channels match those observed, within the ranges allowed by the experimental errors. Since the observed state resembles the Higgs boson of the SM, this is JHEP01(2022)098 guaranteed if the mixing angles in the CP-even scalar sector are sufficiently small. In this case, h is made up mostly by the S H state, which has the properties of a SM Higgs. Thus, decays such as h → ZZ, W W are not modified substantially. Higgs decays into a pair a of muons have been recently searched for by the ATLAS [27] and CMS [28] collaborations. Their data yields the following ratio in terms of the SM predicted value [6], In our model, the mixing in the CP-even scalar sector may induce large deviations from the SM predicted ratio, R hµµ = 1. Even if the mixing is tiny, the large S and σ couplings to muons, required if one wants to address the experimental anomaly in the muon anomalous magnetic moment, induce sizable contributions to R hµµ , which can largely deviate from 1. An approximate analytical expression for R hµµ , valid under some simplifying assumptions, is provided in appendix C. Finally, h can also decay invisibly to a pair of majorons. The current ATLAS limit on the invisible Higgs branching ratio translates into BR(h → JJ) < 0.11 at 95% C.L. [29]. We take this constraints into account in our numerical analysis, using Γ h ≈ Γ SM h = 4.1 MeV for the total Higgs decay width [30].

Numerical results
In this section we present our numerical results. These have been obtained by randomly scanning in the wide parameter space of the model and computing the observables discussed in the previous section. All points in our scans are compatible with current neutrino oscillation data. This is achieved by using a Casas-Ibarra parametrization [31] for the Y ν Yukawa matrix, which can be expressed as Here, we are working in the basis in which the M R matrix, defined in eq. (2.38), is diagonal, In our numerical analysis we assume normal neutrino mass ordering and randomly take neutrino oscillation parameters within the 3 σ ranges obtained by the global fit [32].
Several parameters are chosen randomly in our scans. These are v S , v σ , the vector-like lepton mass M F , the trilinear coupling µ as well as the ρ and Y S Yukawa couplings, where the ρ 2 and the (Y S ) 2 upper bounds are chosen below the non-perturbative regime. In addition, the κ Yukawa matrix has been taken diagonal, with κ ii also chosen randomly. All Yukawa couplings have been assumed to be real for simplicity. The ranges for the parameters that have been chosen randomly in our scans are shown on table 2. We have chosen the VEV v S in the narrow range [0.05, 0.1] GeV. This choice is motivated by the fact that the doublet S does not couple to quarks. A sizable v S VEV would imply a reduction of v H and, as a consequence of this, an increase in the Higgs boson couplings to quarks, already constrained by LHC data. Moreover, v S is also indirectly constrained due to its impact on the diagonal majoron couplings, see eq. (3.10), and a small value is generally required. Furthermore, a small value for v S motivates a similarly small value for the µ trilinear coupling. Otherwise, the heavy CP-even scalars H k and the pseudoscalar A become very heavy and their impact on the phenomenology negligible, see eqs. (2.23) and (2.30). Finally, the scalar potential parameters λ σ = λ (1) have been fixed in all the scans. Note, however, that the exact choice for these quartic couplings is irrelevant for the observables we study. Also, the small values of v S and the quartics assure that SM precision observables are not substantially changed in our model.

JHEP01(2022)098
Our analysis has taken into account several experimental bounds. Starting with the LFV processes, we have checked that the branching ratios of the radiative processes α → β γ as well as the 3- γ satisfy the existing bounds [6]. Regarding the decays α → β J, we have imposed the restrictions on the flavor violating couplings instead of the branching ratios. To do so, we have defined the combination and imposed the bounds derived in [7], that is, On the flavor conserving side, it is well known that astrophysics imposes very stringent constraints on the majoron couplings to charged leptons. These are obtained by considering majoron-induced cooling processes in dense astrophysical media. Regarding the majoron coupling to electrons, ref. [9] finds (at 90% C.L.) The majoron coupling to muons has also been studied, but only very recently [9,10,33].
Using the supernova SN1987A, ref. [10] finds the limit 6 The impact of this bound on the phenomenology of the model will be studied in detail in the discussion that follows. We have also made sure here that the ratios where Γ SM Z → + − is the SM predicted decay width, lie within the 95% CL range, which is estimated to be 0.995 < R Zee < 1.003 and 0.993 < R Zµµ < 1.006 [6]. With respect to Higgs decays, we have considered the very recent measurement of the process h → µµ, discussed in the previous section, and we have rejected the points in our analysis outside the range compatible with eq. (4.6) at 3 σ. Finally, points with BR (h → JJ) > 0.11 [29] have been discarded as well.
Our results for the LFV processes µ → e γ and µ → e J are shown in figure 2, which shows BR(µ → e γ) as a function of BR(µ → e J). The vertical line corresponds to the bound BR(µ → e J) < 10 −5 , already discussed in the previous section, while the horizontal line is the current limit BR(µ → e γ) < 4.2 × 10 −13 , obtained by the MEG experiment [34]. Red points correspond to parameter points that respect all astrophysical bounds, namely the bounds on S ee and S µµ in eqs. (5.4) and (5.5), while the astrophysical bound on the majoron coupling to muons is violated in the blue points. Finally, the clear points are 6 Ref. [10] also gives the more stringent limit Im S µµ < 2.1 × 10 −10 , obtained with more aggressive assumptions in the simulation of the supernova SN1987A. We have explicitly checked that our conclusions would be the same if one imposes this version of the bound. excluded due to one or several of the other experimental constraints mentioned above. First, as can be seen in this figure, our model is able to attain the current experimental bounds on BR(µ → e γ) and BR(µ → e J). Moreover, one finds no difference at all between blue and red points. This implies that the astrophysical bounds on the flavor-conserving couplings S ee and S µµ have no impact on the results for the flavor-violating observables. We also observe that a correlation between BR(µ → e γ) and BR(µ → e J) exists, although these two observables depend on different combinations of parameters. However, it is easy to understand that they are not completely independent. In the limit ρ 1 = (Y S ) 1 = 0, the vector-like fermion F does not couple to electrons. In this case, the only contributions to µ − e LFV observables come from the Y ν Yukawa matrix, which has entries of the size of ∼ 10 −7 − 10 −6 and then leads to tiny LFV branching ratios. Therefore, sizable ρ 1 or (Y S ) 1 couplings are required in order to have observable LFV, and this applies both to µ → e J and µ → e γ. Regarding other LFV processes, our numerical results also show that dipole contributions dominate the amplitude of the 3-body decay µ − → e − e + e − . This leads to strong correlations with µ → e γ, which always has a much larger branching ratio.

JHEP01(2022)098
Furthermore, in the region of parameter space covered by our numerical scan, it is easy to show that BR(µ → e J) clearly correlates with the combination of parameters v σ ρ 1 ρ 2 M −2 F . From the expression of the majoron couplings to charged leptons in eqs. (3.8) and (3.9) and the charged lepton mass matrix in eq. (2.57), derived under the assumption m ρ M F , and assuming v S v H and a large ρ 2 coupling, as motivated by the explanation of the muon g − 2 anomaly, one can obtain the approximation for the off-diagonal e − µ coupling We observe that, as long as the condition m ρ < M F is satisfied, |S eµ | actually grows when the U(1) L symmetry breaking scale v σ increases. This result seems to go against the usual decoupling behavior expected when the new physics scale becomes larger. However, when v σ is increased, the mixing between the SM-like charged leptons and the vector-like lepton F increases as well, hence enhancing the µ − e − J coupling. Eventually, when v σ is pushed above M F , eq. (5.7) becomes invalid and BR(µ → e J) starts to decrease. Several ideas to improve the current limit on µ → e J have been put forward recently. As discussed in detail in [35,36], the limit can be improved by the Mu3e experiment by looking for a bump in the continuous Michel spectrum. According to this analysis, µ → e J branching ratios above 7.3 × 10 −8 can be ruled out at 90% C.L.. Alternatively, reference [9] proposes a new phase of the MEG-II experiment with a Lyso calorimeter in the forward direction, increasing in this way the sensitivity for µ → e J. Therefore, µ → e J already excludes a region of the parameter space of the model, and this region will be substantially enlarged in the future.
In what concerns τ decays, our choice of parameters suppresses all the LFV amplitudes. Since experimental limits in the τ sector are much weaker than for the muon, we do not show plots for LFV τ decays. We note, however, that one can saturate (some of) the experimental bounds also for τ 's in our model, for the appropriate choice of (large) parameters in the 3rd generation. On the other hand, in our model it is not possible to have both, τ → e γ and τ → µ γ, with large rates at the same time, without running into conflict with µ → e γ.
Our model can also induce large contributions to the muon anomalous magnetic moment and address the current experimental anomaly. This is shown in figure 3. This figure displays ∆a µ as a function of the combination of parameters (Y S ) 2 ρ 2 /M F , which enters the Feynman diagram in figure 1. The horizontal dashed line represents the experimental central value, while the green and yellow bands correspond to the 1σ and 3σ ranges, respectively. 7 As in the previous figure, the red points respect the astrophysical bounds on S ee and S µµ , while the blue points only respect the constraint on S ee . Clear points are excluded due to one or several constraints, but are shown for illustration. We have found numerically that all diagrams, with massive scalars or with the majoron in the loop, may have comparable sizes. Interestingly, some points are found within the 1σ interval, hence providing a good explanation for the experimental value of the muon anomalous magnetic moment. These points require relatively light F fermions (with masses of the order of ∼ 1 − 2 TeV) and large (order 1) ρ 2 and (Y S ) 2 Yukawa couplings. However, they violate the bound on S µµ obtained from the supernova SN1987A, since this constraint necessarily implies a low value of ρ 2 . In fact, we note that this figure displays a large concentration on points with low values of ∆a µ in a region where all the red points are found. This region is characterized by ρ 2 1, and hence the dominant contributions to the muon g − 2 do not come from the diagram in figure 1, but are mostly induced by diagrams proportional to (Y S ) 2 2 /M 2 F . These diagrams have an external chirality flip that introduces an m µ suppresing factor and then, as is generically found in a large class of models with this feature, ∆a µ can be at most ∼ 10 −10 .
Complementary information is provided by figure 4, which shows ∆a µ as a function of the vector-like mass M F for three different values of ρ 2 = (Y S ) 2 . Blue points corresponds to ρ 2 = (Y S ) 2 = √ 4π, red points to ρ 2 = (Y S ) 2 = 1 and green points to ρ 2 = (Y S ) 2 = 0.5. This figure has been obtained with a specific parameter scan in which M F ∈ [0.75 , 10] TeV, while the ranges for the other randomly chosen parameters are as in table 2. As expected, all new physics contributions decrease for large M F and strongly depend on the value of the ρ 2 and (Y S ) 2 couplings. When ρ 2 = (Y S ) 2 = 0.5, these are not large enough to address the muon g − 2 anomaly, while when ρ 2 = (Y S ) 2 = 1 this happens in a narrow region of the parameter space characterized by very light vector-like leptons, with masses 1 TeV. Only when ρ 2 = (Y S ) 2 = √ 4π, one can find an explanation for the anomaly in a wide M F range. And even in this case, they eventually become too small to account for the measured muon g − 2. However, this happens for very large vector-like masses. In fact, one finds that vector-like masses as large as 10 TeV still allow for a 3σ explanation of the muon g − 2 anomaly. Such a large mass would make the F fermions unobservable at the LHC.
We turn our attention to Higgs boson decays. As already explained, the mixing in the CP-even scalar sector can induce large deviations from the SM predicted Higgs branching ratios. In particular, a large effective coupling to muons is induced in parameter points in which the muon g − 2 anomaly is explained. This is shown in figure 5. Here we plot the Due to the large value chosen for ρ 2 , the astrophysical bound on S µµ is not respected in this plot. Imposing this constraint would imply R hµµ ≈ 1. We observe that for large v σ and M F the new physics contributions become negligibly small and one finds R hµµ = 1. However, for lower scales one finds many parameter points leading to large deviations from the SM predicted value. In particular, for M F = 1 TeV our scan reveals points with R hµµ as large ∼ 1.4 or as low as ∼ 0.3. These extreme points are of course ruled out by the existing data, but serve as example of how easily Higgs decays into muons can deviate from the SM predictions in our setup.
We finally note that the R hµµ ratio does not correlate with other observables, due to the large number of independent contributions to the h − µ − µ coupling, see appendix C. For this reason, a definite prediction cannot be made. For instance, eqs. (C.7)-(C.9) imply that for vanishing mixing in the scalar sector, c Sσ µµ = c S S µµ = 0 and c S H µµ > c SM µµ , hence predicting R hµµ > 1. However, the α and β angles never vanish and in fact one can find R hµµ < 1 as well.

Summary
In this paper we have proposed a simple model that leads to sizable majoron flavor violating couplings to charged leptons. The particle spectrum is extended with the addition of two new scalar multiplets, as well as three right-handed neutrino singlets and a vector-like lepton. The SM symmetry is also extended with a continuous lepton number global symmetry. As JHEP01(2022)098 a result of spontaneous symmetry breaking, neutrinos acquire non-zero Majorana masses via a type-I seesaw mechanism and a massless Goldstone boson appears in the spectrum, the majoron.
Thanks to large mixings between the SM charged leptons and the vector-like lepton, sizable majoron LFV couplings are generated at tree-level. Therefore, our model constitutes a simple example of a model with tree-level off-diagonal majoron couplings, not suppressed by neutrino masses. This induces plenty of signatures in experiments looking for LFV processes. In particular, we have shown that the decay µ → e J can have large rates, close to the experimental limit. In fact, it already excludes part of the parameter space of the model.
As a by-product of our construction, other interesting phenomenological possibilities emerge: (i) an explanation to the current muon g − 2 discrepancy can be provided in large parts of the parameter space of the model, easily finding points that address the anomaly even within 1σ, and (ii) sizable deviations with respect to the SM predicted Higgs decay rates can be obtained, most notably in h → µµ. These two phenomenological possibilities provide additional handles on the model. We note, however, that an explanation of the muon g − 2 anomaly would lead to tension with recent astrophysical bounds on the majoron coupling to muons.
In this work we have shown that as soon as the lepton sector is extended beyond the minimal models, exotic signatures appear, such as those including a massless majoron in the final state. This motivates the experimental search for processes like µ → e J and estimulates the construction of new theoretical constructions that, in addition to neutrino masses, provide an understanding to other open questions.

JHEP01(2022)098
be a generic complex (3 + n) × (3 + n) complex matrix, given by the blocks m 1 , m 2 , m 3 and m 4 , with dimensions 3 × 3, 3 × n, n × 3 and n × n, respectively. The singular value decomposition of the matrix M is V † R M V L = M , where V R and V L are unitary matrices and M = diag(M 1 , M 2 ), with M 1 and M 2 3 × 3 and n × n real diagonal matrices, respectively, with positive entries. The interaction matrix of the majoron with the charged leptons can be written in the flavor basis as where N is another (3 + n) × (3 + n) matrix and x, y ∈ R. Comparing with eq. (3.4), the majoron coupling matrix in our model is given by N ≡ A and corresponds to m 1 ≡ m e , First of all, we block-parametrize the unitary matrices V R and V L as We also denote where I n is the n × n identity matrix. Therefore, one obtains Analogously, one finds the following relations involving V L : After these preliminaries, we note that the interaction matrix N can be written as Therefore, if we prove that the combinations have real diagonal elements, then also V † R N V L has real diagonal elements and the proof is complete. Let us consider the first term. It is easy to check that

JHEP01(2022)098
Then one can obtain (A.10) The diagonal terms of the resulting matrix are real because AA † and CC † are Hermitian matrices and M 1 , M 2 are real diagonal matrices. We now have to consider the second term in eq. (A.7). It is possible to write Using similar manipulations as for the first term one finds (A.12) and, writing for the sake of brevity only the diagonal blocks of this expression, in the form of a column array, we obtain (A.13) The second terms in the sum cancel for both the upper and lower diagonal blocks, using the unitarity of V L and V R . Then, using again unitarity in the following way 14) we finally end up with and therefore the diagonal components of this matrix are purely real. This concludes the proof.

B Effective coefficients for flavor violating observables
In order to use the analytical results for the flavor violating observables provided in [7,21], one must match the effective Lagrangian in these references to the specific model discussed here. We focus in particular on the 3-body lepton decays − γ . These processes get tree-level contributions in our model from the three CP-even scalars H k (k = 1, 2, 3), the Z-boson, the CP-odd scalar A and the majoron J. The majoron contributions have been computed in [7]. Since all the other mediators are significantly heavier than the SM charged leptons, we can then parametrize their contributions by the effective Lagrangian where we have defined Γ S = 1, Γ V = γ µ and Γ T = σ µν and omitted flavor indices in the effective coefficients for the sake of simplicity. The coefficients A I XY have dimensions of mass −2 . We will now give specific expressions for these coefficients in the model under consideration. In order to do that it proves convenient to define the following 3-component array

H k contributions.
Denoting the matrix in the previous equation as B and transforming the Lagrangian to the mass basis, one can easily perform the matching with eq. (B.1). We get the following expressions for the contributions of A to the effective coefficients A I XY = (A I XY ) βββα .
There are two types of Feynman diagrams contributing to this process. The first class involves a flavor conserving (γγ) and a flavor violating (βα) vertex, while in the second class both vertices violate flavor (βγ and γα). Therefore, the matching with eq. (B.1) would yield non vanishing contributions to both coefficients (A I XY ) γγβα and (A I XY ) βγγα . One can actually Fierz transform the latter flavor structure into the former, thus in the following expressions we set A I XY = (A I XY ) γγβα . The Fierz transformations involved in the matching are the following, where the type of parenthesis indicates the fermion field which is contracted with the gamma matrix in brackets.
In this process both vertices are necessarily flavor violating (γβ and γα). This allows us to easily perform the matching with eq. (B.1) and set in the following expressions A I XY = (A I XY ) γβγα .

C R hµµ analytical expression
In our model, the R hµµ ratio can be approximately written as is the SM Higgs coupling to a pair of muons and c S H µµ , c Sσ µµ and c S S µµ denote the contributions from the gauge eigenstates S H , S σ and S S , respectively. These couplings are given by where we have assumed ρ 1 , ρ 3 ρ 2 and Y S 1 , Y S 3 Y S 2 , as motivated by the explanation of the muon g−2 anomaly and the stringent constraints from lepton flavor violating observables. Furthermore, we have introduced the mixing angles α, β and γ. The CP-even scalar mass matrix M 2 R in eq. (2.23) is diagonalized by the unitary matrix R which, assuming small mixing angles, can be parametrized as where α, β, γ 1. Using now eqs. (2.55)-(2.57) we finally obtain

JHEP01(2022)098
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.