Probing New Physics with Long-Range Neutrino Interactions: An Effective Field Theory Approach

We investigate forces induced by the exchange of two light neutrinos between Standard Model (SM) fermions in the presence of effective operators parametrising physics beyond the SM. We first set up a general framework in which we derive the long-range potential mediated by weakly interacting neutrinos in the SM, retaining both spin-independent and spin-dependent terms. We then derive neutrino-mediated potentials when there are vector, scalar and tensor non-standard interactions present as well as an exotic neutrino magnetic moment. Examining the phenomenology of such long-range potentials in atomic scale laboratory experiments, we derive upper bounds on the Wilson coefficients of the effective operators and compare these to those from processes such as charged lepton flavour violation.


I. INTRODUCTION
Along with the four known types of interaction in nature -the electromagnetic, strong, weak and gravitational -it is possible that there exist additional forces. The exchange of a new spin-zero or spin-1 boson between two fermions may give rise to some exotic new force, e.g., an axion (spin-zero), a dark photon (spin-1) or a light Z are among the potential candidates for such a scenario which might lead to an exotic fifth force, which are being actively searched for at several experiments . Some such forces between elementary particles scale, at large distances, with an inverse power of the distance between the particles; they are often referred to as long-range forces. They may arise from the exchange of a massless mediator between two particles -the inverse-square Coulomb interaction between charged particles being the most common example, arising due to an exchange of a single photon. The prospects of discovering a new long-range force coupling to ordinary matter is highly intriguing from both the theoretical and experimental points of view. In electroweak (EW) theory, the neutrinos are the lightest particles and can be considered as nearly massless on the scale of atoms. The exchange of a single neutrino (and in general any fermion) will change the angular momentum of the exchanging particles involved and therefore cannot give rise to a force between stable matter. However, the exchange of two neutrinos can keep the quantum numbers of the exchanging particles unchanged, and can potentially lead to a long-range force. Historically, the idea of a long-range force mediated by the exchange of two neutrinos has been conceived a long time ago [40][41][42][43][44]. The first explicit computation of the two-neutrino exchange force using the four-Fermi approximation was performed at the leading order in Ref. [41] to obtain a potential of the form V (r) = G 2 F /(4π 3 r 5 ), where G F is the Fermi constant. The neutral-current effects were included in Ref. [42] and the velocity dependence up to first order was included in Ref. [43]. However, in all these calculations neutrinos were assumed to be massless and of a single flavour.
An interesting series of discussions in the literature resulted also from Ref. [44], where it was reported that if neutrinos were massless, the two-neutrino exchange force between neutrons can lead to a large self-energy in a neutron star system through many-body interactions, which far exceeds the order of magnitude of the rest mass of the system itself. It was noted that such a situation can be avoided if the neutrinos are massive, shortening the range of the relevant interaction. In Ref. [45], however, it was argued that the creation and subsequent capture of low-energy neutrinos in the star will fill a degenerate Fermi neutrino sea that can block the free propagation of the neutrinos that are responsible for the neutrino force. In Ref. [46], a calculation of the self-energy using a different technique was reported to obtain a negligible contribution. In Ref. [47] it was stated that the two-neutrino exchange force can be repulsive leading to repulsion among neutrinos instead of filling up the Fermi sea in the neutron star. In Ref. [48], on the other hand, it was argued that the neutrino self-energy does not threaten the stability of the neutron star, to which Ref. [49] differs 1 .
In Ref. [50] the two-neutrino exchange potential was calculated for massive Majorana neutrinos, which was further improved in Ref. [51] to include the effects of flavour mixing for the spin-independent part of the force. In a recent work [52], one can also find a detailed inclusion of flavour mixing effects for the spin-independent part of the neutrino exchange force. In Ref. [53], the second order EW effects have been discussed which is usually ignored in the effective (Fermi theory) approach and the relevant EW second-order shifts are calculated for muonium energy levels. While the first-order EW contribution to the hyperfine splitting in 1S muonium is found to be of order 65 Hz, the second-order corrections are found suppressed by two orders of magnitude, therefore making any new physics corrections of the first-order EW contribution more relevant than higher-order effects in the EW theory.
Recently, atomic and nuclear systems have attracted substantial attention as probes of SM physics and beyond [54][55][56][57][58], with excellent improvements in the experimental precision and a promising future prospect for further improvements on experimental measurements [59]. In particular, in Ref. [55], it was explored how the long-range neutrino exchange force can be probed using atomic and nuclear spectroscopy. In Ref. [58], the possibility of distinguishing between the Dirac vs. Majorana nature of neutrinos is discussed in the context of a violation of the weak equivalence principle, while the same has been explored using the neutrino Casimir force in the plate-plate and point-plate configurations in Ref. [60]. Given the significant amount of progress already made in the literature, it seems desirable to have a robust and systematic analysis of all the possible operator realisations for the long-range neutrino exchange force with equal emphasis on spin-independent and spin-dependent parts. The latter part is particularly relevant given the atomic and nuclear spectroscopy provides sensitivity to both types of longrange neutrino exchange forces, as discussed in Ref. [55]. Furthermore, a clear distinction and comparison of the Dirac vs. Majorana neutrino cases and possible connection of the relevant short-range operators with other relevant observables are also expedient.
In the present work, we consider an effective field theory (EFT) approach to analyse the long-range potential induced by the exchange of two neutrinos in a systematic way, including all the possibilities for the relevant four-fermion contact interactions including the usual SM vector and axial-vector interactions, the scalar and pseudo-scalar interactions and tensor interactions.
The effects of flavour mixing are kept completely general and the possibility of having a righthanded current for the neutrinos are also considered in view of many SM extensions pointing towards such a possibility, see e.g., Refs. [61][62][63][64][65][66][67][68][69][70][71][72][73]. We present both spin-independent and spindependent results for the long-range potential induced by the exchange of two neutrinos and also analyse the effects for a considerable neutrino magnetic moment. We also discuss the possibility of probing spin-independent and spin-dependent components of the long-range potential using state-of-the-art atomic and nuclear spectroscopy experiments. In particular, the muonium atom currently provides the most precise probe providing access to physics at the scale of tens of GeV and is sensitive to the spin-dependent components of the long-range potential, which has prospects of further improvement at J-PARC Muon Science Facility (MUSE) with new highintensity muon beam [59]. In view of the relevant effective operators also inducing charged lepton flavour violating (cLFV) observables, subject to very tight constraints from ongoing and upcoming experiments, we also compare the relevant constraints and comment on their possible complementarity in view of an EFT approach. We also comment on other particle physics probes of these operators, e.g. electron-ν and nucleon-ν scattering, beta decays and ee → ννγ at LEP. This paper is organised as follows. In Sec. II we introduce the low-energy EFT formalism and its connection to the SM gauge-invariant EFT. We also discuss the current bounds on the relevant Wilson coefficients from probes such as cLFV processes, electron-neutrino scattering, nucleonneutrino scattering, beta decays and LEP data. In Sec. III we outline the derivation of a potential associated with the exchange of a virtual particle between two fermions, taking into account their spins. In Sec. IV we derive potentials induced by the exchange of two neutrinos between SM charged-and neutral-current interactions in addition to other non-standard vector, scalar and tensor interactions. We conclude this section by discussing and comparing the potentials in these scenarios. In Sec. V we discuss the prospect of probing beyond the SM effective operators using atomic spectroscopy measurements. We summarise the current experimental measurements and use them to derive limits on the various Wilson coefficients for the Dirac and Majorana neutrino scenarios. In Sec. VI, we discuss the effects of non-vanishing electromagnetic properties of neutrinos and derive the relevant long-range potentials. We also derive the relevant limits on the neutrino electric and magnetic dipole moments using the currently available experimental data. Finally, in Sec. VII, we make our concluding remarks.

A. Low energy EFT
In order to study the effects of new physics interactions of neutrinos in the context of neutrinomediated long-range potentials, we first need to specify the relevant effective field theory (EFT) framework. If the non-standard interactions are the result of some new physics at a high energy scale Λ NP , the general impact of such interactions is to induce operators containing all possible permutations of SM fields respecting the global and gauge symmetries present at a lower scale µ Λ NP , where Λ NP is the cut-off scale of validity for the EFT. This can be written as a series of higher-dimension (d ≥ 5) non-renormalisable operators, where L SM is the SM Lagrangian, O At energies below the EW scale -relevant for the long-range exchange of two neutrinos -the SM gauge group is broken and the operators O This is the so-called low energy effective field theory (LEFT) which has been studied in detail, for example, in Refs. [74][75][76]. In those works a complete basis of operators up to dimension-six is given along with their associated anomalous dimensions, needed to compute the running of the operators from the scale µ up to Λ NP via the renormalisation group equations. Also given are the matching conditions between the LEFT operators and the EFT respecting the SM gauge group (SMEFT) valid at the scale Λ NP . A complete basis and set of anomalous dimensions has also been computed in the SMEFT up to dimension-six [77][78][79][80][81]. However, in general the operators considered in the LEFT can be lepton number violating and all such LEFT operators with d ≥ 6 require SMEFT operators with odd dimension higher than six. The only SMEFT operator at dimension-five is the well-known LNV Weinberg operator [82] where l ρ (ρ = e, µ, τ ) and H are the leptonic and Higgs SU (2) L doublets, respectively, l c ρ = Cl T ρ with the charge-conjugation matrix C andH = iσ 2 H * , where σ 2 is the second Pauli matrix.
Because the LEFT and SMEFT are both constructed out of the SM field content they contain only the left-handed neutrino field ν L , either explicitly for the former and contained in the doublet l for the latter. Technically no assumption is made about the nature of massive neutrinos -whether they are Dirac fermions and have a right-handed component ν R or are selfconjugate Majorana fermions, ν R = ν c L = Cν T L . After the EW symmetry breaking the Weinberg operator in Eq. (2) generates a Majorana neutrino mass term -however, if neutrinos are Dirac and lepton number is strictly conserved then the coefficient C 5 must vanish. 2 The non-standard neutrino interactions relevant to long-range neutrino exchange are those that contain a neutral-current (NC) for the neutrinos (νΓν) and the interacting fermions (f Γf ), where Γ is a product of gamma matrices. We include the right-handed component ν R so that the light neutrinos can be either Dirac or Majorana. The lowest dimension operators containing both neutral neutrino and fermion currents are at dimension-six. In order to compare these with the low-energy Fermi limit of the SM weak interactions, we normalise the Wilson coefficients with respect to the Fermi constant G F . There are ten different Lorentz-invariant operators in the resulting effective Lagrangian, where f = ( , u, d) and the fields are in the flavour eigenstate basis with α, β = e, µ, τ , α, β = u, c, t, α, β = d, s, b, respectively. Likewise, ρ, σ label the neutrino flavours, ρ, σ = e, µ, τ .
So far we have kept the coefficients c XY , g XY and h XX (X, Y = L, R) in the flavour basis of neutrino and fermion fields. The relevant coefficients in the mass eigenstate basis should therefore contain the relevant elements of the CabibboKobayashiMaskawa (CKM) and Pon-tecorvoMakiNakagawaSakata (PMNS) mixing matrices. We will follow the convention that the down-type quark and lepton Yukawa matrices Y d and Y are diagonal and the diagonalisation of the up-type quark and neutrino Yukawa matrices (in the Dirac case) proceed via the bi-unitary where V ,Ṽ , U ,Ũ rotate the left and right-handed up-type quark and neutrino fields according where for clarity the primed and unprimed fields denote flavour and mass eigenstates respectively -the neutrino mass eigenstates are also labelled with the index i. The matrices V and U then correspond to the CKM and PMNS matrices appearing in the SM charged-current, The matricesṼ andŨ do not appear in any SM interaction -while u R is present in the SM NC, the form of the currentū R γ µ u R cancelsṼ , provided it is unitary, and ν R is not present at all in the SM-and so are taken to be unphysical. On the other hand, they will appear for some of the operators in Eq. (3) by rotating the fields to the mass basis.
One can choose to define the coefficients c XY , g XY and h XX (X, Y = L, R) in the mass basis by absorbing the CKM and PMNS matrix elements into the coefficients in the flavour basis. For example, the Wilson coefficient c LL in the mass eigenstate basis is given by where the V γα V * δβ factor is only present for f = u. On the other hand the Wilson coefficient c RR is written in the mass basis as which contains rotation matrices for the right-handed neutrino fields (and possible right-handed up-type quark fields). We immediately see a redundancy in Eq. (9) because there is more than one unknown parameter on the right-hand side. The unknown mixing angles and phases in the mixing matricesṼ andŨ can instead be absorbed back into the parameters of the matrix c RR in the flavour basis -this is equivalent to settingṼ =Ũ = I from the outset. However, we will see thatŨ may contain information about the presence of additional sterile states in the model.
In the SM, the values of the coefficients c LL and c RL are given in the mass basis as where W /3 for s 2 W = sin 2 θ W and θ W is the weak mixing angle. In models such as the type-I seesaw that can introduce additional mass eigenstates, therefore making the 3 × 3 PMNS mixing matrix non-unitary, it is easy to replace δ ij → C ij , where Here U αj corresponds to the generalised PMNS mixing matrix.
To go from effective coefficients at the quark level (e.g. c LL uu;ij and c LL dd;ij ) to the level of nonrelativistic nucleons one must make use of the chiral EFT -viable when the relevant momentum exchange is below the cut-off scale of the EFT, Λ ChEFT ∼ O(1 GeV). At leading order in the EFT, the light pseudoscalar masses are of order m π ∼ O(q) and the neutrinos in the loop only interact with a single nucleon. Interactions of neutrinos with more than one nucleon (for example in deuterium) are suppressed by powers of q/Λ ChEFT .
Following the approach of Ref. [85], the coefficients for effective operators containing nucleon currents can be written in terms of the quark-level coefficients as where the sum is over q = u, d, s and F q/N 1 (q 2 ) and F q/N A (q 2 ) are the NC vector and axial vector form factors for the quark q within the nucleon or nucleus N , respectively. For the proton the following linear combinations at zero-momentum exchange are given in the SM by where where We have used that the vector form factors for the valence quarks in the deuteron are F Here ∆s is the strange axial moment of the deuteron, m D is the deuteron mass, E D is the deuteron binding energy, m π is the neutral pion mass, and f π is the pion decay constant. The renormalisation scale µ is taken to be at the neutral pion mass m π .

B. Inclusion of sterile neutrinos
We now briefly return to the question of matching the LEFT + ν R operators with a SM gauge-invariant EFT. As stated before, the commonly-studied SMEFT does not contain ν R and so can only produce a subset of the LEFT + ν R operators containing just ν L . It may not even be possible to match these operators if there are new particles with masses O(100) GeV or if EW symmetry breaking is non-linear in the new physics sector [103]. Assuming that the matching is possible, in order to produce all low-energy operators, one can introduce n number of sterile states N R to the SMEFT. To this end, a complete basis of lepton number and baryon number conserving and violating operators has been considered in the literature at dimension-five [104], dimension-six [105] and dimension-seven [106,107]. Ref. [86] provides the matching conditions between the LEFT + ν R (in the basis of Appendix A) and the SMEFT + N R .
At dimension-four in the SMEFT + N R there are the terms where M R is the Majorana mass for the sterile states and we have the usual Yukawa term. At dimension-five, in addition to the Weinberg operator of Eq.
(2) we have [104] L (5) where the former is an EW coupling ζ to the U (1) Y field strength operator B µν and the latter is a Majorana mass-like coupling χ to Higgs doublets. After EW symmetry breaking we obtain the following mass terms quasi-Dirac neutrinos such as those studied in Refs. [108,109]. In the former case the mass matrices for the light and heavy neutrinos are approximately respectively, which are then diagonalised to the mass basis via M ν = U T ν · M ν · U ν and M N = U T N · M N · U N . In the mass basis we therefore obtain three light Majorana neutrinos ν = ν c and n relatively heavier Majorana neutrinos N = N c . The weak eigenstates are given in terms of the mass eigenstates by where n = (ν 1 , ν 2 , ν 3 , N 1 , N 2 , ...). As all the light neutrinos are Majorana in this case, the relevant operators with coefficients c LL and c RL (the c LR and c RR operators are equivalent by Eq. (4)) are naively rotated by the 3 × (3 + n) matrix U . However, if any of the N are above the EW scale then they are integrated out of the EFT, leaving only those below. In the presence of light sterile states one can have terms such as where s, s label the weak eigenstates of N R , which are rotated to the mass basis by the n×(3+n) U corresponds to the right-handed analogue of the PMNS mixing matrix diagonalising the right-handed charged-current interactions [71,72]. However, we have already discussed how it multiplies similarly unknown Wilson coefficients and can thus be subsumed.
In this work, we are mainly concerned with the two extreme limits discussed below Eq. (18).
The first is the scenario in which lepton number is conserved -M L = M R = 0 and the righthanded states N R ≡ ν R form three light Dirac neutrinos with the ν L states. The second is if lepton number is violated and any sterile mass eigenstate fields N are integrated out, leaving three light Majorana neutrinos. However, we will see that light sterile fields N (which may or may not be related to the seesaw mechanism generating light left-handed neutrino masses) can be of relevance in Sec. VI.

C. Bounds from other probes
In the SM, the left-handed neutrinos are part of SU (2) L doublet with the charged leptons as their partners. Therefore, for a given new physics model at a scale higher than the EW symmetry breaking scale, the SM gauge-invariant operators that mediate the long-range neutrino interactions can also mediate charged lepton flavour changing processes [110] which are constrained stringently from experimental non-observation of various cLFV observables [111][112][113].
The cLFV radiative decays and µ → e conversion are particularly relevant in this context as they are subject to intensive searches at various ongoing and upcoming experiments. The decays of tau into light mesons accompanied by a lepton are also relevant since the relevant bounds are expected to be improved significantly in Belle II [114]. To illustrate the relevance of cLFV processes let us consider the contact interaction of Eq. (3) For f L = e L , this operator can be generated by the dimension six SMEFT [77][78][79][80][81] operator wherel is the SU(2) doublet (ν L , e L ). It is easy to see that such an operator can also induce the cLFV interaction (e ρ γ µ P L e σ )(eγ µ P L e) (where ρ = σ) with the relevant Wilson coefficient constrained by the experimental limits on cLFV decays, e.g., e σ → e ρ eē. However, such a scenario can be avoided by instead constructing such operators at higher dimension, e.g., dimension eight where pQ is the usual antisymmetric SU(2) contraction and Λ NP is the heavy new physics scale.
where X (X )corresponds to the bases of operators with Lorentz structure γ µ × γ µ , and ζ corresponds to the flavour indices ρσf f . A complete list of all the relevant dimension six operators in the "Warsaw" basis can be found in [78], while the relevant dimension eight operators when the external fermion is a SU (2) L doublet can be found, for example, in [115]. A particularly interesting basis has also been proposed recently in [113]; following this basis in the case where the external fermion (f ) is a SU (2) L doublet quark or lepton (f = q, l α with α = ρ, σ), the relevant dimension six SMEFT operators are where the SU (2) L contractions are understood to be inside the parentheses. At dimension eight, the relevant SMEFT operators are given by In the case when the external fermion f = l α with α = ρ, the relevant operators reduce to In the case where the external fermion (f ) is SU (2) L singlet quark or lepton, the relevant dimension six and eight SMEFT operators are given by After the EW symmetry breaking, at Z-pole (µ = m Z ) the SM gauge group invariant operators are matched onto the LEFT operators given in Eq. (3). However, it is important to include the RG running induced mixings via W and Higgs loops exchange between various SMEFT operators discussed above. This is particularly relevant because even if the relevant interactions in a given new physics model may not induce cLFV at the tree level, such operators get induced via the mixings at the one-loop level. A detailed discussion of the relevant matching and mixing effects is beyond the scope of the current work and can be found for example in [113]. In Table   I   Other than the cLFV processes discussed above, the Wilson coefficients relevant to first-and second-generation leptons are also subjected to direct bounds from the experimental data on various scattering processes such as ν µ e scattering in CHARM-II [116,117] (which is supposed to be improved by an order of magnitude at the DUNE near detector [118]), neutrino-nucleon scattering data at CHARM and CDHS [116,117,119]. The Wilson coefficients relevant to tau are constrained from eē → ννγ data at LEP [120]. However, these bounds are orders of magnitude weaker as compared to the bounds from cLFV processes. Some relevant discussion about these bounds can be found for example in [120][121][122][123]. In addition, the observation of coherent elastic neutrino-nucleus scattering at COHERENT [124,125] and beta decays [126] are also relevant for deriving bounds on the relevant Wilson coefficients [86].

III. LONG-RANGE POTENTIALS AND SCATTERING AMPLITUDES
Within Quantum Field Theory (QFT), a force acting between particles can be interpreted as the exchange of a virtual particle (or multiple particles) between external on-shell states. As depicted to the left of Fig. 1, a mediator or mediators are necessary to exchange the momentum q = p α − p α = p β − p β between the two interacting particles f α and f β with initial momenta p α and p β and final momenta p α and p β respectively.
In the Feynman diagrammatic approach it is possible to derive a long-range potential V (r, v) for an interaction -most generally a function of the relative displacement between the particles r and the average velocity of the system, v = 1 2 by taking the Fourier transform of the invariant amplitude of the scattering process [127], i.e., where the invariant amplitude M(s, t) is an analytic function of the Mandelstam variables The potential is time independent in the static limit of momentum transfer, q ≈ (0, q) and t ≈ −q 2 , which is an accurate approximation for particles interacting at a distance. Furthermore, one can also exploit the analyticity properties of M(s, t) which enable the spectral decomposition [42] M where ρ(s, t ) is the so-called "spectral function" of the process. The spectral function is related to the imaginary part of the discontinuity on the real t-axis of M(s, t), where for → 0. One can now insert the decomposition of Eq. (32) into Eq. (31) and evaluate the angular integral dΩ = dθ dφ sin θ contained in d 3 q. This integration is non-trivial if M(s, t) and therefore ρ(s, t) depend on θ and φ -for example if there are spin-dependent terms containing the dot product of q and a particle spin σ. In fact, such terms arise naturally when taking the non-relativistic limit of the scattering amplitude.
We follow the approach of Ref. [128] and divide the spectral function ρ(t) (omitting the dependence on s) according to a basis of 16 spin operators, where f k (v 2 ) are polynomials in powers of v 2 corresponding to higher order terms in the nonrelativistic expansion. The operators O k form a complete basis constructed from the relevant three-momenta (q and P) and the interacting particle spins (s α = σ α /2 and s β = σ β /2), where (α, β) is a shorthand for (α ↔ β).
Combining Eqs. (31), (32) and (35), the potential can also be split up as where The functions V k (r, v) can be computed by first evaluating the integral without the factor O k and multiplying by a single power of r, where for the second equality we have integrated over q, θ and φ and relabelled the dummy variable t as t. As outlined in Ref. [128], the functions V k (r) can be readily computed by applying derivatives to the V k (r) functions. We have for example the following operations for

IV. LONG-RANGE POTENTIALS FROM TWO-NEUTRINO EXCHANGE
In this section we will derive, using Eq. In Sec. IV A we derive the potential V LL αβ (r) when only the SM CC and NC interactions are present. In Sec. IV B we include right-handed vector-type neutrino currents and derive the potentials V LR αβ (r) and V RR αβ (r) when one or both of the neutrino currents are right-handed. In Sec. IV C we introduce scalar interactions and derive the vector-scalar and scalar-scalar potentials V V S αβ (r) and V SS αβ (r). In Sec. IV D we consider tensor interactions, determining the vector-tensor potential V V T αβ (r). We derive each potential for Dirac and Majorana neutrinos, examining the dependence on the distance in the short-and long-range limits and on the spins of the external states. We finally plot and compare the potentials in Fig. 4 of Sec. IV E.

A. Standard Model charged and neutral currents
We begin by deriving the potential V LL αβ (r) arising from the SM diagrams in Fig. 2. For simplicity we determine the amplitude M αβ (and the corresponding spectral function ρ αβ ) by integrating out the W ± and Z boson propagators and using the LEFT interaction Lagrangian of Eq. (3). We see that W ± exchange can only occur for charged leptons while Z exchange is possible for both leptons and quarks within a nucleon/nucleus N . Both W ± and Z exchange contribute to the coefficient c LL αβ;ij while only Z exchange contributes to c RL αβ;ij -the values for these are given in Eq. (10). The external fermion currents are therefore either left-or righthanded while the neutrino currents are strictly left-handed.
Applying the appropriate Feynman rules from the interaction Lagrangian of Eq. (3), we can write the invariant amplitude of the scattering process in Fig. 1 where 1/(4m α m β ) is a normalisation factor convenient in the non-relativistic limit [42]. The amplitude firstly contains the sum over the neutrino mass eigenstates, i, j, which we allow to run from 1 to N = 3 + n to allow for the presence of n additional Dirac or Majorana states. It also contains the sum over the possible chiralities (X, Y = L, R) of the external fermion currents. As we are focussing on scattering processes in which the flavours of the interacting fermions do not change, the coefficients c XY αβ;ij will always be diagonal in the flavour of the external fermions (α = β). We therefore relabel c XY αα;ij ≡ c XY α;ij in Eq. (41) and the following discussion. The amplitude in Eq. (41) is also split conveniently into two Lorentz tensor factors. The first is the product of external fermion bilinears where u sα (p α ) and u s β (p β ) are four-component Dirac spinors for the fermions f α and f β (or nucleon N ) and P X and P Y are the usual chirality projection operators P L = (1 − γ 5 )/2 and The second factor N µν ij integrates the product of massive neutrino propagators over the loop momentum k, We now use the method from Eq. (37) onwards to calculate the potential. Using Eq. (33) we first determine the spectral function by taking the discontinuity the amplitude. The discontinuity only needs to be taken for the neutrino loop factor, where Θ(x) is the Heaviside step function, m 2 ij = (m 2 i + m 2 j )/2 is the average of the squares of the neutrino masses, ∆m 2 ij = m 2 i − m 2 j is the difference in the squares and Λ(x, y, x) is the Källén function, To compute the spectral function we contract H αβ µν with disc(N µν ij ). The Lorentz indices of H αβ β . An assumption we now make is that the external fermions are non-relativistic. In this limit with the lowest-order terms in the non-relativistic expansion -Appendix B lists the lowest-order terms for bilinear products such as The terms that dominate are proportional to 4m α m β , cancelling the 1/(4m α m β ) normalisation factor in the amplitude. Higher-order terms in the expansion are suppressed by powers of q/m α and can be neglected.
The discussion has so far been valid for Dirac neutrinos. For Majorana neutrinos only the axial part contributes to the left-handed current and is a factor two larger than the Dirac axial vector current. The neutrino loop factor N µν is instead given by where an additional factor of 1/2 is required due to the permutation symmetry of the Majorana states in the loop.
Using Eq. (33) we can now write the spectral function as where we choose the LL superscript to indicate the presence of two left-handed neutrino currents. Now inserting Eqs. (44) and (42) into Eq. (47) and taking the non-relativistic limit, we obtain to lowest-order which retains one spin-independent and two spin-dependent terms. The factors X LL αβ;ij and Y LL αβ;ij are given by the following combinations of the LEFT coefficients, where we have used the shorthand notation (c LL ± c RL ) α;ij ≡ c LL α;ij ± c RL α;ij . Taking the SM values of the coefficients in Eq. (10) and assuming a unitary light neutrino mixing matrix U αi such that C ij = δ ij , these factors for leptons (α, β = e, µ, τ ) are for example The functions F D ij , F M ij and F V ij in Eq. (48) are given by The difference between the Dirac (D) and Majorana (M ) cases is reflected in the function F To determine the overall potential V LL αβ (r) we evaluate the integral in Eq. (39) for each of the three parts of the spectral function multiplying these operators. We then take the appropriate derivatives in Eq. (40) to derive the three components of the potential V LL k (r) (k = 1, 2, 3) and add these to obtain wherer = r/|r| is the unit displacement between the interacting states and the integral functions We define these functions to be dimensionless in order to take the dimensionful factor G 2 F /(4π 3 r 5 ) out of the sum. The potential therefore scales naively as 1/r 5 though we will see that this behaviour changes in the long-range limit. The difference between the Dirac and Majorana cases is now a difference in the functions The neutrino-mediated potential in Eq. (52) simplifies when only a single massive neutrino is considered. Firstly, the mixing factors in X LL αβ;ij and Y LL αβ;ij are replaced as U * αi U αj → 1 and δ ij → 1 and the summation is now over a single state i = j ≡ ν. For the potential between two charged leptons we have for example X LL αβ;νν = (1 + g V ) 2 and Y LL αβ;νν = (1 + g A ) 2 . Secondly, the functions I D ij (r), I M ij (r) and I V ij (r) take the closed-forms where K n (x) are modified Bessel functions of the second kind and G m,n p,q is the Meijer G-function. Using the relations in Appendix C we can also determine the functions J D νν (r), J M νν (r) and J V νν (r). For interacting leptons the spin-independent parts of the Dirac and Majorana potentials become respectively, in agreement with previous results [52].
The functions in Eq. (53) depend on the product 2m ν r -given the behaviour of the modified To verify this quantitatively we expand the functions in Eq. (53) and therefore the singleneutrino potential in the opposing limits. For r r ν we find to lowest-order in both the Dirac and Majorana cases, as expected. The potentials therefore decrease with the distance as 1/r 5 up to half the neutrino Compton wavelength. Eq. (55) is not just valid for a single neutrino -it can be obtained for three (or in general, N ) neutrinos by neglecting the neutrino masses m i and m j appearing in the functions I X ij (r) and J X ij (r) in Eq. (52). In this limit the functions tend to the constant values as outlined in Appendix C. It is now possible to pull these constants out of the sum in Eq. (52) and identify X LL αβ;νν = N i,j X LL αβ;ij and Y LL αβ;νν = N i,j Y LL αβ;ij . Expanding in the opposite limit r r ν gives in the Dirac case while in the Majorana case In the Dirac case both spin-independent and spin-dependent terms scale as e −2mν r /r 5/2 , while in the Majorana case the spin-independent and σ α · σ β terms scale as e −2mν r /r 7/2 . The term containing (σ α ·r)(σ β ·r) however also scales as e −2mν r /r 5/2 in the Majorana case.
Finally, we note that the operators with coefficients c LL α;ij and c RL α;ij may include the effects of new physics, which can be parametrised as small corrections δc LL α;ij and δc RL α;ij to the SM values of c LL α;ij and c RL α;ij . Deviations from the SM potential V LL αβ (r) therefore arise as corrections to the factors X LL αβ;ij and Y LL αβ;ij , where the correction can either be at the vertex with fermion f α or f β .

B. Right-handed vector non-standard interactions
Motivated by theories such as the Left-Right Symmetric Model (LRSM) we now introduce a right-handed neutrino current. We will first derive the neutrino-mediated potential V LR αβ (r) induced when there is a SM CC or NC interaction at one vertex and a right-handed neutrino current at the other, depicted in Fig. 3. In the LEFT interaction Lagrangian of Eq. (3) we now allow the coefficients c LR αβ;ij and c RR αβ;ij to be non-zero along with c LL αβ;ij and c RL αβ;ij . The spectral function ρ LR αβ (t) in this scenario is the same as Eq. (48) but with one coefficient replaced as c XL α;ij → c XR α;ij and one chirality projection operator replaced as P L → P R in the neutrino loop factor N µν ij . We also add an identical contribution with (α ↔ β) to account for the right-handed current being either at the vertex with the external fermion f α or f β . If the external fermions are identical (α = β) we must multiply the spectral function by an additional factor of 1/2 to avoid double counting -this gives a factor 1/(1 + δ αβ ). The discontinuity of the neutrino loop factor N µν ij for Dirac neutrinos is which is suppressed by the factor m i m j /q 2 because the process is helicity-suppressed. A negative helicity neutrino ν i created by the left-handed current will be annihilated by the right-handed current with an associated factor m i /q -for both neutrinos this results in the m i m j /q 2 factor.
The physics is identical to the helicity-suppressed contribution to Majorana neutrino exchange in the previous subsection.
Contracting the g µν factor in Eq. (60) with the product of external fermion bilinears H αβ µν we obtain [γ µ P X ] α [γ µ P Y ] β . The non-relativistic limit can now be taken to obtain the spectral function -in the Dirac case where the pre-factors X LR αβ;ij and Y LR αβ;ij are Using the same method as the previous section to derive the potential from the spectral function, we find in the Dirac case where the dimensionless function I LR ij (r) is given in Appendix C The Majorana case is different due to the symmetry relations of Eq. (4) -the right-handed current operator with coefficient c LR αβ;ij is equivalent to the left-handed current operator with coefficient c LL αβ;ij and thus the coefficients are related by c LL αβ;ij = −c LR αβ;ji . This is equivalent to the vector current vanishing for Majorana neutrinos. The potential we derive from the righthanded current operator is therefore identical to Eq. (52) and the coefficient c LR αβ;ij gets the same contributions from the SM CC and NC interactions as c LL αβ;ij . If on the other hand we were to introduce additional light sterile Majorana states N R with right-handed interactions as in Eq.
where δX LR αβ;ij = −δX LL αβ;ji and δY LR αβ;ij = −δY LL αβ;ji are given in Eq. (59). This gives the correction to the SM potential We again us the single neutrino simplification to study the properties of the potentials in Eqs. (63) and (65) -the function I LR ij (r) takes the closed form In the short-range limit, or r r ν , we expand the Dirac potential Eq. (63) as This potential is also valid in the three (or N ) neutrino picture by using that I LR ij (r) ≈ 1 in the limit m i ≈ 0 and identifying X LR αβ;νν = In the single neutrino simplification for the correction in Eq. (65) to the SM Majorana potential takes the same form as Eqs. (55) and (58) in the short-and long-range limits respectively.
We finish this subsection by considering the case where there are two right-handed neutrino currents at the interaction vertices. Now the potential takes the same form as Eq. (52), where Consequently, the short-and long-range potentials are given by Eqs. (55) and (58) respectively with the replacements X LL αβ;ij → X RR αβ;ij and Y LL αβ;ij → Y RR αβ;ij .

C. Scalar non-standard interactions
We now derive the neutrino-mediated potential in the presence of a scalar non-standard interaction. In our framework these are the operators in Eq. (3) with the coefficients g LL , g RL , g LR and g RR normalised to the Fermi constant G F . We first focus on the case of a scalar interaction at one vertex and a SM CC or NC interaction at the other, as shown Fig. 3.
The spectral function can be determined in this scenario according to where the sum is over the possible chiralities of the external fermion currents (X and Y ) and the neutrino current of the scalar operator (Z). We have taken into account that the scalar interaction may be at either vertex by adding an identical contribution with (α ↔ β). The Majorana case is treated in the same way by retaining only twice the axial vector current and dividing by a factor of two due to the permutation symmetry of the neutrinos in the loop.
The discontinuity of N µ ij in the Dirac case for example is where the minus (positive) sign is for a left-handed (right-handed) neutrino current at the scalar interaction and the product of external fermion bilinears is Contracting H αβ µ with disc(N µ ij ) and making use of the non-relativistic limits of the fermion bilinear products given in Appendix B, we obtain in the Dirac case while in the Majorana case we obtain The factor X V S αβ;ij containing the scalar coefficients is and the functions F ∆ ij and F S ij are given by These spectral functions only contain terms proportional to the parity-violating spin operators O 9 = σ α · q and O 10 = σ β · q proportional to linear combinations of the operators O 9 and O 10 in Eq. (36). Taking the components of the spectral function multiplying these operators we compute the functions V 9 (r) and V 10 (r) in Eq. (39) and from these the components of the overall potential using from Ref. [128].
Thus we derive the following vector-scalar potentials for the Dirac and Majorana cases respectively, where the dimensionless functions J ∆ ij (r) and J S ij (r) are given in Appendix C. The first thing to note about these potentials is that they depend on the distance as 1/r 4 and contain a single power of the neutrino mass m ν in the numerator. This is more suppressed than the SM-SM potential in Eq. (52) which scales as 1/r 5 but less suppressed than the right-handed current potential for Dirac neutrinos in Eq. (63) which scales as 1/r 3 but is suppressed by m 2 ν . The second point to note is that the potentials written in Eqs. (81) and (82) retain a factor of i -this can simply be absorbed into the factor X V S αβ;ij after a suitable redefinition of the scalar coefficients g XY α;ij . We now consider the case where both interactions are scalar. We now obtain the potential via the spectral function where the discontinuity of the neutrino loop factor N ij is given for example in the Dirac case by if the chirality of the neutrino currents are the same (X = Z = L, R) and if the chiralities of the neutrino currents are opposite (X = Z). The external fermion bilinear product is now H αβ = [P W ] α [P Y ] β and we obtain the following potential in the Dirac case where the combination of scalar coefficients are given by X SS αβ;ij = (g LL + g RL ) α;ij (g LL + g RL ) β;ij + (g LR + g RR ) α;ij (g LR + g RR ) β;ij , Y SS αβ;ij = (g LL + g RL ) α;ij (g LR + g RR ) β;ij + (g LR + g RR ) α;ij (g LL + g RL ) β;ij , and the dimensionless functions I SD ij (r) and I LR ij (r) are given in Appendix C. For Majorana neutrinos we instead obtain where the combination of scalar coefficients is and the function I SM ij (r) is also given in Appendix C. We see that the Dirac potential depends on the distance as 1/r 5 only when the neutrino currents are of opposite chirality -when they are the both left-or right-handed the potential becomes suppressed as m 2 ν /r 3 . This suppression does not occur for Majorana neutrinos -the potential scales as 1/r 5 for any combination of the coefficients g XY α;ij .

D. Tensor non-standard interactions
We now derive the neutrino-mediated potential in the presence of a tensor non-standard interaction. In our framework these are the operators in Eq. (3) with the coefficients h LL and h RR normalised to the Fermi constant G F . We first focus on the case of a tensor interaction at one vertex and a SM CC or NC interaction at the other, as shown Fig. 3.
The spectral function can be determined in this scenario according to where the sum is over the possible chiralities of the external fermion currents (X and Y ) and the neutrino current of the tensor operator (Z). We have again taken into account that the tensor interaction may be at either vertex by adding an identical contribution with (α ↔ β).
The Majorana case is treated in the same way as previous subsections.
The discontinuity of N µνρ ij in the Dirac case is and the external fermion bilinear product is H αβ Contracting these factors and using the non-relativistic limits in Appendix B, we obtain the spectral function in the Dirac case while in the Majorana case we obtain The coefficients X V T αβ;ij and Y V T αβ;ij containing the dependence on the tensor coefficients are The function F ∆ ij is the same as Eq. (78) and F T ij is given by The spectral functions above contain terms proportional to the parity-violating spin operators O 9 = σ α · q, O 10 = σ β · q and O 11 = (σ α × σ β ) · q. We can again take the components of the spectral functions multiplying these operators and evaluate the functions V 9 (r), V 10 (r) and V 11 (r) of Eq. (39). From these we use Eq. (80) and to derive the full vector-tensor potential in the Dirac case and in the Majorana case where the dimensionless functions J ∆ ij (r) and J T ij (r) are given in Appendix C. We note that these potentials, like the vector-scalar potentials of the previous section, scale as m ν /r 4 . They are similarly contain only parity-violating spin operators. The difference between the potentials for Dirac and Majorana neutrinos arises from the different r-dependence of the functions J ∆ ij (r) and J T ij (r). Finally, we see that the diagonal elements in the i, j sum vanish for Majorana neutrinos.

E. Comparison of potentials
In Fig. 4 we compare a selection the potentials derived in the previous subsections. To the left of Fig. 4 we plot the spin-independent parts of the vector-vector potentials V LL αβ (r), V LR αβ (r) and V RR αβ (r) for positronium (e − e + ) and for either three Dirac or Majorana neutrinos. The potential V LL ee (r) is calculated using SM values for the factors X LL ee;ij and Y LL ee;ij in Eq. (50). The potential V LR ee (r) has a single SM vertex and is interpreted as a correction to V LL ee (r) in the Majorana case, though we plot it separately. The potential V RR ee (r) is derived from two non-standard right-handed neutrino currents. We set m 1 = 0.1 eV and take normal-ordered (NO) values of the mixing angles θ 12 = 33.8 • , θ 23 = 48.6 • and θ 13 = 8.6 • , the CP phase δ = 221 • and the mass splittings ∆m 2 21 = 7.55 · 10 −5 eV 2 and ∆m 2 31 = 2.5 · 10 −3 eV 2 . We set the non-standard coefficients to be c LR ee;ij ≡ c LR e δ ij = 10 −2 δ ij , i.e. only non-zero for diagonal i, j. We first note the small difference between Dirac and Majorana potentials V LL ee,D (r) and V LL ee,M (r) . In the short-range limit r 1/2m 1 the potentials are identical while in the long-range limit r 1/2m 1 the Majorana potential is slightly smaller the Dirac potential, in agreement with the results of Ref. [58]. The potentials are generally seen to fall off as 1/r 5 until the neutrinos become non-relativistic around r ∼ 1/2m 1 and the potentials are exponentially suppressed.
We see that the potentials are many orders of magnitude smaller than the gravitational potential V grav ee (r) between the electron and positron. We also note the large difference between the Dirac and Majorana potentials V LR ee,D (r) and V LR ee,M (r) -while the Dirac potential is slightly larger than the Majorana potential in the long-range limit, in the short-range limit the former scales as 1/r 3 and is suppressed by two powers of the neutrino masses while the latter scales as 1/r 5 and is unsuppressed. This is because the Majorana potential is interpreted as a correction to the SM potential V LL ee,M (r) and thus scales in the same way. V LR ee,M (r) is around two orders of magnitude smaller than V LL ee,D (r) ≈ V LL ee,M (r) due to the suppression from c LR e = 10 −2 . The potential V RR ee (r) is shown just for the Dirac case -because it contains two factors of c RL e = 10 −2 it is seen to be below V LR ee (r). To the right of Fig. 4 we plot the scalar-scalar potentials for Dirac and Majorana neutrinos V SS ee,D (r) and V SS ee,M (r) and compare them to the spin-independent part of V LL ee (r) and the gravitational potential V grav ee (r). We choose a scalar coefficient g XY ee;ij ≡ g XY e δ ij = 10 −2 δ ij to be non-zero for a single choice of the chiralities X, Y -looking at Eqs. (86) and (89) we see that the surviving terms of the Dirac potential scale in the short-range limit as 1/r 3 while for the Majorana potential as 1/r 5 , as can be seen in the diagram.

V. ATOMIC SPECTROSCOPY
There are a number of ways to probe exotic long-range forces over a range of distances.
Starting at the macroscopic scale, precision torsion balance experiments adopt the method originally used by Cavendish to measure the gravitational constant G. Theories looking to resolve the discrepancy between the observed dark energy density ρ d ≈ 3.8 keV/cm 3 and the theoretical prediction from QFT (a factor of ∼ 10 120 larger) predict Yukawa violations or powerlaw modifications of the gravitational force at length-scales of r ∼ 1 µm − 1 mm [129]. These and other theories involve extra time [130] and space [131] dimensions and new scalar and vector mediators such the axion [132], dilaton [133], dark photon and Z [134], all of which can alter the typical 1/r scaling of the gravitational potential and break the weak equivalence principle.
As can be seen in Fig. 2, the neutrino-mediated potentials fall off exponentially for r 1 µm, roughly corresponding to the Compton wavelength of the lightest neutrino with m 1 = 0.1 eV.
For point sources such as an electron and positron the associated forces are many orders of magnitude smaller than their gravitational attraction. In theory this can be overcome by using neutral aggregate matter with a coherent weak charge, boosting the strength of the neutrinomediated force with respect to the gravitational force [58]. It remains to be seen if torsion balance experiments can overcome the strong effect of the Earth's gravity to measure this. Another method is to measure the pressure exerted on two parallel plates by the Casimir force induced by the neutrino potential [60]. Current experiments are however ∼ 20 orders of magnitude below the required sensitivity to measure the neutrino contribution.
To attain a greater sensitivity to the neutrino-mediated potentials one must therefore go to smaller distances where the potentials can be seen to exceed the gravitational potential in Fig. 4.
The most stringent measurements come from nuclear and atomic spectroscopy probing r ∼ 1 fm and r ∼ 1Å, respectively. We outline some of the methods explored in the literature.
Atomic spectroscopy of heavy atomic species (Z 1) might appear to be the most suitable method for probing the spin-independent part of the neutrino-mediated potential thanks to the coherent scaling of the nucleus -going up roughly with the number of neutrons N 1. The spin-dependent part on the other hand acts incoherently because nuclear pairing interactions leave the ground state nucleus with at most two unpaired nucleon spins. However, the complexity of many-electron interactions in heavy atoms makes the theoretical predictions for transitional frequencies inadequate for the current experimental precision. One can instead measure the isotope shift -the difference in atomic splittings for different isotopes -in systems such as Ca + by observing a non-linearity in the King plot [149]. This has been used to constrain models with Z bosons, exotic Higgs bosons and chameleon particles [150][151][152][153] and more recently the neutrino-mediated potential [55].
A relevant probe at nuclear length scales is the binding energy of the deuteron D + , a bound state of a proton and a neutron. One can model the binding energy with a spherical potential well with an infinitely repulsive inner hard core in order to find the radial wave-function of the system. This in turn can be used to calculate the expectation value of the neutrino-mediated potential and the shift to the binding energy. The difference in the measured [154] and predicted [155,156] binding energies has been used to constrain the neutrino-mediated potential [55].
The sensitivity of simple atomic-like systems such as positronium (e − e + ) and muonium (e − µ + ) to the neutrino-mediated potential may be more promising than the deuteron and other nuclear probes. As we will see, for these leptonic systems the characteristic cut-off scale (below which the r-dependence of the potential no longer holds) is provided by the cut-off of validity of the EFT, not the charge-radius of the nucleon or nucleus in semi-leptonic systems like hydrogen (e − p), deuterium (e − D + ) or their muonic counterparts (µ − p and µ − D + ). At present the best measured splittings of these systems are the 1S − 2S and ground state hyperfine splittings.
These splittings have also been predicted to high accuracy and used as precision tests of QED.
For example, the dominating Dirac, radiative, recoil and radiative-recoil QED corrections to the Fermi expression of the ground state hyperfine splitting E F have been calculated up to orders [157][158][159]. Smaller weak [53,160] and hadronic corrections [161] have also been calculated. The EW corrections have been calculated for the muonium hyperfine splitting to next-to-leading-order [53].
We will follow the same approach as Ref. [55] which derives the shifts to energy level splittings using the expectation value of the position-space potential V (r). Using the experimental and SM-predicted values for the 1S − 2S and hyperfine splittings of positronium and muonium, we will use the predicted shifts from the exotic neutrino-mediated potentials to put upper bounds on the non-standard coefficients c XY , g XY and h XY .

A. Shifts to energy levels
The small shift to an atomic energy level due an exotic force can be calculated to first order in perturbation theory by taking the expectation value of the associated potential V (r), where n 2S+1 L J labels the unperturbed atomic state with n the principal quantum number n, S the total spin, L = {S, P, D, ...} the total orbital angular momentum and J the total angular momentum. Shifts to the 1S − nS and n-hyperfine splittings are respectively The average of the potential over the atomic quantum numbers is the position-space integral where Ψ n, ,m (r) is the atomic wave-function. For the two-body systems we are considering, Ψ n, ,m (r) = R n, (r)Y ,m (θ, φ) is the separable solution of the Schrödinger equation with the Coulomb potential V C (r) = −Zα/r.
We will be comparing the shifts induced by exotic potentials depending differently on r. For example, the SM CC and NC induced potential V LL αβ (r) in Eq. (52) scales as 1/r 5 in the shortrange limit, while the right-handed current induced potential V LR αβ (r) for Dirac neutrinos scales as 1/r 3 . Assuming that V (r) is only a function of r (and not θ and φ) the integration over the spherical harmonic component Y ,m (θ, φ) is unity and the average over the hydrogen-like radial wave-function R n, (r) for general r-dependence is, where r c is a lower cut-off on the distance corresponding to an upper cut-off scale of validity for the four-fermion EFT. For SM CC and NC interactions this distance is around the inverse Z boson mass and we define r c = 1/m Z = 1.097 × 10 −11 eV −1 . We can write the Fermi coupling in terms of this length scale using This distance scale could be different for a non-standard effective interaction mediated by a particle with a mass above or below the EW scale -a Z for example. In this case the distance cut-off is r c = 1/m Z . This mediator may also interact with the SM via a coupling g . Comparing this to the normalisation of the effective interaction to the Fermi coupling, Depending on whether the new physics is above or below the EW scale, or strongly or weakly coupled, the lower distance scale of validity r c compares to the SM Fermi cut-off r c as While this discussion is valid for an EFT with point-like particles, for a semi-leptonic system the cut-off r c must take into account the finite size of the nucleon or nucleus -e.g. for a proton We can now integrate Eq. (106) using the hydrogen-like radial wave-function, where L j k (x) is the associated Laguerre function andã 0 is the reduced Bohr radius of the system with reduced mass m r ,ã For hydrogen this is the standard Bohr radiusã 0 ≈ a 0 = 1/(m e α). For different values of d in Eq. (106) and expanding in r c we obtain Here the parameter A n is given by with To compute the average in Eq. (105) we must also take the angular average of the spindependent terms in V (r) -for example the factors σ α · σ β and (σ α ·r) (σ β ·r) in V LL αβ (r) and V LR αβ (r). Firstly, as we will be only considering n 2S+1 S J states for the 1S − nS and n-hyperfine splittings, the following equality holds for = 0 In order to determine the hyperfine splitting between singlet and triplet configurations of external particle spins we must also evaluate the spin dot-product σ α · σ β s in these cases. These are σ α · σ β s=0 = −3 (singlet) and σ α · σ β s=1 = 1 (triplet).
The averages of the parity-odd potentials V V S αβ (r) and V V T αβ (r) -which depend on the spin operators σ α ·r, σ β ·r and (σ α × σ β ) ·r -vanish. However, the potentials can induce transitions between different states similar to an electric dipole moment. While not the focus of this section, atomic and molecular EDM experiments have been used to constrain spin-dependent, P -and T -violating potentials induced by axion exchange in Ref. [33]. In the context of the neutrino-mediated force, Ref. [57] has suggested probing atomic parity violation by measuring the optical rotation of light as it passes through vaporised atoms.
The expectation value of the SM-induced potential V LL αβ (r) can now be written as Recall that the functions I giving the average for the potential Computing the average of the potential V LR αβ (r) for Dirac neutrinos in Eq. (63) requires evaluating the average of the factor m i m j I LR ij (r)/r 3 . For distances r 1/(2m i ), which illustrates that the potential is too suppressed to be a useful probe of the non-standard coefficients c LR and c RR . In the Majorana case the potential V LR αβ (r) has the same r-dependence as V LL αβ (r) and so is added on top of the SM contribution as in Eq. (59).
We see in all cases that the expected shift from V LL αβ (r) is much smaller than this experimenttheory discrepancy. We see that the leptonic systems provide larger shifts in relation to the experiment-theory difference compared to the semi-leptonic systems. This is mainly due to the cut-off r c = 1/m Z being two orders of magnitude smaller than the charge radii of the proton and deuteron. Of the leptonic systems we see that the experimental measurements of the muonium splittings are the most precise -the predicted shift due to neutrino-exchange δE n−hfs αβ ≈ −150 mHz is around three orders of magnitude smaller than the experiment-theory difference. The hyperfine splitting of muonium is therefore the most stringent probe. H (e, p) / D (e, D + ) c LR e < 1.5 · 10 9 , c LR p < 5.5 · 10 9 g XY e · (g XY D − 6.48g XY p ) < 1.6 · 10 10 TABLE III. Upper limits on the non-standard coefficients c RL and c RR probed by the right-handed current potential V LR αβ (r) and coefficients g XY (for X, Y = L, R) probed by the scalar-scalar potential V SS αβ (r). c LR and c RR are constrained from the hyperfine splittings of the systems (f α , f β ), while the g XY are constrained from the 1S − 2S splittings. To avoid a helicity-suppression we assume three light active Majorana neutrinos with m 1 = 0.1 eV and NO masses and mixings. For simplification we take found from Eq. (122) to be which depends linearly on the coefficient c LR ββ . This potential relies on two effective interactions -one from SM CC and NC interactions and the other from a non-standard interaction -which may possess different cut-offs r c and r c . The cut-off appearing in Eq. (129) must therefore be the larger of these two scales. For simplicity we assume that the new physics arises around the EW scale m Z and therefore r c ≈ r c regardless of the exotic coupling strength g . This allows us to rewrite Eq. (129) as We now use Eq. (130) to compute the predicted shift as a function of the non-standard coefficient c LR . To simplify the sum over mass eigenstates (i, j) we take the coefficients to be diagonal in the mass basis, i.e. c LR α;ij = c LR α δ ij . We now write the inequality relating this predicted shift to the difference between experimental and theoretical values, and rearrange to put an upper bound on the value of c LR α . We note that c LR α;ij gets a contribution from the SM for Majorana neutrinos -however, even if we include this contribution it is too small to affect the upper bound derived for the non-standard coefficient. In Table III  We now consider the scalar-scalar potential V SS αβ (r) in Eq. (76) which does not depend on the the external particle spins -we must instead use the 1S − 2S splitting to derive upper bounds on the coefficients g XY . For this potential this splitting is found to be Taking again the differences in the experimental and theoretical values for the splittings, we derive the upper bounds on the coefficients in Table III . We now see that the constraints from positronium and muonium are roughly comparable while those from hydrogen/deuterium remain less stringent.

VI. NEUTRINO ELECTROMAGNETIC PROPERTIES
In this final section we will derive long-range potentials induced by possible non-standard electromagnetic properties of the neutrinos. The long-range potential induced by a neutrino magnetic dipole moment has been studied before, for example in Ref. [51].
Interactions between the neutrino mass-eigenstate fields and the electromagnetic field can be written generically as where Λ µ ij is a 4×4 matrix in spinor space which may contain space-time derivatives. To calculate an amplitude for the ννγ vertex one must take the matrix element of the neutrino current j µ eff (x) between initial and final neutrino states, where the vertex function Γ µ ij depends only on the momentum-transfer q. It is parametrised as The functions f Q ij (q 2 ), f A ij (q 2 ), f M ij (q 2 ) and f E ij (q 2 ) are the real charge, anapole moment, magnetic and electric dipole moment form factors, respectively. When coupling to a real photon with = µ ij and f E ij (0) = ij are the neutrino millicharge, anapole moment, magnetic and electric dipole moments, respectively. The above discussion is valid for Dirac neutrinos -for Dirac antineutrinos the form factors becomef For Majorana neutrinos we remember that the same electromagnetic process is described by two terms in the Lagrangian, and therefore the matrix element becomes This enforces the constraints on the form factors: The diagonal elements of the real charge, magnetic and electric dipole form factors therefore vanish for Majorana neutrinos -only the anapole moment form factor retains non-zero diagonal elements. All off-diagonal elements, or transition moments, can be non-zero depending on the relative CP phases η i of the neutrino mass eigenstates. If η i = η j then the off-diagonal elements of the real charge and magnetic dipole moment vanish, whereas if η i = −η j the offdiagonal elements of the anapole moment and electric dipole moment vanish.
In the low-energy effective field theory discussed in Sec. II one can generate electromagnetic properties for both Dirac and Majorana neutrinos. For Dirac neutrinos the following operator arises at dimension-six, which is written in the flavour-basis and where F µν is the electromagnetic field strength tensor.
If they are instead Majorana, ν ρR = ν c ρL and Eq. (138) must instead arise at dimension-seven, because the two ν L contained within SU (2) L doublets must contract with Higgs doublets.
It is nonetheless possible for a magnetic moment term to arise at lower dimension when introducing the right-handed Majorana states N R to the SMEFT. This is the LNV operator at dimension-five in Eq. (17), which is written in the EW basis, where ζ is an n × n matrix with n the number of sterile states.
We can rotate Eqs. (138) and (139) to the mass basis in a similar way to Eq. (3), where we re-iterate that ν = (ν 1 , ν 2 , ν 3 ) and n = (ν 1 , ν 2 , ν 3 , N 1 , N 2 , ...). In the mass basis, µ D and µ M are 3 × 3 matrices given by and ζ is a (3 + n) × (3 + n) matrix given by Splitting these in general complex dipole moments into real and imaginary parts as we obtain  139) and (141) we are again considering the type-I seesaw with N = 3 + n massive neutrinos, we can split for example the magnetic dipole moment into where the antisymmetric matrices µ M ν and µ M N contain the transition dipole moments for the light and heavy neutrino mass eigenstates respectively, while µ M νN contains the transition dipole moments between active and sterile states. From the form of the mixing matrixŨ in Eq. (20) in the seesaw limit we can write L ζ explicitly in the mass basis as It is clear from this that µ M ν and µ M νN are suppressed by the factors ε 2 and ε compared to µ M N , respectively.
We now move on to consider the long-range potential for the processes shown in the left two Feynman diagrams of Fig. 5. In these diagrams a pair of mass eigenstate neutrinos interacts via a SM CC and NC process at one vertex and via a photon at the other, coupled to the neutrino magnetic or electric dipole moment. The amplitude for this process is given by where the neutrino loop factor N µν ij is given by and the product of external fermion bilinears is H αβ µν = [γ µ P X ] α [γ ν ] β , where the X = L, R depends on the presence of a SM CC or NC interaction.
Taking the discontinuity of Eq. (149) and using Eq. (37), we obtain in the Dirac case, where we have normalised by the Bohr magneton µ B = e/2m e . Here, which take into account that the SM current can be at the interaction vertex of fermion f α and magnetic (or electric) dipole moment the interaction vertex of f β , or vice versa. The 1/(1 + δ αβ ) factor again takes into account double counting if α = β.
In the Majorana case we have where It can be seen that the magnetic moment does not contribute to the potential in the Majorana case -this is simply a case of the whole amplitude vanishing when only the magnetic moment and axial part of P L contribute. The first thing to observe in these potentials is that the rdependence, 1/r 3 , is the same as for the right-handed current potential V LR αβ (r) in the Dirac case. However, there is now a factor of αG F instead of G 2 F and the potential is now proportional to one power of the neutrino masses instead of two. For µ ij ∼ µ B and noting that G 2 F m 2 ν αG F m ν /m e for m ν ∼ 0.1 eV, we see that the potential is far less suppressed than V LR αβ (r) in the Dirac case. We can instead consider the process depicted by the Feynman diagram to the right of Fig. 5, where the two mediated neutrinos are coupled to both the external fermions by their magnetic or electric dipole moment. The dipole-dipole potential obtained in this case (valid for both Dirac and Majorana neutrinos) is where X γγ ij = µ . We see that there are two terms -one for the presence of two magnetic dipole moments and the other for two electric dipole moments.
The cross-term for a magnetic and electric dipole moment vanishes.
In Fig. 6 we compare (for positronium) the spin-independent potentials V V γ αβ (r) and V γγ αβ (r) to the spin-independent part of SM potential V LL αβ (r). We take a non-zero value of the magnetic moment, µ ij = µ ν = 10 −12 µ B , and let the electric dipole moment vanish. We see, as expected, that the potentials scale as 1/r 3 in the short-range limit r 1/(2m i ). However, unlike the potentials V LR αβ (r) and V SS αβ (r) they dominate over the SM potential for a wide range of distances. Because the vector-dipole potential V V γ αβ (r) is proportional to the neutrino masses and therefore is suppressed in the short-range limit, we focus instead on the shifts induced by the dipoledipole potential V γγ αβ (r). Using the same procedure outlined as in Sec. V to calculate the expectation value of the potential, we find We also derive upper limits on the heavy sterile neutrino magnetic moments µ ij = µ N (1 − δ ij ) for (i, j) = 4, 5, i.e. introducing two heavy Majorana neutrinos in the type-I seesaw. Active neutrino magnetic moments µ ν ∝ ε 2 and active-sterile transition magnetic moments µ νN ∝ ε are neglected.
where X γγ νν = In the second column of Table IV we show the upper bounds on the magnetic moments when we assume µ ij = µ ν such that X γγ νν = 9µ 2 ν (for three light Dirac neutrinos) derived from the positronium, muonium and difference in the deuterium and hydrogen 1S − 2S splittings. In brackets is the upper bound when we assume there to be three light Majorana neutrinos, for which µ ij = µ ν (1 − δ ij ) and X γγ νν = 6µ 2 ν . These limits also apply for the electric dipole moments. In the third column of Table IV we consider the scenario where two heavy sterile Majorana neutrinos are introduced in the type-I seesaw. Here the magnetic moments for the light active neutrinos and the active-sterile transition magnetic dipole moments are suppressed as ε 2 and ε respectively. We therefore take µ ij ≈ 0 for all (i, j) = 1 . . . 5 apart from µ ij = µ N (1 − δ ij ), where µ N is the transition dipole moment between the two sterile states. We again use the 1S − 2S splittings of the different systems to put an upper bound on this parameter, shown in Table IV.

VII. CONCLUSIONS
The exchange of force carriers are the basis of our understanding how structures are formed in nature. On scales larger than nuclei, the SM of particle physics incorporates the photon being largely responsible at the scales of atoms and larger. Beyond the SM, searches for massless or very light exotic mediators are being carried out in a large number of experiments and at different length scales, ranging from nuclear and atomic precision spectroscopy to the effects of fifths forces in astrophysics. Still within the SM, we already have more light particles, neutrinos, that are already known to be lighter than 0.1 eV but as fermions they cannot act as single exchange particles between matter particles. It is nevertheless possible that two neutrinos are simultaneously exchanged between two matter particles. In a Feynman diagram representation, this corresponds to a contribution at the first-loop order, see e.g. Fig. 2 for those arising in the SM. Due to the weakness of neutrino interactions and the loop suppression, the effect is small but nevertheless it may be possible to probe the effect of such a SM neutrino exchange in (exotic) atoms, especially muonium [55].
We have here considered an EFT approach parametrising the effect of potential New Physics to analyse the long-range potential induced by the exchange of two neutrinos. This includes all possibilities for the relevant four-fermion contact interactions between two neutrinos and two charged leptons or quarks. We have calculated both the spin-independent and spin-dependent long-range potentials as well as for a neutrino magnetic moment. Using our results, we discuss the potential of probing the exchange potentials and the underlying effective operator couplings using state-of-the-art atomic and nuclear spectroscopy experiments high precision QED calculations. Normalising the operator coefficients G eff of the relevant four-fermion contact interactions relative to the standard Fermi constant G eff = G F c XY , we have found that the current precision in atomic spectroscopy is sensitive to coefficients as low as c XY = O(10 2 ) for muonium. If the exchange is accompanied by a neutrino magnetic or electric dipole moment µ ν , values of order µ ν = O(10 −2 ) µ B are being probed.
We have worked in the low energy effective field theory approach to model both the exchange of SM EW bosons as well as any exotic contributions. In Ref. [53], the importance of the second order EW effects have been discussed which is usually ignored in the effective (Fermi theory) approach. Therein, the 1S hyperfine splitting energy shift of muonium was calculated in momentum space and in a gauge invariant fashion, including all relevant EW loop contributions, also those arising from electrons in the loop. For the SM case, they determine the same overall energy shift as in the present work as well as in [55], and we therefore do not consider any other ultraviolet-complete scenarios. The limits on effective neutrino operators extracted from atomic spectroscopy can be used to constrain New Physics scales. For example, the muonium 1S hyperfine splitting energy shift in the SM is of the order |δE 1-hfs | ≈ 0.14 α 4 m 3 e G F ≈ 6 × 10 −16 eV ≈ 150 mHz, and New Physics scales close to the EW scale are currently being probed. Future advancements in experimental muonium spectroscopy [53] and QED precision calculations [180,181] are expected to improve the sensitivity to |δE 1-hfs | ≈ 10 Hz ≈ 5 × 10 −15 eV. 3 While this will not improve on the existing limits from other processes as discussed in Sec. II C, atomic scale probes have the advantage that the effective operator treatment is valid down to very low energy scales corresponding to the Bohr radius, Λ NP αm e ≈ 3 keV.
g LL = S + P = 1 2 C S − iD S − C P + iD P , g RL = S − P = 1 2 C S + iD S + C P + iD P , g LR =˜ S +˜ P = 1 2 C S − iD S + C P − iD P , where flavour indices have been suppressed. We note that the flavour indices for fermions and neutrinos are swapped in our convention, i.e. L αβγδ = c LL γδ;αβ .
Appendix B: Spinor identities and non-relativistic limit A crucial step to take in deriving the spectral functions or absorptive parts of the invariant scattering amplitudes M (α,β) is taking the non-relativistic limit of the external interacting fermion bilinears. This can be done by expanding the bilinears to first order in both the 3momentum transfer q = p α − p α = p β − p β and the sum of 3-momenta P = p α + p β = p α + p β , where Γ a = {I, γ 5 , γ µ , γ µ γ 5 , σ µν } is one of the 16 irreducible products of γ matrices and u sα (p α ) and ξ sα are respectively the 4-component Dirac spinor and 2-component Weyl spinor for a fermion f α with mass m fα , 3-momentum p α and spin s α .
This expansion must be made for the external fermion bilinear at each of the interaction vertices. Hence the bilinears only appear as the products [Γ a ] α [Γ b ] β . We retain the higher orders terms in P and q arising from this product for comparison with the basis of 16 operators in Ref. [128], a complete set of scalar operators constructed from two spins and two momenta.
The products of scalar-like fermion bilinears are which are proportional to the O 1 , O 3 and O 9 ±O 10 operators in Ref. [128] respectively. Throughout this work however we consider a SM weak vector interaction at one vertex and an arbitrary scalar, vector or tensor-like interaction at the other. These fermion bilinears are therefore not used in this work, but are relevant for axion-mediated long-range potentials [33,38,182].
The products of vector-like fermion bilinears are The first of these products contains terms proportional to The relevant products of scalar-like and vector-like (which must be contracted with the momentum exchange q µ ) fermion bilinears are proportional to O 2 and O 3 respectively. These are used for the case of a scalar interaction at one vertex and a CC or NC interaction at the other.
A second set of dimensionless integrals J N ij (r) is derived from the above by performing the