The Majoron at two loops

We present singlet-Majoron couplings to Standard Model particles through two loops at leading order in the seesaw expansion, including couplings to gauge bosons as well as flavor-changing quark interactions. We discuss and compare the relevant phenomenological constraints on Majoron production as well as decaying Majoron dark matter. A comparison with standard seesaw observables in low-scale settings highlights the importance of searches for lepton-flavor-violating two-body decays $\ell \to \ell' +$Majoron in both the muon and tau sectors.

The Standard Model (SM) has emerged as an incredibly accurate description of our world at the particle level. Even its apparently accidental symmetries, baryon number B and lepton number L, are seemingly of high quality and have never been observed to be violated. One could, however, argue that the established observation of nonzero neutrino masses is not only a sign for physics beyond the SM but also for possible lepton number violation by two units. This argument is based on an interpretation of the SM as an effective field theory (EFT) and the observation that the leading non-renormalizable operator * Julian.Heeck@uci.edu † hpatel6@ucsc.edu is Weinberg's dimension-five operator (LH) 2 /Λ [1]. This operator violates lepton number (∆L = 2) and leads to Majorana neutrino masses of order H 2 /Λ after electroweak symmetry breaking, which gives the correct neutrino mass scale for a cutoff Λ ∼ 10 14 GeV. Besides explaining neutrino oscillation data, an EFT scale this high has little impact on other observables and thus nicely accommodates the absence of non-SM-like signals in our experiments. An ever-increasing number of renormalizable realizations of the Weinberg operator exist in the literature, the simplest of which arguably being the famous type-I seesaw mechanism [2] that introduces three heavy righthanded neutrinos to the SM field content. While the Weinberg operator explicitly breaks lepton number, underlying renormalizable models could have a dynamical origin for ∆L = 2 via spontaneous breaking of the global U (1) L symmetry. This leads to the same Weinberg operator and thus Majorana neutrino masses, but as a result of the spontaneous breaking of a continuous global symmetry also a Goldstone boson appears in the spectrum. This pseudo-scalar Goldstone boson of the lepton number symmetry was proposed a long time ago and was dubbed the Majoron [3,4].
The Majoron is obviously intimately connected and coupled to Majorana neutrinos, but at loop level also receives couplings to the other SM particles. This makes it a simple renormalizable example of an axion-like particle (ALP), defined essentially as a light pseudo-scalar with an approximate shift symmetry. Although not our focus here, by coupling the Majoron to quarks it is even possible to identify it with the QCD axion [5][6][7][8][9], thus solving the strong CP problem dynamically. The main appeal of the Majoron ALP is that its couplings are not free but rather specified by the seesaw parameters, which opens up the possibility to reconstruct the seesaw Lagrangian by measuring the Majoron couplings [10]. This is aided by the fact that the loop-induced effective operators that couple the Majoron to the SM are only suppressed by one power of the lepton-number breaking scale Λ, whereas right-handed neutrino-induced operators without Majorons are necessarily suppressed by 1/Λ 2 [11][12][13][14][15][16][17], rendering it difficult to reconstruct the seesaw parameters in that way.
In this article we complete the program that was started in the inaugural Majoron article [3] and derive all Majoron couplings to SM particles. The tree-level and one-loop couplings were obtained a long time ago; here we go to the two-loop level in order to calculate the remaining couplings, which include the phenomenologically important couplings to photons as well as to quarks of different generations. Armed with this complete set of couplings we then discuss various phenomenological consequences and constraints on the parameters. This includes a discussion of Majorons as dark matter (DM).
The rest of this article is structured as follows: in Sec. II we introduce the singlet Majoron model and reproduce the known tree-level and one-loop couplings. In Sec. III we present results of our novel two-loop calculations necessary for the Majoron couplings to gauge bosons and to quarks of different generations. The phenomenological aspects of all these couplings are discussed in Sec. IV. Finally, we conclude in Sec. V.

II. MAJORON COUPLINGS AT TREE LEVEL AND ONE LOOP
In this article we consider the minimal singlet Majoron model [3], which introduces three right-handed neutrinos N R coupled to the SM lepton doublets L and Higgs doublet H, and one SM singlet complex scalar σ carrying lepton number L = −2, minimally coupled to the righthanded neutrinos proportional to the Yukawa matrices y and λ. We do not specify the scalar potential V(H, σ) but simply assume that σ = (f + σ 0 + iJ)/ √ 2 obtains a vacuum expectation value f , which then gives rise to the right-handed Majorana mass matrix M R = f λ/ √ 2. J is the Majoron, σ 0 is a massive CP-even scalar with mass around f , assumed to be inaccessibly heavy in the following. Both M R and the charged-lepton mass matrix are chosen to be diagonal without loss of generality, effectively shifting all mixing parameters into y. Electroweak symmetry breaking via H = (v/ √ 2, 0) T yields the Dirac mass matrix M D = yv/ √ 2. The full 6 × 6 neutrino mass matrix in the basis (ν c L , N R ) = V n R is then where V is the unitary 6 × 6 mixing matrix to the states n R , which form the Majorana mass eigenstates n = n R + n c R . The diagonal mass matrix M n = diag(m 1 , . . . , m 6 ) consists of the physical neutrino masses arranged in ascending order. Throughout this article, we denote mass matrices with capital letters M x and individual mass eigenvalues with small letters m i . In the mass eigenstate basis, the tree-level neutrino couplings to J, Z, W − , and h, take the form [18,19] where g w = e/ sin θ W with Weinberg angle θ W and In the last equation we used V αk = 1 1 αk since we work in the basis where the charged-lepton mass matrix is diagonal. The 6 × 6 matrix C and the 3 × 6 matrix B satisfy a number of identities [20,21] that are particularly important in order to establish ultraviolet (UV) finiteness of amplitudes involving neutrino loops: So far we have not made any assumption about the scale of M R . In the following we will work in the seesaw limit M D M R , resulting in a split neutrino spectrum with three heavy neutrinos with mass matrix M R and three light neutrinos with seesaw mass matrix naturally suppressed compared to the electroweak scale v. This hierarchy permits an expansion of all relevant matrices in terms of the small where U is the unitary 3 × 3 Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix. Parametrically this corresponds to an expansion in the scale hierarchy v f which we refer to as the seesaw expansion. To leading order, the matrices take the form Note that AM R A T is diagonal, which imposes constraints on A and provides an implicit definition of U . These constraints may be automatically satisfied using the Casas-Ibarra parametrization [22]; however, more useful for our purpose is the Davidson-Ibarra parametrization [23], which uses M ν = −AM R A T and M D M † D as the independent matrices containing all seesaw parameters. Since M ν is essentially already fixed by neutrino oscillation experiments (modulo the phases, hierarchy, and overall mass scale), the next step is to experimentally determine M D M † D . As we will see, this could in principle be achieved by measuring Majoron couplings without ever observing the heavy right-handed neutrinos.
To this effect let us point out some interesting properties of the hermitian matrix M D M † D [10]: its determinant is simply det M D M † D = det M n = 6 j=1 m j , which is strictly positive in the model at hand even if one of the light neutrinos were massless at tree level [24]. Thus, M D M † D is positive definite, which yields a chain of inequalities for the off-diagonal entries (M D M † D ) ij , i = j (see e.g. Ref. [25]): This provides a useful way to constrain magnitudes of the elements of M D M † D since its trace appears in many couplings of the Majoron.
From Eq. (3) all loop-induced Majoron couplings are necessarily proportional to 1/f . But many couplings contain additional powers of M −1 R ∝ 1/f , which make them higher order in the seesaw expansion. We will neglect these suppressed couplings and focus on those that are down by only one power of 1/f . For the sake of generality, we determine the Majoron couplings assuming an explicit shift-symmetry-breaking Majoron mass term − 1 2 m 2 J J 2 , making J a pseudo-Goldstone boson. This mass could be explicit [26,27] or arise from quantumgravity effects [28][29][30].

A. Neutrino couplings
By inserting Eq. (6) into Eq. (3), the tree-level Majoron coupling to the light active Majorana neutrinos in the seesaw limit is These diagonal Majoron couplings to neutrinos are formally second order in the seesaw expansion since The omitted offdiagonal Jn i n j couplings are determined by the ma- 3 , which are further suppressed, and lead to irrelevantly slow active-neutrino decays n i → n j J [4].
Assuming for simplicity m J m 1,2,3 , the Majoron's partial decay rate into light neutrinos is For sufficiently large f the Majoron becomes a long-lived DM candidate [27,29,[31][32][33][34][35][36], discussed in Sec. IV C. As mentioned earlier, the Majoron couplings to all other SM particles are leading order in the seesaw expansion, i.e. proportional to 1/f , and may easily dominate the phenomenology despite the additional loop suppression [10]. Therefore, a thorough discussion of the Majoron requires knowledge of all loop-induced couplings that are leading order in the seesaw expansion. Using the tree-level couplings of Eq. (3) we calculate the loopinduced Majoron couplings to the rest of the SM particles and provide them below.

B. Charged fermion couplings
The leading order couplings to charged fermions are obtained from the one-loop diagrams in Fig. 1. These were calculated long ago, both in the one-generation case [3] and in the three-generation case, which leads to off-diagonal Majoron couplings to leptons [18]. At leading order in the seesaw expansion, these couplings take a simple form [10], with (diagonal) quark couplings and charged lepton couplings where M ,u,d denote the diagonal mass matrices of the appropriate SM fermions. In addition to exhibiting decoupling in the seesaw limit M R ∼ f → ∞, these couplings vanish in the electroweak symmetric limit v → 0 as expected since J is an electroweak singlet. The quark couplings can be used to derive the Majoron couplings to nucleons N = (p, n) T , using the values from Ref. [37]: At this point let us make some remarks about CP violation. Already in the one-loop processes above one encounters loop-induced Majoron mixing with the Brout-Englert-Higgs boson h, which would result in Majoron couplings to the scalar bilinearf f as opposed to the pseudo-scalarf iγ 5 f . It was noted in Ref. [18] that the relevant J-h mixing diagrams vanish for m J = 0. For m J = 0 the J-h amplitude is of order (v/f ) 2 in seesaw and hence negligible. This can be understood by noting that CP-violating phases in the Davidson-Ibarra parametrization reside both in the active-neutrino mass

C. Couplings to gauge bosons
At one-loop order, the only non-vanishing couplings to gauge bosons are to W W and ZZ, with typical diagrams shown in Fig. 2. However, they are higher order in the seesaw expansion, which can be understood as follows: the amplitudes come with a factor of M D M † D /f in order to achieve the necessary N R -ν L mixing to close the loop; on dimensional grounds there is an additional M −2 R suppression since this is the only high mass scale in the loop. Explicit one-loop formulae can be found in Refs. [38,39]. The leading seesaw behavior for coupling to gauge bosons without the M −2 R suppression starts at two-loop order. At this point it is appropriate to discuss the connection to anomalies in the minimal Majoron model. There is some confusion in the literature regarding the question of whether the Majoron is the Goldstone boson of the anomaly-free U (1) B−L or the anomalous U (1) L . Both choices seem equally valid because baryon number remains unaltered by the Lagrangian in Eq. (1). Since according to common lore Goldstone couplings to gauge bosons are determined by the anomaly structure of the theory, this leads to a paradox when attempting to guess the form of Majoron couplings to W and Z. The resolution was recently presented in Ref. [40], where it was explained that Goldstone couplings to gauge bosons are driven entirely by non-anomalous processes. Anomalies still serve as a useful bookkeeping device for the couplings to vector-like gauge bosons such as gluons and photons, but fail for chiral gauge bosons. Disregarding anomalies it is then necessary to calculate Goldstone couplings to gauge bosons in perturbation theory, the results of which we present in the next section. Additionally, we emphasize that the non-vanishing Majoron coupling to electroweak gauge bosons (Eqs. (24), (28), and (35) below) does not lead to nonperturbative violation of the shift symmetry beyond m J due to the absence of electroweak instantons without B + L violation [41][42][43]. Therefore, electroweak instantons cannot generate Majoron mass.

III. MAJORON COUPLINGS AT TWO LOOPS
In this section we present Majoron couplings for which the leading seesaw behavior arises at two loops. To automate the evaluation of some O(100) Feynman diagrams contributing to each effective coupling we implemented this model in Feynman gauge including all Goldstone bosons [19] in FeynRules [44], and generated the necessary amplitudes with FeynArts [45]. We validated our implementation by reproducing the tree-level and oneloop couplings above.
The Feynman diagrams naturally divide into two sets (see Fig. 3). Set I diagrams contain the one-particle irreducible (1PI) two-loop diagrams. Set II diagrams contain the reducible diagrams that are dominated by J-Z mixing. We used an in-house Mathematica implementation of expansion by regions as described in Ref. [46] to carry out a double asymptotic expansion m 4,5,6 → ∞, m 1,2,3 → 0 of the two-loop vertex integrals, and to algebraically reduce the one-loop [47] and two-loop [48] tensor integrals in dimensional regularization. We treated γ 5 naively, such that it anti-commutes with all other Dirac matrices while also preserving the cyclic property of traces. Finally, after expanding around four spacetime dimensions, we summed over fermion generations to extract the leading seesaw behavior of each Majoron coupling. We found the couplings to be expressible as simple sums of one-loop functions and rational terms. In principle, two-loop self-energy and vacuum integrals may be present, but they cancel away in the course of reduction, leaving behind rational terms.
We have checked our results by confirming that all amplitudes are proportional to the expected tensor structures, are UV finite upon using the relations of Eq. (5), and have the expected limiting low-energy/small-mass behavior. Additionally, we confirmed that our results are insensitive to the treatment of γ 5 by reevaluating them in several different ways, including projecting the integrals onto form factors and also starting from cyclically reordered Dirac traces, and finding the same answer upon expanding around four spacetime dimensions.
We pause to comment on how we quote our results for couplings to gauge bosons {V V } = {gg, γγ, Zγ, ZZ, W + W − }. We phrase our results in terms of on-shell decay amplitudes It is commonplace to see these amplitudes interpreted as effective couplings as they appear to match onto EFT operators of the form [49] where V µν is the appropriate field-strength tensor and V µν its dual. However, we caution the reader that the identification with effective couplings in this way is somewhat clumsy for the following reasons. First, matching onto local operators should be carried out for off-shell Green functions which have been expanded in the external momenta. Second, our effective couplings g JV V cannot be viewed in a Wilsonian sense, since degrees of freedom lighter than the Majoron contribute in certain mass ranges, nor can it be viewed in the 1PI sense since the couplings include Set II diagrams that are not oneparticle irreducible. Therefore, interpreting our results as coefficients of effective operators should be done with care.

A. Coupling to two gluons
Assuming a sufficiently heavy Majoron m J > ∼ Λ QCD , the coupling to free gluons comes entirely from J-Z mixing diagrams of Set II in Fig. 3. A straightforward evaluation of the decay amplitude J → gg at leading order in the seesaw expansion yields the simple expression with T u 3 = −T d 3 = 1/2 and the loop function For small m J , the amplitude vanishes as g Jgg ∼ m 2 J and indicates that at leading order in the derivative expansion the amplitude matches onto an operator (∂ 2 J)G aµνGa µν instead of JG aµνGa µν as in Eq. (14). This implies that the Majoron does not solve the strong CP problem, as this operator is insensitive to a constant shift J → J + c that could otherwise be used to cancel the strong CP θ term [40]. Furthermore, contrary to the claim in Ref. [38], Majorons without tree-level couplings to quarks cannot solve the strong CP problem even at higher loop order [50,51].

B. Coupling to two photons
The Majoron coupling to photons at two-loop level receives contributions from both sets of Feynman diagrams in Fig. 3, and yields the partial decay rate into two photons, The contributions from Set II were calculated in Ref. [10] with the result already simplified with the help of the electroweak anomaly cancellation condition Just as for the gluon coupling, the amplitude vanishes as g Jγγ ∼ m 2 J for small Majoron masses, implying that the leading effective operator this amplitude matches onto in the derivative expansion is (∂ 2 J)F µνF µν rather than the typically occurring JF µνF µν . For m J m f we can relate the total Majoron-photon coupling g Jγγ to the dimensionless diagonal fermion couplings of Eqs. (10) and (11), g Jf f Jf iγ 5 f , as which agrees with the EFT result of Ref. [52]. Since the g Jf f couplings can have different signs and magnitudes, the Jγγ coupling for m J < m e could be heavily suppressed. The key point and crucial result of this full two-loop calculation is that the Jγγ coupling has a richer structure than anticipated in Ref. [10] based on the evaluation of g II Jγγ alone. This is illustrated in Fig. 4 where we show |g Jγγ | × f for a variety of hierarchies of the diagonal entries (M D M † D ) . The SMfermion mass thresholds together with the different signs in g Jf f potentially suppress g Jγγ by orders of magnitude. The typical size of the coupling for m J > MeV is |g Jγγ | ∼ 10 −5 f −1 (M D M † D )/(100 GeV) 2 , simply due to the unavoidable suppression factor α/(8π 3 ). In Fig. 4 we used the current-quark masses to evaluate g Jγγ ; for m J < ∼ Λ QCD they should be replaced by hadronic loops. We have not attempted this, but we refer the interested reader to standard axion literature on the topic [53][54][55].

C. Coupling to Z and photon
Next we present the Z-photon coupling g JZγ = g I JZγ + g II JZγ , which receives contributions from Set I and Set II diagrams in Fig. 3. The results are with c W ≡ cos θ W and s W ≡ sin θ W . We have used to simplify the formula. In the limit m J , m Z → 0, the amplitude is non-vanishing, and matches onto the effective operator JZ µνF µν , as in Eq. (14).

D. Coupling to two Z bosons
The Majoron coupling g JZZ = g I JZZ + g II JZZ to two Z bosons receives contributions from Sets I and II: where and C 0 is the scalar three-point Passarino-Veltman function. Despite the appearance of (m 2 J −4m 2 Z ) −1 , the amplitude is regular at threshold m J → 2m Z . The coupling g JZZ is nonvanishing in the limit m J , m Z → 0, and matches onto JZ µνZ µν .

E. Coupling to two W bosons
Finally, we present results for the two-loop amplitude for J → W + W − in order to extract the coupling g JW W , which receives contributions from diagrams in Set I and Set II where we have separated the Set II J-Z mixing contributions based on the type of SM fermions running in the loop. The Set I diagrams give the Set II J-Z mixing diagrams with two charged leptons in the loop give the Set II J-Z mixing diagrams with two down quarks in the loop give and the Set II J-Z mixing diagrams with two up quarks in the loop give Here, V q is the unitary Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. In the last two formulae, the twoargument loop function is The coupling is regular at threshold m J → 2m W , and nonvanishing in the limit m J , m W → 0, and matches onto JW + µν W − µν .

F. Coupling to γ-Higgs and Z-Higgs
Besides the usually considered pseudo-scalar couplings to two gauge bosons discussed above, CP-invariance also allows couplings of J to hZ and hγ. The former arises already at one-loop level but is seesaw suppressed; the dominant contributions to both couplings then arise at two-loop level. Because of the large number of diagrams and low phenomenological relevance of decays such as h → γJ compared to the processes derived above we will not, however, present the results here.

G. Flavor changing quark couplings
At the two-loop level we find off-diagonal Majoron couplings to quarks, which can lead, for example, to s → dJ or K → πJ at hadron level. Such flavor-changing couplings have long been advocated to search for light bosons and axions [56,57] and have enjoyed increased attention in recent years [58][59][60][61], partly because of an improved reach at existing and upcoming experiments such as NA62 and Belle II.
In the Majoron model the relevant flavor-changing quark-level couplings arise at two loops and involve a large number of diagrams; see Fig. 5. The leading logarithmic contribution to coupling to down quarks is arising from the Set I diagrams in Fig. 5. The couplings are of minimal-flavor-violating (MFV) type [62], as expected from the fact that J is quark-flavor blind. The subleading contribution from the remaining diagrams is given by where the diagonal matrix F u has entries For sub-TeV Majoron mass and large-enough righthanded neutrino masses, log (M R /m W ) > ∼ 1 as is assumed in our seesaw expansion, the leading logarithmic contribution L Jdd dominates over L sub Jdd , which we neglect in the following for simplicity.
The MFV matrix M d V † q M 2 u V q relevant for L Jdd makes it clear that off-diagonal terms would vanish if all up quarks were degenerate, so these terms are necessarily proportional to up-quark mass differences. Numerically this matrix evaluates to (39) keeping only the largest entries. The biggest amplitude is therefore b → sJ, which is, however, experimentally less clean than s → dJ, further discussed in Sec. IV A.
From the dominant flavor-changing down-quark couplings in L Jdd we immediately obtain the corresponding flavor-changing up-quark couplings as with the markedly smaller MFV coupling matrix Taken together with the weaker experimental limits on u → u J we can ignore these couplings in practice.
From the point of view of the seesaw expansion, the quark-flavor changing couplings are actually the dominant Majoron couplings. They only decouple as log(M R )/f , whereas all other couplings decouple at least as 1/f . It was noted before that Goldstone bosons with effective diagonal couplings Jm qq iγ 5 q yield flavorchanging quark couplings at one loop that depend logarithmically on the UV scale [63,64], whereas an initial coupling JW µνW µν does not have such a dependence [58]. In our case JW µνW µν gives only a seesawsuppressed contribution to Jqq , and the log(M R ) terms originate from an effective coupling to Goldstone bosons, JG + G − .

IV. PHENOMENOLOGY
Having obtained all Majoron couplings to leading order in the seesaw expansion we can discuss existing constraints and signatures.

A. Light Majorons
We start with the simplest case of a massless Majoron, which most importantly gives a vanishing coupling to photons. It proves convenient to phrase our discussion in terms of the dimensionless parameters as they capture the Majoron couplings in most cases [10]. The off-diagonal entries of M D M † D are directly constrained by the lepton-flavor-violating (LFV) decays → J [18,65,66]. For m m , the partial widths read and involve a left-handed final-state lepton, leading to an anisotropic decay [10]. The constraints in the tau sector are Br(τ → J) < O(10 −3 ) [67] and lead to [10] |K τ e | < 6 × 10 −3 , which can be improved by Belle and Belle-II [68][69][70]. In the muon sector, the best constraints on a Majoron with anisotropic emission come from µ → eJ [71] (to be improved with Mu3e [72]) and µ → eJγ [73]. The latter is also sensitive to m J = 0 and provides a limit The diagonal couplings K of a massless Majoron are constrained by astrophysics, stellar cooling in particular, and imply from the electron [74] and nucleon coupling [75], respectively. The bound on the trace tr(K) is particularly powerful since it provides upper bounds on all entries of the positive-definite K [10] by means of the inequality of Eq. (7). This then puts an upper bound on the Majoron coupling to muons and taus, which is far better than any direct bound on these couplings. It also ensures that rare decays such as K → πJ and Z → γJ [54,55,76,77] are unobservably suppressed for a massless Majoron. Two-loop couplings are hence irrelevant for massless Majorons.
Overall we see that a massless Majoron gives seesawparameter constraints of order While this is far off the "natural" value , it can be realized by assuming certain matrix structures in M D that suppress to be discussed in more detail in Sec. IV B. As we have seen, the relevant couplings of a massless Majoron are those to nucleons and electrons, but even µ → eJ could be observable.
The phenomenology becomes more interesting for nonzero Majoron mass, specifically values above ∼ 10 keV in order to kinematically evade the stellar cooling constraints. This is shown in Fig. 6 for Majoron masses above 1 MeV. In addition to the one-loop lepton-flavorviolating decays that probe K αβ we now also have relevant constraints from the two-loop quark-flavor-violating Upper limits on combinations of K αβ = (MDM † D ) αβ /(vf ) for Majoron masses above MeV. The shaded regions exclude |K αβ | or tr(K) by non-observed rare decays, the dashed lines show the potential future reach; see text for details. K → πJ and B → KJ further scale with log(MR/mW ), which has been set to 1 1 here. The black region is a very naive estimate of SN1987 constraints on the di-photon coupling, setting for simplicity K = tr(K)/3. The yellow region is the SN1987 constraint on the JN N coupling [78]. The off-diagonal entries have to satisfy |K αβ | < tr(K)/2; see Eq. (7). decays, especially K → πJ and B → KJ (limits and future sensitivity taken from Ref. [61], based on Refs. [58,79,80]). These decays probe the quantity tr(M D log(M R /m W )M † D ), but to simplify comparison with the other limits we set log(M R /m W ) = 1 1 to obtain a limit on tr(K). It should be kept in mind, however, that a larger log enhancement can make these rare processes even more relevant. Also potentially relevant are the twoloop Majoron couplings to photons via the effective coupling g Jγγ from Sec. III B. Astrophysical limits on this coupling are extremely strong for m J < ∼ 100 MeV [81]; since g Jγγ is m 2 J suppressed for m J < ∼ MeV, the region where the photon coupling is important is between MeV and 100 MeV. This is unfortunately precisely the mass region where the light quarks that run in the J-γ-γ loops should be replaced by hadrons; as a very naive way to incorporate this we simply set m u = m d = m π and m s = m K in Eq. (19). The g Jγγ limit in Fig. 6 should therefore not be taken too seriously. In light of these uncertainties we do not discuss how the g Jγγ coupling depends on the various diagonal K , but rather set them all equal to tr(K)/3 to allow a comparison to the other limits. It is clearly possible to suppress g Jγγ significantly in the region of interest by choosing hierarchical K , as shown in Fig. 4.
As can be appreciated from Fig. 6, even the strong astrophysical constraints on Majorons do not rule out flavor-violating rare decays, with significant experimen-tal progress expected in the near future. Even the twoloop suppressed d → d J decays provide meaningful constraints.
µ → eJ is well constrained already, and we expect it to eventually become the most sensitive probe of the K matrix entries for m J < m µ , even beating out stellar cooling limits. τ → J, on the other hand, is mainly relevant for m J between ∼ 100 MeV and m τ . For smaller m J the flavor-conserving constraints on K from g Jγγ and g JN N become stronger, which suppresses the LFV modes via the inequality of Eq. (7). For m J > m τ , the main rare decays are B → KJ and Z → Jγ [49], the former is shown in Fig. 6, and the latter gives irrelevant constraints on the K entries of order 10 4 .

B. Comparison with seesaw observables
So far we have discussed the interactions of the Majoron, but, of course, the right-handed neutrinos N R also mediate non-Majoron processes, discussed at length in the literature [11][12][13][14][15][16][17]. Assuming again that the N R are heavy enough to be integrated out, the relevant dimension-six operators involving SM fields all depend on the matrices which drive LFV processes such as → γ as well as lepton-universality violating effects such as Γ(Z → ¯ )/Γ(Z → ¯ ), recently discussed thoroughly in Ref. [17]. In comparison, we have seen above that all Majoron operators depend on the matrices For f ∼ M R , this makes the Majoron operators potentially dominant, while f M R suppresses them to an arbitrary degree [10]. To properly compare Majoron and non-Majoron processes it is necessary to pick a structure for M D , which is guided by our experimental reach.
As we have seen above, even future limits on Majoron production cannot reach the natural seesaw scale M R ∼ 10 14 GeV, and the same is true for other N Rmediated processes [17]. By no means does this preclude observable effects, since it is possible to use the ma- [82,83], potentially realizing a lepton-number symmetry [84] as in the inverse seesaw [85][86][87]. Following Refs. [17,88] we can solve M ν = 0 for M D , which then requires only tiny perturbations δM D to produce the observable neutrino masses via The key observation is that δM D is negligible in M D M † D . M ν = 0 requires the low-scale seesaw structure with complex z and real ξ α without loss of generality [17]. 1 This product structure of M D implies for any function b. The off-diagonal entries are then real and entirely determined by the diagonal ones: Notice that this violates the strict inequality derived earlier in Eq. (7), which assumed non-vanishing neutrino masses. Equation (53) drastically simplifies the discussion of the seesaw parameter space, seeing as all observables now only depend on the three real diagonal entries of M D b(M R )M † D instead of the nine parameters it could contain in general. Furthermore, all Majoron and non-Majoron parameter matrices have the same flavor structure and only differ in their absolute magnitude. Equation (53) also shows that any low-scale seesaw texture automatically maximizes the off-diagonal flavor-violating entries of the relevant coupling matrices (Eq. (48) or Eq. (49)).
We compare Majoron and non-Majoron limits in Fig. 7, the latter adopted from Ref. [17]. We set M R /TeV = 1 1 as well as f = 1 TeV and stress again that the Majoron limits can be suppressed arbitrarily by increasing f .  [17], the non-Majoron limits are completely dominated by flavor-conserving observables such as Z → e + e − and other electroweak precision data. LFV such as τ → eγ and τ → eee reside far in the already excluded region and even future improvements, e.g. in Belle-II, do not reach the allowed parameter space. In comparison, Majorons with masses between ∼ 100 MeV and m τ do give relevant constraints from τ → eJ and will probe significantly more parameter space with upcoming Belle and Belle-II analyses [70]. Standard LFV in the τ e (and τ µ in complete analogy) sector are hence doomed to be unobservable in the seesaw model, but the Majoron LFV channels τ → J could be observable and deserve more experimental attention.  [90], provides important constraints on the parameter space, and all future µe LFV will probe uncharted terrain. For m J < m µ , the Majoron channel µ → eJ already sets better limits than µN → eN and can continue to dominate over µ → eγ and µ → 3e in the future. Ultimately, µN → eN conversion in Mu2e [91] and COMET [92] has the best future reach.
In this region of parameter space Majorons can form DM [27,29,[31][32][33][34][35][36], with a production mechanism that can be unrelated to the small decay couplings [10,27,93,94]. The defining signature of Majoron DM is a flux of neutrinos from DM decay with E ν m J /2 and a known flavor composition [10]. For m J > ∼ MeV these neutrino lines could potentially be observable via chargedcurrent processes in detectors such as Borexino or Super-Kamiokande [10], while lower masses are more difficult to probe [95].
The loop-level couplings generate the much more constrained decays J → ff , γγ, which, however, depend on the matrix K and are hence complementary to the neutrino signature, as discussed in detail in Ref. [10]. This analysis used only the Set II of diagrams to calculate J → γγ, namely the expression proportional to tr(K), Eq. (19). The new full expression presented in Sec. III B leads in general to a suppression of the diphoton rate and a more involved dependence on the K matrix entries (Fig. 4). We omit a full recasting of existing DM→ γγ limits onto our J → γγ expression since it is not very illuminating, but stress that this diphoton suppression makes the neutrino modes even more dominant.

V. CONCLUSION
Majorons, the Goldstone bosons of spontaneously broken lepton number, were proposed in the early 1980s in models for Majorana neutrino masses. Since then experiments have indeed found evidence for non-zero neutrino masses, although it is not yet clear whether they are of Majorana type. With the motivation for Majorons as strong as ever, we have set out in this article to complete the program that was started almost 40 years ago and calculate all Majoron couplings to SM particles. The couplings to neutrinos (tree level) as well as charged leptons and diagonal quarks (one loop) were known previously. Here we presented the two-loop couplings to gauge bosons (Jγγ, JγZ, JZZ, JW W , Jgg) and flavorchanging quarks (Jdd , Juu ). Phenomenologically relevant of these are currently only the Majoron coupling to photons as well as the Jdd couplings behind the rare decays K → πJ and B → KJ.
Standard seesaw effects in an EFT approach are encoded in the matrix M D M −2 R M † D , which drives, for example, → γ. Majoron couplings, on the other hand, depend on the matrix M D M † D /f , which is parametrically larger in the seesaw limit and can indeed give better constraints in parts of the parameter space. For example, while τ → γ and other τ LFV are unlikely to be observable in the seesaw model, τ → J can be observably large and deserves more experimental attention. The singlet Majoron model together with the coupling texture of Eq. (51) implied by low-scale seesaw is a very minimal UV-complete realization of an axion-like particle and thus a well-defined benchmark model. The dominant theoretical challenge not addressed here is the replacement of quarks by hadrons in loops, which we leave for future work. We expect future studies to elucidate additional aspects of this model, in particular when the Majoron is used as a portal to dark matter.