All orders factorization and the Coulomb problem

In the limit of large nuclear charge, $Z\gg 1$, or small lepton velocity, $\beta \ll 1$, Coulomb corrections to nuclear beta decay and related processes are enhanced as $Z\alpha/\beta$ and become large or even non-perturbative (with $\alpha$ the QED fine structure constant). We provide a constructive demonstration of factorization to all orders in perturbation theory for these processes and compute the all-orders hard and soft functions appearing in the factorization formula. We clarify the relationship between effective field theory amplitudes and historical treatments of beta decay in terms of a Fermi function.

Factorization theorems underlie much of our ability to retain theoretical control in precision measurements involving nucleons, nuclei, and other hadrons [52][53][54][55][56][57].Factorization arises from the separation of different energy scales involved in a physical process, with the components in the factorization formula identified with contributions from each scale [58,59].In terms of a sequence of effective field theories (EFTs), the components are identified as the corresponding sequence of matching coefficients, and the final low-energy matrix element.Historically, Coulomb corrections have been understood not in terms of EFT, but by appealing to wavefunction methods i.e., solutions of the Dirac or Schrodinger equation [2,[60][61][62].Such wavefunction descriptions contain the correct long-distance behavior, which is however, intertwined with model-dependent short-distance behavior.Separating scales allows us to systematically resum logarithms and study higher order radiative corrections using standard tools of effective field theory.The interplay of high order Coulomb corrections with other subleading effects is crucial for precision measurements, and in particular for nuclear beta decays [50].
In this paper, we demonstrate factorization for radiative corrections induced by photon exchange between charged leptons and a static Coulomb field, and compute explicit all-orders expressions for the components of the factorization formula.We describe how traditional wavefunction methods are related to dimensionally regulated Feynman integrals order by order in perturbation theory.Using this correspondence, and a new all-orders calculation of the short-distance region, we extract the universal MS Coulomb corrections to the matrix element for a contact interaction (as is relevant for nuclear beta decays) to all orders in perturbation theory.
The remainder of the paper is organized as follows.Section 2 introduces notation for Coulomb corrections from a diagrammatic perspective.Section 3 considers the Schrodinger-Coulomb problem and establishes the correspondence between wavefunctions and the diagrammatic expansions.Section 4 considers the Dirac-Coulomb problem and extracts the relevant EFT matrix element to all orders in perturbation theory.Section 5 highlights new and interesting features of the preceding analysis and comments on phenomenological applications.

Coulomb corrections and contact interactions
Consider a reaction that takes place via an effective contact interaction in the vicinity of a heavy particle with charge Z.The outgoing charged particles ("leptons") can exchange photons with the heavy particle ("nucleus").For neutral current processes (i.e., when the initial and final nuclear states have the same charge Z), QED radiative corrections can be straightforwardly organized as a series in powers of α, Zα and Z2 α, with each power being separately QED gauge invariant. 1We will consider the static limit, in which the particle of charge Z in the initial and final state is heavy, so that recoil corrections can be neglected.For low momentum probes satisfying |p| ≪ 1/R with R the charge radius of the heavy particle, the point-like limit is applicable and universal corrections to the amplitude can be computed using Feynman rules for a static external Coulomb field [64].In this static limit, terms ∼ Z m α n vanish for m > n [65].In the following we consider the leading series of terms, ∼ (Zα) n , for n ≥ 0.
As an explicit example, consider di-lepton production via a short-range neutral current in some nuclear decay: where states A and B have charge Z, and v µ B = v µ A = v µ = (1, 0) which defines the static limit.As discussed above, this can be reduced to an external field problem describing the production of a di-lepton pair in a static Coulomb field, For a contact interaction, the Coulomb corrections on each leg factorize.For example, In general, the Coulomb corrected matrix element can be represented as We may therefore, without loss of generality, study the Coulomb corrections on a single leptonic leg: any process with multiple leptons that is mediated by a hard current reduces to a product of Coulomb corrections on each leg.
The perturbative series for a relativistic lepton contains non-trivial Dirac structure that must be inserted between external polarization spinors.For the process (1), Eq. ( 4) becomes Traditional analyses of beta decay are expressed in terms of position-space Coulomb wavefunctions for the leptons evaluated at the origin of coordinates, ψ(r = 0) [2,[60][61][62].As we discuss in Appendices B and C, the diagrammatic series represented in Eq. ( 4) is equivalent to a wavefunction solution.However, starting at twoloop order the wavefunction ψ(r = 0) is UV divergent. 3The amplitude must be renormalized and matched consistently with the underlying contact interaction.In order to execute this program, we phrase the problem in terms of factorization of momentum space amplitudes, using dimensional regularization in the MS scheme.Coulomb corrected amplitudes can then be matched consistently to underlying quark-level Lagrangians, and model-dependent position-space wavefunctions are replaced by a systematic expansion in EFT operators.

Schrodinger-Coulomb problem
Consider the quantum mechanical corrections to a tree-level process for a final-state particle of mass m and electric charge (−e) scattering from a Coulomb potential with source charge (+Ze): (we suppress the overall tree-level amplitude factor) Here, D = 3 − 2ϵ is the spatial dimension with dimensional regularization parameter ϵ, and λ is a photon mass regulator.The Schrodinger-Coulomb problem describes the limit p ≪ Λ UV ≪ m, where Λ UV ∼ R −1 denotes the scale of nuclear or hadronic structure.
The amplitude (6) can be evaluated at each order in perturbation theory.With the photon mass regulator in place, the integrals in Eq. ( 6) are UV and IR finite at ϵ → 0. By convention M (0) = 1 and at one-loop order, where the final expression denotes the limit ϵ → 0, λ → 0.

Factorization
Two momentum regions [58,59] are relevant in the integrals (6): the soft region with |L| ∼ λ; and the hard region with |L| ∼ p. Neglecting power corrections in λ/p, the amplitude may be written In the language of effective operators, M H represents a matching coefficient and M S represents a low-energy operator matrix element, when the full theory represented by Eq. ( 6) is matched onto a low-energy theory containing only soft degrees of freedom. 4

Soft factor
The soft limit of Eq. ( 6) is readily seen to exponentiate, yielding the soft factor to all orders [68,69], where the one-loop result is

Hard factor
The hard factor can similarly be evaluated explicitly, order by order in perturbation theory.The hard momentum region is isolated by expanding at L 2 ≫ λ 2 .At first order, M where5 β = p/m and the MS coupling α is related to the bare charge e in D = 3 − 2ϵ dimensions as6 At ϵ → 0, it is readily seen that H .
At second order M where the integral is evaluated in Appendix A. At third order, M The evaluation of this integral is also performed in Appendix A. At higher-loop order, direct evaluation of integrals becomes increasingly difficult.We will see how wavefunction methods provide a closed-form expression for arbitrary loop order.

Renormalization
Before turning to the all-orders discussion, we present the renormalized hard matching coefficient through three-loop order in the MS scheme.Identifying the above amplitudes as bare matching coefficients, M H ≡ M bare H , writing and requiring that Z(µ) has only 1/ϵ terms when expressed in terms of α, we find with The renormalized matching coefficient (at ϵ = 0) is then where α reduces to the on-shell QED coupling α at ϵ → 0 (recall that there are no dynamical leptons in the nonrelativistic theory).Since the product M S M H is UV and IR finite (at λ ̸ = 0), the quantity Z is identical to the (MS) operator renormalization constant for the soft operator, From the explicit form of Eqs. ( 9) and ( 10), the renormalization constant to all orders is given by in agreement through three-loop order with Eq. ( 18).The renormalized soft function is

Wavefunction solution and all-orders hard function
We recognize Eq. ( 6) as the perturbative expansion of the position-space wavefunction evaluated at r = 0 for a particle scattered by a Coulomb source and described by the Hamiltonian, The all-orders solution at leading power is (see Appendix B) , where ψ (−) denotes the scattering solution that matches asymptotically to a plane wave plus an ingoing spherical wave.Combining Eqs. ( 22) and ( 24) we obtain the closed form result This result reproduces the above results, cf.Eq. ( 19), through three-loop order.

Dirac-Coulomb problem
In place of Eq. ( 6), consider the amplitudes for a relativistic fermion in the Coulomb field of an extended object with a charge form factor The Dirac-Coulomb problem corresponds to the hierarchy p ∼ m ≪ Λ UV .For F (L 2 ) = 1, E = m, and / p− / L i +m → 2m, the amplitude reduces to the Schrodinger Coulomb problem (6).The fermionic case represented by Eq. ( 26) involves nontrivial Dirac structure, and a dependence on UV momentum scales |L| ≫ p.In the limit of a point-like source we have F (L 2 ) = 1.Similar to the Schrodinger-Coulomb case, we first consider the low-order contributions.At one-loop, for λ → 0 and ϵ → 0, Similar to Eq. ( 8), we can express the result, up to λ/E power corrections as the product of soft and hard factors, with M S as in Eq. ( 10), and M H now containing two different Dirac structures, At tree level, the hard factor is given by and at one loop, M At two loop order, using integrals from Appendix A, M

Factorization
The integrals in Eq. ( 26) are UV divergent by power counting when F (L 2 ) = 1.The explicit computations above show that M S M H is UV divergent beginning at two loop order, indicating sensitivity to short distance physics.
Regulating UV divergences with F (L 2 ) introduces a new UV scale, and a corresponding momentum region in loop diagrams with |L| ∼ Λ UV ≫ p.The factorization formula is In the following, we compute the explicit form of M UV using an illustrative charge form factor.We then introduce an alternative finite-distance regulator that permits an all-orders solution of M UV .Combined with an all-orders solution for the total amplitude M using the same UV regulator, and the all-orders solution of M S , we then extract M H to all orders in perturbation theory.

UV contribution from a charge form factor
In dimensional regularization, the factor M UV is computed by setting λ = p = 0.For simplicity, we take At one loop order, M Nontrivial contributions begin at two-loop order, M Let us compute renormalized expressions through two-loop order.In the MS scheme, the renormalized soft function is again given by Eq. ( 22), The renormalized hard function through two loop order is and the renormalized UV function for the form factor in Eq. ( 33) is It is readily checked that with the explicit results ( 36), ( 37) and ( 38), the product (32) is independent of µ S and µ H through two loop order.

UV contribution with finite distance regulator
Consider the series of amplitudes representing the perturbative expansion of the Dirac wavefunction at finite distance: For loop momentum |L| ≫ 1/|r| the rapid oscillations of the exponential regulate the integral, and the finite distance r acts as UV regulator.In the limit 1/r ≫ p, the amplitudes are described by the factorization theorem Eq. ( 32).The finite distance regulator is convenient since regulated amplitudes correspond to coordinate space solutions of the Dirac equation, which for |p| ≪ 1/r have a closed form solution (cf.Appendix C).We may relate the finite distance regulator scheme to a conventional MS-regulated amplitude by applying the method of regions [58,59].The finite-distance regulated amplitude M r satisfies the factorization theorem (32), where the UV matching coefficient depends on r.
We will now show that M UV (r) can be computed to all orders in perturbation theory.This fact is related to the structure of the loop integrals with a finite distance regulator, Eq. ( 39), as compared to a charge form factor, Eq. ( 26): the regulator affects only the final ( d D L n ) loop integration, so that all of the preceding integrals are recursively one-loop.Details are presented in Appendix D, with the results for bare amplitudes at arbitrary even and odd orders in perturbation theory respectively: The quantity α is given in terms of the MS coupling α in Eq. ( 12), by As discussed in Appendix D, both series can be expressed in closed form for arbitrary nonzero ϵ in terms of Bessel functions.The MS renormalization constant can also be computed in closed form.A careful treatment of the small-ϵ asymptotics of the bare amplitudes, cf.Appendix D, then yields the all-orders result, where η = 1 − (Zα) 2 .The result (44) is renormalized in the MS scheme using the coupling defined in Eq. ( 12).

Wavefunction solution and all-orders hard function
The amplitude ( 39) is related to the perturbative expansion of a solution to the Dirac equation, where ψ (−) (r) denotes the solution that is asymptotically a plane wave plus incoming spherical wave.The solution, ignoring power corrections in λ/p and p/r −1 , is Here F (Z, E, r) is the Fermi function, the phase factor e iϕ is given by and the quantity M is given by In the Dirac (i.e., "Bjorken and Drell") basis for γ µ and with relativistic normalization u(p) † u(p) = 2E, the expression is where and χ is a two-component spinor.Using Eq. ( 46), the explicit all-orders results for M S in Eq. ( 22), and M UV in Eq. ( 44), the hard function appearing in the factorization formula (32) is The amplitude has been explicitly decomposed into separate factors depending on a single scale, λ, p, or r −1 (here we are not distinguishing the scales p, m, and E).We remark that the explicit appearance of exp(γ E ) accompanying 2p/µ in Eq. ( 53) may seem unexpected since the hard amplitude must match conventional MS renormalized amplitudes order by order in perturbation theory.However, these factors cancel against implicit factors 7 of γ E from the expansion of Γ(η − iξ)/Γ(1 + 2η).Given Eq. ( 53) we can extract the anomalous dimension for contact operators to all orders in Zα.We differentiate M H with respect to µ H and obtain This is the contribution to the anomalous dimension from each light-particle leg.For example the operator mediating Eq. ( 1) has an anomalous dimension of 2γ O .For an operator mediating a beta decay, 54) is the leading-Z contribution to the anomalous dimension [50,51].Including the one-loop beta function with n f dynamical fermions, the scale dependence of contact operators can be obtained in closed form where we have introduced the notation η L,H = η(α L,H ).This expression is useful when analyzing QED radiative corrections for the beta decays of heavy nuclei [50].
The hard function (53), describes the limit p ∼ m ≪ Λ UV , where Λ UV ∼ R −1 denotes the scale of nuclear or hadronic structure.When the lepton is non-relativistic, p ≪ m ≪ Λ UV , it is convenient to expand the hard function as where Allowing for arbitrary values of ξ, we find through second order in β, where ψ denotes the digamma function, ψ(x) = Γ ′ (x)/Γ(x).At each order in β 2 , the expressions (57) sum an infinite series of terms involving powers ξ n .At β → 0, the leading term for the "large" upper component M + H reduces to the Schrodinger-Coulomb result (25) which corresponds to a "non-relativistic Fermi function", cf.Refs.[49,60].

Discussion
The formula (32), and its non-relativistic analog (8), provides an all-orders explicit demonstration of factorization for the Coulomb problem.We find that Coulomb corrections factorize among different legs for a contact interaction (see Section 2).The universal hard matching coefficient in this formula, M H in Eq. ( 53), can be applied to different processes, and large logarithms can be summed to all orders using renormalization group methods.The non-relativistic limit for p ≪ m ≪ Λ UV is given by Eq. (57).By identifying the amplitudes as quantum field theory objects in a standard regularization scheme (i.e., MS scheme in dimensional regularization), we can systematically compute subleading perturbative contributions and match to hadronic and nuclear matrix elements.More detailed discussions of these points are presented elsewhere [9,50,51].It is interesting to note that for unpolarized observables to beta decay, the spin-summed matrix element squared,8 differs from the historically defined Fermi function even when evaluated at r −1 H = µ H e γE .We observe that finitedistance regulated amplitudes have special algebraic properties that allow for explicit all-orders expressions, for both bare and renormalized matrix elements as shown explicitly in Eqs. ( 41), ( 42) and (44).This example of all-orders renormalization may be of formal interest.
As an illustration of how the formalism applies to different processes, let us return to Eq. ( 1).For definiteness, suppose that the neutral current reaction is mediated by exchange of a vector boson of mass m B .The tree-level amplitude depicted in Eq. ( 2) takes the form where Γ tree = γ 0 (A + Bγ 5 ) for some numbers A and B. When Λ ≫ m B ≫ p, the boson mass plays the role of UV regulator. 9The factorization formula describing the infinite sum in Eq. ( 2), is Here the conjugate amplitude is denoted M = γ 0 M † γ 0 .It is straightforward to compute M UV from the diagrams in Eq. ( 2), neglecting charged lepton masses and momenta.Through two-loop order, after MS renormalization, It is readily seen that the scale dependence of M UV (µ H ) cancels against the product of M H (µ H ) for the charged leptons.
An important application of the formalism presented above is to the description of precision nuclear beta decay, e.g., for |V ud | determination [36] and tests of first row CKM unitarity [23, 24, 26, 36, 40-44, 70, 71].Consider the decay of a heavy atom to a negatively charged ion, a positron, and a neutrino [22,25,30,[72][73][74][75][76][77][78][79][80][81][82][83][84][85][86][87], Beta decays are a complicated multi-scale problem, involving energies from the weak scale ∼ 100 GeV, down to scales set by atomic screening ∼ 100 eV.Structure dependent corrections e.g., due to nuclear charge distributions, can be subsumed into a short-distance Wilson coefficient in the point-like theory considered here.The embedding of Coulomb corrections in a broader EFT framework is crucial for the systematic separation of physical scales and computation of QED radiative corrections.For charged current processes such as beta decays, the chargemismatch between the initial and final heavy particle (i.e., nucleus) introduce sub-leading effects whose analysis can be substantially simplified using eikonal algebra [63].Systematic evaluation of these subleading corrections differ from previous phenomenological approaches and lead to numerical differences that are larger than the existing estimated error budget for outer radiative corrections [36,50]; detailed calculations are presented elsewhere [50,51]. the

A Loop integrals
We collect here some results for loop integrals that are used in the main text.Integrals are defined in Euclidean D-dimensional space with D = 3 − 2ϵ.

A.1.1 Scalar Integrals
Consider the two-loop integral Using that the integral of a total derivative vanishes in dimensional regularization, and inserting (∂/∂L i 2 )L i 2 and (∂/∂L i 2 )L i 1 under the integral, yields the following "integration by parts" [88] relation, where we use the shorthand m ± to denote the raising or lowering indices in J, e.g., 2 ± J(a 1 , a 2 , a 3 , a 4 , a 5 ) = J(a 1 , a 2 ± 1, a 3 , a 4 , a 5 ).In particular, for the two-loop integral appearing in Eq. ( 14), where the integrals on the right-hand side are recursively one-loop and are readily evaluated:

A.2 Three loop integrals
Consider the three-loop integral, Integration by parts identities are (cf.Eq. ( 64) at so that the integral of interest in Eq. ( 15) is The first integral in Eq. ( 75) is given by the product of two-and one-loop integrals, where J(0, 1, 1, 1, 1) is evaluated above.The second integral in Eq. ( 75) is recursively two-loop, To evaluate J 0, 2, 1, 1, 1 2 + ϵ , we first perform the L 3 integral in Eq. ( 77), where we introduce Integration by parts for K yields so that As a function of the integration variable x in Eq. ( 79), the terms on the right side of Eq. ( 82) are Each integral may be evaluated as a series in ϵ, yielding,

B Wavefunction solution: Schrodinger-Coulomb
Consider the Lippmann-Schwinger equation and its related Born series for the solution of the Schrodinger equation where Ĥ0 = p2 /(2m) is the free Hamiltonian and V = V (x) is the potential.For a finite range potential, the +i0 (−i0) prescription in Eq. ( 85) corresponds to a plane wave plus outgoing (incoming) spherical wave at large distance.Inserting a complete set of momentum eigenstates we arrive at where Ṽ (L) = d 3 x e iL•x V (x) is the potential in momentum space.In particular, for a Yukawa potential, V (x) = (−Ze 2 ) exp(−λ|x|)/(4π|x|), we have Ṽ (L) = −Ze 2 /(L 2 + λ 2 ).Setting x → 0 and choosing the outgoing +i0 prescription, the wavefunction ψ (+) p (0) provides an all-orders solution for the amplitude Eq. ( 6).
Let us solve the Schrodinger equation, in the limit where λ ≪ |p| (but to all orders in Zα).Here r = |x|.Let us write ψ p (x; λ) = e ip•x F p (x, λ).Choosing p along the ẑ direction, p = pẑ, we look for the solution that reduces to F = 1 at z → −∞ to obtain ψ (+) , and the solution that reduces to F = 1 at z → +∞ for ψ (−) .The differential equation for F is We may now apply boundary layer theory [89], solving for solutions at short and long distances and matching the solutions in their common domain of validity p −1 ≪ r ≪ λ −1 .For r ≪ λ −1 , the Schrodinger equation is with solution where 1 F 1 (a, b, c) is the confluent hypergeometric function.For r ≫ p −1 , the Schrodinger equation is with solution In the overlap region p −1 ≪ r ≪ λ −1 , the respective solutions can be expanded as Identifying in the overlap region, and using that 1 F 1 (a, b, 0) = 1, we have, up to λ/p power corrections, The incoming solution ψ

C Wavefunction solution: Dirac-Coulomb
The Dirac equation can be similarly shown to have a Lippmann-Schwinger solution and associated Born series.Let us define Φ(x) = u(p)e ip•x , where u(p) is a Dirac spinor.The solution of the Dirac equation with a potential can be written as The amplitude of interest, Eq. ( 39), is given by ū(p)M r = ψ (−) (−r) = [ψ (−) (−r)] † γ 0 .We require the solution ψ (−)   with a small but non-zero photon mass λ. References [90,91] present the angular momentum components for the strict λ = 0 solution, which is related to our problem by a normalization that must be computed.
To determine the complete solution including λ dependence, we identify this solution with ψ (−) < , up to a normalization that is fixed by matching to ψ (−) > in the overlapping region of validity.For simplicity we perform the matching by projecting onto the S-wave component of the outgoing spherical wave.
Let us consider the upper components of ψ (±) in the Dirac basis for γ µ , and introduce where χ is a 2-component spinor.Similar to the Schrodinger-Coulomb problem, we look for solutions F > when r ≫ p −1 , and F < when r ≪ λ −1 .The large-distance solution obeys an identical equation to the Schrodinger-Coulomb problem (with ξ = Zα/β and β = p/E representing the relativistic velocity).The solution for F (−) > is given in Appendix B, and for the matching we require the small-r limit.Considering the outgoing spherical wave component, the S-wave projection is The relevant component of the small-r solution involves the quantity [90,92] where exp(iκ) = (1 + imξ/E)/(η + iξ).From the large-r limit of this expression, taking the outgoing spherical wave component, we have Comparison of Eqs.(97) and (99) in the overlap region p −1 ≪ r ≪ λ −1 determines C(p, λ).Using 1 F 1 (a, b, 0) = 1, and taking the r → 0 limit of the complete solution, we have (cf.Eq. ( 16) of Ref. [92]) where

D All orders UV function with a finite-distance regulator
We can compute the UV matching coefficient introduced in Eq. ( 40) by setting λ = p = 0 and evaluating the remaining integrals using dimensional regularization.Examining the perturbative series we find that the (bare, unrenormalized) UV matrix element has the following structure where In particular, all even orders of perturbation theory contribute to F 1 and all odd orders to F 2 .The lowest order loop integrals are given by I and for higher orders, These integrals are recursively one-loop, and can be evaluated by repeated use of the following identity: where ν 1 = 2 and ν j+1 = ν j + 3 − D so that ν j = 2 + 2(j − 1)ϵ.The final integral involving e −iL2n•x for I In particular, when expressed in terms of α, the coefficients in the perturbative expansion of F bare i are expressible entirely as rational functions of ϵ.Choosing µ = (re γE ) −1 so that α can be identified with the MS coupling, we find that the MS operator renormalization constant can be written as exp 1 ϵ n a n gn for some numbers 10 a n .The sequence of coefficients can be related to the Catalan numbers C(n) = (2n)!/(n!(n + 1)!). 11The series in the exponent then converges, and is given by log The series in Eq. ( 108) can also be summed, and converges for any nonzero ϵ.The answer is given by Using Eqs. ( 113) and (114) we can see how renormalization works at all orders in the coupling.We require the ϵ → 0 asymptotic behavior of the Bessel functions.The relevant identity is (cf.Eq. ( 10.20.4) of Ref. [94]) where we adopt the notation of Ref. [94] and use ∼ to denote "asymptotic to".Using this identity, Sterling's approximation, and the large argument limit of the Airy function it is straightforward to show that For F 1 it is convenient to introduce 1/2ϵ ′ = 1 − 1/2ϵ and g′ = g(1 + 2ϵ ′ ).We then find Notice that the form of Z that we obtained from recognizing the infinite sequence using our perturbative result is precisely what is needed for all-orders renormalization, cf.Eqs. ( 112) and (116).We find (119) The µ dependence of the renormalized coefficient functions F i is governed by the anomalous dimension, and γ O is determined by the coefficient of 1/ϵ in the corresponding MS operator renormalization constant: (121) 10 The leading orders obtained by direct evaluation from Eq. ( 108) are We have checked explicitly to 16 th order in g that the renormalization constant can be written as exp 1 ϵ n ang n , consistent with the explicit all-orders expressions in Eqs. ( 116) and (117). 11We were able to identify this sequence with help from the Online Encyclopedia of Integer Sequences [93].

F 2 |
µ=(re γ E ) −1 = lim ϵ→0 Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.This work is supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DE-SC0011632 and by the Walter Burke Institute for Theoretical Physics.RP acknowledges support from the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DE-SC0011632 and the Neutrino Theory Network Program Grant under Award Number DE-AC02-07CHI11359.