Improved Constraints on the Variation of the Weak Scale from Big Bang Nucleosynthesis

We analyze the recent EMPRESS data on the 4 He abundance in Big Bang nucleosynthesis as a function of the Higgs vacuum expectation value v and compare our calculation to the recently published work of Burns et al. [1]. The EMPRESS result for the 4 He abundance can be explained within 2 σ by 0 . 03 ≤ δv/v ≤ 0 . 07. However, as first noted in [1], a Higgs VEV in this range would worsen the discrepancy between theory and experiment for the deuterium abundance significantly.


Introduction
Big Bang nucleosynthesis (BBN) is a fine laboratory for the possible variations of the fundamental constants of nature, see e.g.[2][3][4][5][6][7].In a recent publication [1], Burns et al. studied possible constraints on a variation of the Higgs vacuum expectation value (VEV) from BBN simulations.They use the new PRyMordial code [8] to compute the Helium-4 ( 4 He) and deuterium (d) abundances for different values of the Higgs VEV and compare their results to the newest observation by the EMPRESS collaboration [9] that differs from the values given by the Particle Data Group (PDG) [10] and from current theoretical predictions.The possible interval for the Higgs VEV that they found reproduces the EMPRESS 4 He abundance would, however, worsen the deviation for deuterium significantly.This work triggered the investigations here, which improve the methodology of earlier studies, such as [11][12][13][14][15][16][17], and thus gives more reliable bounds on the possible variations of the weak scale.Here, we mostly consider the element abundances from the PDG as the EMPRESS data are presently not included in the PDG listings.
Changing the Higgs VEV v influences all elementary particle masses as they scale linearly with v, assuming that the Yukawa couplings remain unchanged.This influences not only the neutronproton mass difference, which is particularly significant for the 4 He abundance, but also the QCD scale Λ QCD .Burns et al. take the latter into account [1].They also use some general scaling arguments for the Higgs VEV dependence of the pertinent nuclear interactions.While this is certainly a legitimate procedure, we believe that the calculation can be improved in some parts, as shown here.In particular, the treatment of the neutron-proton mass difference, see Sect. 2 is updated, and the quark mass dependence of the deuteron binding energy, see Sec. 3, will be the subject of our work.In the most sophisticaed treatments, this has been considered using chiral effective field theories either in a pionless [16] or pionfull [17] approach.However, there is a sizeable uncertainty related to the quark mass variation of the four-nucleon contact interactions, which have so far been estimated using naive dimensional analysis [18,19] and also the so-called resonance saturation approach [17,20].Here, we overcome these limitations by combining lattice QCD data with the low-energy theorems for nucleon-nucleon scattering [21,22] (and references therein).The weak n ↔ p rates are discussed in Sec. 4 and the n + p → d + γ reaction is considered in Sec. 5. Our results are displayed and discussed in Sec. 6. Appendix A contains some details on a one-boson-exchange model that can be used to study the deuteron for a varying Higgs VEV.In Appendix B, we compare in detail the differences between our work and previous publications on this topic.

Neutron-proton mass difference
The quark masses change linearly with the Higgs VEV and thus also the mass difference between up-and down-quark scales with the relative difference δv/v.Throughout this work, we assume that the corresponding Yukawa couplings are constant.The neutron-proton mass difference has a QCD part coming from the quark mass difference and a QED part coming from Coulomb effects between the charged quarks in the nucleon.To a good approximation, the latter should stay constant under variation of the Higgs VEV.The QCD contribution to the n − p mass difference is not, however, equal to just the current quark mass difference1 but has also contributions from the matrix element ⟨p|uū − d d|p⟩ [23].This can be understood easily from the two-flavor QCD Hamiltonian, where the explicit chiral symmetry breaking originates from the mass term m u ūu + m d dd, and the neutronproton mass difference is given by its isovector I = 1 component, In [24], the electromagnetic contribution to the neutron-proton mass difference was precisely determined using the Cottingham sum rule and the strong part can then be deduced from the physical mass splitting, leading to The neutron-proton mass difference as a function of the Higgs VEV should therefore be taken as We note in passing that the separation of the strong from the electromagnetic contribution is afflicted with some problems, for a pedagogical discussion see [25].
3 Pion mass dependence of the deuteron binding energy and the nucleon-nucleon scattering lengths Changing the quark masses influences, of course, also the pion mass.To first order in chiral perturbation theory the average pion mass can be found from the Gell-Mann-Oakes-Renner relation where the quark masses scale linearly with δv/v and the constant B 0 , which is related to the scalar singlet quark condensate, is taken to depend on δv like the QCD scale Λ QCD , which [1] found to be proportional to (1 + δv/v) 0.25 .The pion mass at a given value of the Higgs VEV is therefore Here, M phys π = 138.03MeV is the physical average pion mass for the present value of v (meaning δv = 0).For the Higgs VEV varying by ±10%, that is |δv/v| ≤ 0.1, the pion mass varies between 129 MeV and 147 MeV.A change in the pion mass affects also other physical quantities that are very relevant for BBN simulations -above all the average nucleon mass m N .The axial vector coupling g A as well as the pion decay constant (which will become relevant later) do also depend on M π .One could use chiral perturbation theory for calculating the pion mass dependence of these parameters, as was done, e.g., in [17].However, for already moderate changes in the pion mass, chiral perturbation theory becomes inaccurate for certain observables in the one-baryon sector.More precisely, most lattice data can be fit well with a linear function in M π , see e.g.[26], and the pion mass dependence of g A exhibits strong cancellations of one-and two-loop contributions [27].We therefore fit a rational function to lattice QCD data for these quantities at different values of M π .Quite differently, the pion mass dependence at one-loop of the pion decay constant matches the lattice data for not too large quark masses.Consider first the nucleon mass m N .We fit the data from Ref. [28] with a rational function where m phys N = 938.94MeV is the value of the physical nucleon mass, that is for δv = 0. We note, however, that in the range for the Higgs VEV variation discussed here, the third order chiral perturbation theory (ChPT) representation of m N gives a pion mass variation consistent with the one from fitting to the lattice data.Still, we prefer to work with this fit function as it allows for larger variations in M π and thus δv.Similarly, for the axial-vector coupling, the data from [29] are fitted with a polynomial Including only the more recent and more precise lattice QCD data from [30], one arrives instead at Figure 1: Rational function fit to lattice QCD data for the nucleon mass m N ( [28]) and the axial vector coupling g A ( [29,30]) as a function of the pion mass squared.We forced the fits to go through the physical point (black filled square).The bands correspond to 1σ-errors of the fit.The dashed line and hatched band display the fit for g A using only data from [29].
where g phys A = 1.267 denotes the physical value at δv = 0.In all these equations, M π is given in GeV.The fit results for m N (left panel) and g A (right panel) are displayed together with the used lattice data in Fig. 1.The bands correspond to 1σ-errors of the fit.Comparing both fits for the axial vector coupling, it is apparent that there is still quite some uncertainty when it comes to calculations on the lattice.We choose to work with the latter fit for g A that includes the more recent data.Choosing one parametrization or the other has, however, very little effect on the following calculations, because the range of pion mass variations is much smaller than what is shown in the figure.Note also that the one-loop O(p 3 ) ChPT result strongly deviates from the fit already in the small range of pion mass variations considered here.
The axial-vector coupling enters as a parameter when computing the weak neutron-proton interaction rates (see below).A change in the average nucleon mass has, of course, wide-ranging effects in the BBN simulation.We change the neutron and proton mass so that they fulfill both the new value for the mass difference and the average mass: In [1], the deuteron binding energy was computed using a one-boson exchange model taken from [31].More precisely, these authors used three Yukawa functions corresponding to the pion-, omega-and sigma-exchange and essentially tuned the σ-mass to get the corrrect deuteron binding energy at δv = 0.The pion mass was scaled with the varying Higgs VEV as given above and the other meson masses as (δv/v) 0.25 .It should be noted, however, that the quark mass dependence of the σ is more intricate than the one of a genuine quark-antiquark meson as worked out in detail in [32].In the Appendix, we discuss a somewhat more sophisticated one-boson-exchange model to better understand the approximations made in [1].Here, we use lattice QCD results for the dependence of the deuteron binding energy on the pion mass.In Ref. [33], five data points for the nn and deuteron binding energy or scattering length were compiled and the so-called low-energy theorems (LETs), see e.g., [21,22], were used to calculate the scattering length at a given pion mass from the binding energy obtained from lattice data at different quark masses without having information on the corresponding effective range.Note that at these pion masses, the nn system is also bound.Using this method and adding two recent lattice QCD data points from [34,35], we get the pion mass dependence of the nucleon-nucleon scattering lengths in the S -wave channels and the deuteron binding energy by fitting a cubic polynomial to the obtained data.The corresponding results can be found in Figs. 2 and 3.The scattering lengths become relevant later for the n + p → d + γ reaction.In more detail, the deuteron binding energy as a function of x = M π − M phys π is given by where all masses are given in units of MeV.We note that the situation about the deuteron being bound is still a disputed issue in lattice QCD, see e.g [36], thus we forced the fit to go to the physical value for the physical pion mass.Further, the inverse scattering length in the 1 S 0 channel is .9) while the inverse scattering length in the 3 S 1 channel (corresponding to the deuteron) is given by (3.10) Calculating the deuteron binding energy in this way for different values of the Higgs VEV yields a similar result to what was found in [1] as can be seen from comparing Fig. 4 to Fig.The data points from Refs.[34] and [35] are direct lattice QCD results.The fit was forced to go through the physical point.  of [1].It is unclear to us why these results are so similar.Using a more sophisticated one-bosonexchange model, we indeed find a different result as discussed in appendix A. 4 Weak n ↔ p rates Before the universe is cool enough for deuterium fusion to be efficient and for BBN to actually start, the weak interactions between neutrons and protons are dominant.Because especially the 4 He abundance is very sensitive to the number of neutrons existing at the beginning of BBN, it is quite important to accurately calculate the neutron-proton interaction rates.Here, one has to consider effects from a changing electron mass (as noted before, m e also scales linearly with δv) on the thermal background as well as the change in the W-boson mass that modifies the Fermi constant G F : The relevant Higgs VEV dependent part of the neutron lifetime is given by where f R is the integral over the Fermi energy spectrum.For small temperatures one can ignore temperature-dependent effects from the neutrino and electron Fermi distribution functions when calculating the neutron β-decay rate.So, in this approximation, f R is given by [43]  F(Z, p e ) is the Fermi function including Coulomb effects between the charged electron or positron and the proton in the final state.Of course, one could also include QED radiative corrections, finite-mass, electroweak magnetism and finite-temperature effects, all of which are considered in the calculation of the weak n ↔ p rates in PRyMordial [8].Because these are higher order effects, we will not include them in the Higgs VEV correction to the neutron lifetime.At high temperatures the Fermi functions of the neutrino and the electron cannot be ignored, thus one has to recompute the weak neutron-proton rates every time one changes the relevant parameters (see also [44] for a detailed discussion).Because the neutron-proton mass difference as well as the electron mass changes with δv, the variation in f R influences the neutron lifetime significantly.As mentioned above, the axial-vector coupling depends on the pion mass and therefore also on the Higgs VEV which should be included in the calculation.We therefore recalculate the n ↔ p weak rates for every value of δv/v changing the electron, neutron and proton mass as mentioned above as well as G F and g A .The neutron lifetime changes according to In Fig. 5 the fractional change in the properly normalized n ↔ p rates obtained from our calculation is displayed.Fig. 6 shows the neutron lifetime as a function of the Higgs VEV.We observe a significant variation from the experimental value at δv = 0, which is τ n = 878.4s [10].
5 The n + p → d + γ reaction Finally, we can use the pionless effective field theory (EFT) approach by [45] for calculating the n + p → d + γ reaction rate, which is the first nuclear reaction taking place in BBN and therefore quite important.In [1], it is assumed that all nuclear reaction rates scale like Λ QCD , so a factor of (1 + δv/v) 0.25 is multiplied to every rate.While this might be a reasonable approximation, one can make use of the analytic expression for the n + p → d + γ rate from [45] and plug in the nucleon mass, the 1 S 0 scattering length, the deuteron binding energy and the effective range ρ d at a given value of the Higgs VEV.The deuteron effective range can be found from the binding energy and Figure 7: Reaction rate for n + p → d + γ for δv/v = 0, ±0.1 calculated from pionless EFT [45] with nuclear parameters changed as described in the text.The shaded green band gives the original rate mutliplied by (1 ± 0.1) 0.25 , so the range in which the n+ p → d +γ rate varies in the calculation of [1].
the triplet scattering length a t via: The effective range of the 1 S 0 channel only comes in at higher orders and we can safely ignore changes in this parameter as well as possible scale effects in the isoscalar contributions to the cross section.The resulting reaction rate for δv/v = 0, ±0.1 is shown in Fig. 7. Changing the Higgs VEV therefore has a huge impact on the n + p → d + γ reaction.We also see that the functional form is not consistent with a simple scaling with (1 + δv/v) 0.25 .Before presenting our final results, let us summarize the differences in our calculation to the one in [1] and in [17].First, we use an updated value for the neutron-proton mass difference, as given by Eq. (2.3).Second, the quark mass dependence of the deuteron binding energy is obtained from lattice QCD data with the help of the nucleon-nucleon scattering LETs [21,22], whereas in [1] a one-boson exchange model was used with the appropriate scalings of the various meson exchanges.Interestingly, the resulting Higgs VEV dependence of the deuteron binding energy comes out very similar in both approaches.We also used lattice QCD data to model the pion-massdependence of the nucleon mass and the axial-vector current, which was not taken into account in [1].In [17], the pion mass dependence of these quantities was derived from chiral perturbation theory and the deuteron binding energy was found using chiral nuclear EFT at leading order, given by the one-pion-exchange potential together with the short-range four-nucleon contact interactions without derivatives, the latter exhibiting sizeable uncertainties.Then, we normalized the weak n ↔ p rates using the appropriately varied neutron lifetime including the change in the Fermiintegral, see Eqs. (4.2) and (4.3).The main difference to [17] is that only the linear dependence of the weak rates on the neutron lifetime was considered.This does not include temperaturedependent effects.Finally, we made use of the pionless EFT approach by [45] to include effects of a varying Higgs VEV on the n + p → d + γ rate, as explained in this section.

Results
Finally, putting everything together, our results can be found in Fig. 8. What is most noticeable is the behaviour of the deuterium abundance: it increases drastically for negative values of δv/v, which is a result of the large changes in the n + p → d + γ rate.We expect the variation of the 4 He abundance to be less pronounced than in [1] due to the different treatment of Q N , see Sect. 2.
Note that especially the deuterium abundance is sensitive to the choice of the rates for the 12 key reactions: in PRyMordial, one can choose between reaction rates collected by PRIMAT [46] or by the NACRE II collaboration.As the former was chosen in [1]2 , we made the same choice so as to make the results more comparable.Using another set of rates would result in a shift of the abundances and alter the constraints for δv/v accordingly.We did not use the improved rates from Ref. [44] since we are presently in the process of updating them and supplying theoretical uncertainties.This would go beyond the scope of this investigation.
The improved constraints for the Higgs VEV variation that can be derived by our calculation are collected in Tab. 1.We obtain these constraints by comparing our results to the observed values for the 4 He and d abundances given by the PDG [10].For 4 He we also present constraints coming from the recently published values by the EMPRESS collaboration [9], which was also done in [1].The drastic increase of the deuterium abundance for negative δv/v results in a very narrow range of possible variations of the Higgs VEV.This range is, however, included in the interval of possible δv/v that can explain the PDG 4 He abundance within 2σ.We find for the linear dependence of the abundances on the Higgs VEV We do note that the deuterium abundance can, of course, not be described by a linear function in δv, as is evident from Fig. 8.In Ref. [17] the values of 1.9±3.4 and −4.0±0.3 were found for deuterium and 4 He, respectively.This discrepancy comes mainly from the treatment of the n + p → d + γ rate and from including temperature-dependent effects in the n ↔ p rate.The dependences of the deuteron binding energy and the 1 S 0 scattering length are also very different.Note, however, that our treatment goes beyond the linear approximation and also, we have made clear improvements in the calculations of the Higgs VEV dependence of a number of pertinent quantities.Also, our bounds are dominated by the deuteron abundance and are more stringent than found before.In Fig. 9 we show the evolution of the abundances for different values of the Higgs VEV not including effects from the n + p → d + γ rate as would be implied by [45] but scaling it like all other rates proportional to the QCD scale Λ QCD , i.e, as ∼ (1+δv/v) 0.25 .The result for the deuterium Table 1: Upper panel: Constraints on the Higgs VEV variation δv/v derived from 2σ bounds of the observed values given in the PDG [10] and by the EMPRESS collaboration [9] for the abundances of 4 He and deuterium.Lower panel: Constraints obtained from a calculation where the n+p → d+γ rate is treated like in [1], see Fig. 9.   abundance is markedly different from the more realistic calculation and now the constraints for 4 He and d contradict each other.
For the constraints coming from the EMPRESS data, we can confirm the result of [1] that tun-ing the Higgs VEV so that the predicted 4 He abundance matches the new EMPRESS result would worsen the deviation between theory and experiment for the deuterium abundance significantly.In the future, we will use the methodology developed in Refs.[44,47] and here to upgrade the constraints on combined variations of the Higgs VEV and the electromagnetic fine-structure constant α from the BBN abundances of all light elements.Also, a fresh look on possible variations of the Yukawa couplings in the spirit of [48] would be an interesting venue.

A Deuteron binding energy from a one-boson-exchange potential
In Sect. 3 we used lattice QCD data to model the pion-mass dependence of the deuteron binding energy which was similar to the dependence found in [1] from calculating the binding energy using a one-boson exchange potential (OBEP).The authors of [1] use a superposition of three Yukawa functions corresponding to pion-, sigma-and omega-meson exchanges, where the first two are attractive and the last one supplies the required repulsion.The mass of the σ is varied so as to get the appropriate B d at the physical pion mass.Note that no meson-nucleon form factors are employed in that calculation.We try to confirm their results by implementing the OBEP potential described in [49], where also contributions from ρ-meson exchange are considered and a common form factor for all meson-nucleon vertices is utilizied.For more details, we refer to the textbook [50] (and references therein).In [49], the small contribution from the D-wave for the 3 S 1 -channel is neglected, that is one works strictly with angular momentum L = 0.The potentials corresponding to the different meson exchanges are then given in coordinate space as: and they are scaled with the Higgs VEV as explained in [1]: the pion mass with (δv/v) 1.25/2 and the other masses as (δv/v) 0.25 .The coupling constant g πNN is taken from the Goldberger-Treiman discrepancy [51] g πNN (δv) = g A (δv)m N (δv) where m N are g A given in Eqs.(3.3) and (3.5) as a function of the pion mass.F π (M π ) can also be found from a fit to lattice data by [52]: ) with M π given in GeV.Finally, we use the coupling constants given in [49] as parameter set I: = 14.17 and g ρ T /g ρNN = 6.1, and the tiny ω tensor coupling is neglected.We note that this model just describes the deuteron binding energy and not any more the full S-wave phase shifts.For that, one would have to vary the cut-off for each meson exchange.However, for the calculation performed here such a simplified treatment is sufficient.The resulting change in the deuteron binding energy, when we vary the Higgs VEV by ±10%, is much smaller than what was given in [1] and also has the opposite sign, see Fig. 10.Note further that keeping g πNN constant instead of varying it as in Eq.(A.7) has only a mild influence on the results.
Plugging in this deuteron binding energy dependence, we find that the n + p → d + γ rate does not vary as much as before and resembles more the (1 + δv/v) 0.25 -scaling that was used in [1].This is presented in Fig. 11.
The resulting 4 He and d abundances for this calculation can the found in Fig. 12.They are close to the results obtained by simply scaling the n + p → d + γ rate as described in [1], which is evident when comparing Fig. 12 to Fig. 9.
Figure 11: Reaction rate for n+ p → d+γ for δv/v = 0, ±0.1 calculated from pionless effective field theory [45] with nuclear parameters changed as described in the text, but the binding energy taken from a OBEP calculation.The shaded green band gives the original rate mutliplied by (1 ± 0.1) 0.25 , so the range in which the n + p → d + γ rate varies in the calculation of [1].In both papers (and all that we are mentioning in this section), the calculation for the neutron-proton mass difference uses the values from [23] for the electromagnetic and strong contribution, while we use the more recent ones from [24].It is not clear to us if the authors included temperaturedependent effects in the weak n ↔ p rates and effects from a changing deuteron binding energy or other nucleon-nucleon scattering parameters on the n + p → d + γ rate, as well as a varying nucleon mass, were not considered.In [11], the pion mass dependence on the Higgs VEV does not include any variation of the strong scale and for the deuteron binding energy dependence on the pion mass the authors made a simple estimation.In [12] a changing pion mass or deuteron binding energy was not mentioned at all.Yoo and Scherrer [13] then improved the treatment of the deuteron binding energy in 2003, assuming a linear dependence on the pion mass.The authors also derived constraints on the variation of the Higgs VEV from CMB measurements in addition to BBN simulations.Li and Chu (2006) [14] employed much of the methodology from [13], but included now the variation of the deuteron binding energy (but not other scattering parameters) in a pionless EFT approach for the n + p → d + γ rate.
Campbell and Olive studied the dependence of primordial abundances on the VEV of some dilation field in 1994, including effects of a varying Higgs VEV [53].In our understanding they did not, however, implement their methods into a BBN code but used a simple approximated expression for the Helium-4 abundance.Similarly, Dmitriev and Flambaum (2003) [54] and Flambaum and Shuryak (2003) [55] considered only changes in the quark masses and the strong scale and do not use a BBN code but find constraints for example from the Helium-5 binding.Dmitriev, Flambaum and Webb (2004) [56]  for the n + p → d + γ rate.Although Bedaque, Luu and Platter [16] found this scaling to be rather accurate for the deuteron binding energy, it does not include effects from changing other nucleon-nucleon scattering paramters.
A Higgs VEV variation was also part of the work of Coc et al. in 2006 [48].The employ a scaling of the n+p → d+γ rate like in [56] and consider a linear dependence of the deuteron binding energy either on the σand ω-meson mass and the nucleon mass or the same linear dependence on the pion mass as in [13].These authors considerered combined variations of the Higgs VEV and the Yukawa couplings.Using some approximate scaling relations, as they appear in the discussed unified theories, limits on the variations of the Yukawa couplings are displayed.Thus, a direct comparison with our work is not possible.
In their extensive work on the variation of fundamental constants in the BBN framework, Dent, Stern and Wetterich [15], carried out the Higgs VEV variation independently of varying the fermion masses.One would have to adjust the response matrices accordingly to properly compare our results.The authors otherwise used a very similar approach to [13].
Finally, the differences to the works of Bedaque, Luu and Platter [16], Berengut et al. [17] and, of course, to the newest publication of Burns et al. [1] were explained in detail in the text.
The ways we have improved the methodology of the publications listed here is summarised at the end of section 5.

Figure 2 :
Figure 2: Cubic polynomial fit (black line) to lattice QCD data for the deuteron binding energy in the 3 S 1 channel of nucleon-nucleon scattering.The fit was forced to go through the physical point.

Figure 3 :
Figure 3: Cubic polynomial fit (black line) to inverse scattering length resulting from the LET calculations in the 1 S 0 (left panel) and the 3 S 1 channel (right panel) of nucleon-nucleon scattering.The data points from Refs.[34] and[35] are direct lattice QCD results.The fit was forced to go through the physical point.

2 x 2
dx , where x = p

Figure 5 :Figure 6 :
Figure 5: Fractional change in the weak n ↔ p reaction rates including all relevant effects.

Figure 8 :
Figure 8: 4 He (left) and d (right) abundance as a function of δv/v.

Figure 9 :
Figure 9: 4 He (left) and d (right) abundance as a function of δv/v not including changes in the n + p → d + γ rate as would be implied by pionless EFT but scaling it like Λ QCD .

Figure 10 :
Figure10: Deuteron binding energy as a function of δv/v calculated using the OBEP from[49].

Figure 12 :
Figure 12: 4 He (left) and deuterium (right) abundance as a function of δv/v for the deuteron binding energy calculated using an OBEP.

2 d(
consider only the change in the deuteron binding energy and use the simple scaling ⟨σ(n + p → d + γ)v⟩ ∝ B 5/