diff --git a/.github/workflows/check+deploy.yaml b/.github/workflows/check+deploy.yaml index a7631a6..b835e40 100644 --- a/.github/workflows/check+deploy.yaml +++ b/.github/workflows/check+deploy.yaml @@ -63,7 +63,7 @@ jobs: python3 -m pip install dist/wilson-*.whl - name: Upload build as artifact - uses: actions/upload-artifact@v1 + uses: actions/upload-artifact@v4 with: name: wilson-dist-${{ github.sha }} path: dist diff --git a/wilson/classes.py b/wilson/classes.py index 379d3dd..83b86ed 100644 --- a/wilson/classes.py +++ b/wilson/classes.py @@ -1,7 +1,7 @@ """Main classes used at the top level of the wilson package: `Wilson`: main interface to the wilson package, providing automatic running -and matching in SMEFT and WET +and matching in SMEFT and WET and running within nuSMEFT. `RGsolution`: Class representing a continuous solution to the SMEFT and WET RGEs to be used for plotting. @@ -17,6 +17,9 @@ import voluptuous as vol import warnings +from wilson.run.nusmeft import EFTevolve +from wilson.run.nusmeft import beta_nusmeft + class ConfigurableClass: """Class that provides the functionality to set and get configuration options. @@ -107,6 +110,8 @@ class Wilson(ConfigurableClass): # default config options: # dictionary with option name as 'key' and default option value as 'value' + + # dictionary with option name as 'key' and default option value as 'value' _default_options = {'smeft_accuracy': 'integrate', 'qed_order': 1, 'qcd_order': 1, @@ -115,6 +120,8 @@ class Wilson(ConfigurableClass): 'mb_matchingscale': 4.2, 'mc_matchingscale': 1.3, 'parameters': {}, + 'yukawa_scale_in': {}, + 'gauge_higgs_scale_in': {} } # option schema: @@ -128,7 +135,8 @@ class Wilson(ConfigurableClass): 'mb_matchingscale': vol.Coerce(float), 'mc_matchingscale': vol.Coerce(float), 'parameters': vol.Schema({vol.Extra: vol.Coerce(float)}), - }) + 'yukawa_scale_in' : vol.Schema( { vol.In(['Gu', 'Gd', 'Ge', 'Gn']): vol.All(vol.All(vol.Length (min=3, max=3)), vol.Length(min=3, max=3)) } ), + 'gauge_higgs_scale_in': vol.Schema({vol.In(['g', 'gp', 'gs', 'Lambda', 'm2']): vol.Coerce(float)}) }) def __init__(self, wcdict, scale, eft, basis): """Initialize the `Wilson` class. @@ -149,6 +157,12 @@ def __init__(self, wcdict, scale, eft, basis): self.wc.validate() self._cache = {} + def update_dim4_running(self): + par = {} + par.update(self.get_option('yukawa_scale_in')) + par.update(self.get_option('gauge_higgs_scale_in')) + return par + def __hash__(self): """Return a hash of the `Wilson` instance. The hash changes when Wilson coefficient values or options are modified. @@ -239,7 +253,7 @@ def match_run(self, scale, eft, basis, sectors='all'): wc_out = smeft.run(scale, accuracy=smeft_accuracy).translate(basis) self._set_cache('all', scale, 'SMEFT', wc_out.basis, wc_out) return wc_out - else: + elif eft in ['WET', 'WET-4', 'WET-3']: # if SMEFT -> WET-x: match to WET at the EW scale wc_ew = self._get_from_cache(sector='all', scale=scale_ew, eft='WET', basis='JMS') if wc_ew is None: @@ -250,6 +264,24 @@ def match_run(self, scale, eft, basis, sectors='all'): wc_ew = smeft.run(scale_ew, accuracy=smeft_accuracy).match('WET', 'JMS', parameters=self.matching_parameters) self._set_cache('all', scale_ew, wc_ew.eft, wc_ew.basis, wc_ew) wet = WETrunner(wc_ew, **self._wetrun_opt()) + else: + raise ValueError(f"Output EFT {eft} unknown or not supported") + + # nuSMEFT + elif self.wc.eft == 'nuSMEFT': + smeft_accuracy = self.get_option('smeft_accuracy') + if eft == 'nuSMEFT': + # if input and output EFT just nuSMEFT, just run. + nusmeft = EFTevolve(self.wc.translate('Warsaw', sectors=translate_sectors, parameters=self.parameters),beta_nusmeft.nubeta) + nusmeft.ext_par_scale_in(self.update_dim4_running()) + wc_out = nusmeft.run(scale, accuracy=smeft_accuracy).translate(basis) + self._set_cache('all', scale, 'nuSMEFT', wc_out.basis, wc_out) + return wc_out + else: + raise ValueError(f"Output EFT {eft} unknown or not supported") + + + elif self.wc.eft in ['WET', 'WET-4', 'WET-3']: wet = WETrunner(self.wc.translate('JMS', parameters=self.parameters, sectors=translate_sectors), **self._wetrun_opt()) else: diff --git a/wilson/run/__init__.py b/wilson/run/__init__.py index d1b495e..7883b22 100644 --- a/wilson/run/__init__.py +++ b/wilson/run/__init__.py @@ -2,4 +2,5 @@ from . import smeft +from . import nusmeft from . import wet diff --git a/wilson/run/nusmeft/__init__.py b/wilson/run/nusmeft/__init__.py new file mode 100644 index 0000000..5e55ae9 --- /dev/null +++ b/wilson/run/nusmeft/__init__.py @@ -0,0 +1,11 @@ +"""SMEFT and nuSMEFT renormalization group evolution. + +Based on SMEFT: arXiv:1308.2627, arXiv:1310.4838, and arXiv:1312.2014. + nuSMEFT: arxiv:2103.04441, arxiv:2010.12109 + +SMEFT part Ported from the [DsixTools](https://dsixtools.github.io) Mathematica package. +""" +from . import beta_nusmeft +from . import classes +from . import definitions +from .classes import EFTevolve diff --git a/wilson/run/nusmeft/beta_nusmeft.py b/wilson/run/nusmeft/beta_nusmeft.py new file mode 100644 index 0000000..702d75f --- /dev/null +++ b/wilson/run/nusmeft/beta_nusmeft.py @@ -0,0 +1,380 @@ +#new beta function for nuSMEFT + +from wilson.run.smeft.beta import * +from wilson.util import nusmeftutil +import numpy as np + +def nubeta(C, HIGHSCALE=1, newphys=True): + + HIGHSCALE=1 + g = C["g"] + gp = C["gp"] + gs = C["gs"] + m2 = C["m2"] + yn = 0 + Lambda = C["Lambda"] + + Gu = C["Gu"] + Gd = C["Gd"] + Ge = C["Ge"] + Gn = C["Gn"] + +# Gu = np.zeros((3,3)) +# Gd = np.zeros((3,3)) +# Ge = np.zeros((3,3)) +# Gn = np.zeros((3,3)) + + yn = 0 + ye = -1 + yl = -1/2 + yd = -1/3 + yu = 2/3 + yq = 1/6 + yh = 1/2 + Nc = 3 + + CF2 = 3/4 + b02 = 19/6 + b01 = -41/6 + + + + GammaH = np.trace(3*Gu @ Gu.conj().T + 3*Gd @ Gd.conj().T + Ge @ Ge.conj().T) + Gammaq = 1/2*(Gu @ Gu.conj().T + Gd @ Gd.conj().T) + Gammau = Gu.conj().T @ Gu + Gammad = Gd.conj().T @ Gd + Gammal = 1/2*Ge @ Ge.conj().T + Gammae = Ge.conj().T @ Ge + Gamman = Gn.conj().T @ Gn # + + + Xie = 2*my_einsum("prst,rs", C["le"], Ge) \ + - 3*my_einsum("ptsr,rs", C["ledq"], Gd) \ + + 3*my_einsum("ptsr,sr", C["lequ1"], np.conj(Gu)) \ + - my_einsum("vwpr,vw", C["lnle"], np.conj(Gn)) + + Xid = 2*(my_einsum("prst,rs", C["qd1"], Gd) \ + + 4/3*my_einsum("prst,rs", C["qd8"], Gd)) \ + - (3*my_einsum("srpt,sr", C["quqd1"], np.conj(Gu)) \ + + 1/2*(my_einsum("prst,sr", C["quqd1"], np.conj(Gu)) \ + + 4/3*my_einsum("prst,sr", C["quqd8"], np.conj(Gu)))) \ + - my_einsum("srtp,sr", np.conj(C["ledq"]), Ge) \ + - my_einsum("vwpr,vw", C["lnqd1"], np.conj(Gn)) + + Xiu = 2*(my_einsum("prst,rs", C["qu1"], Gu) \ + + 4/3*my_einsum("prst,rs", C["qu8"], Gu)) \ + - (3*my_einsum("ptsr,sr", C["quqd1"], np.conj(Gd)) \ + + 1/2*(my_einsum("stpr,sr", C["quqd1"], np.conj(Gd)) \ + + 4/3*my_einsum("stpr,sr", C["quqd8"], np.conj(Gd)))) \ + + my_einsum("srpt,sr", C["lequ1"], np.conj(Ge)) \ + - my_einsum("vwrp,vw", np.conj(C["lnuq"]), Gn) + + Xin = 2*my_einsum("pwvr,wv", C["ln"], Gn) \ + - 3*my_einsum("prvw,vw", C["lnqd1"], np.conj(Gd)) \ + - 3*my_einsum("prvw,wv", C["lnuq"], Gu) \ + - my_einsum("prvw,vw", C["lnle"], np.conj(Ge)) + + nuBeta = beta(C, HIGHSCALE, newphys) + + if not newphys: + # if there is no new physics, generate a dictionary with zero + # Wilson coefficients (i.e. zero beta functions) + BetaSM = nusmeftutil.C_array2dict(np.zeros(5000)) + BetaSM.update(nuBeta) + return BetaSM + + nuBeta["Gn"] = np.zeros((3,3)) #FIXME? + + nuBeta["ll"] = nuBeta["ll"] + (1/2)*my_einsum("pr,st", Gn @ Gn.conj().T, C["phil1"]) \ + + (1/2)*my_einsum("pr,st", Gn @ Gn.conj().T, C["phil3"]) \ + + (1/2)*my_einsum("st,pr", Gn @ Gn.conj().T, C["phil1"]) \ + + (1/2)*my_einsum("st,pr", Gn @ Gn.conj().T, C["phil3"]) \ + - (1/2)*my_einsum("vs,tw,prvw", np.conj(Gn), np.conj(Gn), C["ln"]) \ + - (1/2)*my_einsum("vp,rw,stvw", np.conj(Gn), np.conj(Gn), C["ln"] ) \ + - (1/2)*my_einsum("rv,tw,pvsw", np.conj(Gn), np.conj(Ge), C["lnle"]) \ + - (1/2)*my_einsum("tw,rv,swpv", np.conj(Gn), np.conj(Ge), C["lnle"]) \ + - (1/2)*my_einsum("pv,sw,rvtw", Gn, Ge, np.conj(C["lnle"])) \ + - (1/2)*my_einsum("sw,pv,twrv", Gn, Ge, np.conj(C["lnle"])) + + nuBeta["lq1"] = nuBeta["lq1"] + my_einsum("pr,st", Gn @ Gn.conj().T, C["phiq1"])\ + - my_einsum("pv,rw,stvw", Gn, np.conj(Gn), C["qn"]) \ + + (1/4)*my_einsum("rv,ws,pvwt", np.conj(Gn), np.conj(Gu), C["lnuq"]) \ + + (1/4)*my_einsum("pv,wt,rvws", Gn, Gu, np.conj(C["lnuq"])) \ + - (1/4)*my_einsum("rv,tw,pvsw", np.conj(Gn), np.conj(Gd), C["lnqd1"]) \ + - (1/4)*my_einsum("pv,sw,rvtw", Gn, Gd, np.conj(C["lnqd1"])) \ + + 3*my_einsum("rv,tw,pvsw", np.conj(Gn), np.conj(Gd), C["lnqd3"]) \ + + 3*my_einsum("pv,sw,rvtw", Gn, Gd, np.conj(C["lnqd3"])) + + + nuBeta["lq3"] = nuBeta["lq3"] + my_einsum("pr,st", Gn @ Gn.conj().T, C["phiq3"]) \ + + (1/4)*my_einsum("rv,ws,pvwt", np.conj(Gn), np.conj(Gu), C["lnuq"]) \ + + (1/4)*my_einsum("pv,wt,rvws", Gn, Gu, np.conj(C["lnuq"])) \ + + (1/4)*my_einsum("rv,tw,pvsw", np.conj(Gn), Gd, C["lnqd1"]) \ + + (1/4)*my_einsum("pv,sw,rvtw", Gn, Gd, np.conj(C["lnqd1"])) \ + - 3*my_einsum("rv,tw,pvsw", np.conj(Gn), np.conj(Gd), C["lnqd3"]) \ + - 3*my_einsum("pv,sw,rvtw", Gn, Gd, C["lnqd3"]) + + nuBeta["ld"] = nuBeta["ld"] + my_einsum("pr,st", Gn @ Gn.conj().T, C["phid"]) \ + - my_einsum("vp,wr,vwst", np.conj(Gn), Gn, C["nd"]) \ + + (1/2)*my_einsum("ws,rv,pvwt", np.conj(Gd), np.conj(Gn), C["lnqd1"]) \ + + (1/2)*my_einsum("wt,pv,rvws", Gd, Gn, np.conj(C["lnqd1"])) \ + + 6*my_einsum("ws,rv,pvwt", np.conj(Gd), np.conj(Gn), C["lnqd3"]) \ + + 6*my_einsum("wt,pv,rvws", Gd, Gn, np.conj(C["lnqd3"])) + + nuBeta["lu"] = nuBeta["lu"] + my_einsum("pr,st", Gn @ Gn.conj().T, C["phiu"]) + - my_einsum("vp,wr,vwst", np.conj(Gn), Gn, C["nu"]) \ + - (1/2)*my_einsum("wt,rv,pvsw", Gu, np.conj(Gn), C["lnuq"]) \ + - (1/2)*my_einsum("ws,pv,rvtw", np.conj(Gu), Gn, np.conj(C["lnuq"])) + + nuBeta["le"] = nuBeta["le"] + my_einsum("rs,pt",np.conj(Ge), Xie) \ + + my_einsum("pt,rs", Ge, np.conj(Xie)) \ + + my_einsum("pr,st", Gn @ Gn.conj().T, C["phie"]) \ + - my_einsum("vp,wr,vwst", np.conj(Gn), Gn, C["ne"]) \ + + (1/2)*my_einsum("ws,rv,pvwt", np.conj(Ge), np.conj(Gn), C["lnle"]) \ + + (1/2)*my_einsum("wt,pv,rvws", Ge, Gn, np.conj(C["lnle"])) + + nuBeta["qu1"] = nuBeta["qu1"] + (1/Nc)*my_einsum("rs,pt", np.conj(Gu), Xiu) \ + + (1/Nc)*my_einsum("pt,rs", Gu, np.conj(Xiu)) + + nuBeta["qu8"] = nuBeta["qu8"] + 2*my_einsum("rs,pt", np.conj(Gu), Xiu) \ + + 2*my_einsum("pt,rs", Gu, np.conj(Xiu)) + + nuBeta["qd1"] = nuBeta["qd1"] + (1/Nc)*my_einsum("rs,pt", np.conj(Gd), Xid) \ + + (1/Nc)*my_einsum("pt,rs", Gd, np.conj(Xid)) + + nuBeta["qd8"] = nuBeta["qd8"] + 2*my_einsum("rs,pt", np.conj(Gd), Xid) \ + + 2*my_einsum("pt,rs", Gd, np.conj(Xid)) + + nuBeta["ledq"] = nuBeta["ledq"] - 2*my_einsum("ts,pr", np.conj(Gd), Xie) \ + - 2*my_einsum("pr,ts", Ge, np.conj(Xid)) \ + + 2*my_einsum("pv,wr,wvst", Gn, Ge, np.conj(C["lnqd1"])) \ + + 2*my_einsum("pv,wt,vrsw", Gn, np.conj(Gu), C["nedu"]) + + nuBeta["lequ1"] = nuBeta["lequ1"] + 2*my_einsum("st,pr", Gu, Xie) \ + + 2*my_einsum("pr,st", Ge, Xiu) \ + + 2*my_einsum("pv,wr,wvts", Gn, Ge, np.conj(C["lnuq"])) \ + - 2*my_einsum("pv,sw,vrwt", Gn, Gd, C["nedu"]) + + nuBeta["lequ3"] = nuBeta["lequ3"] \ + + (1/2)*my_einsum("pv,sw,vrwt", Gn, Gd, C["nedu"]) + + nuBeta["quqd1"] = nuBeta["quqd1"] - 2*my_einsum("pr,st", Gu, Xid) \ + - 2*my_einsum("st,pr", Gd, Xiu) + + #SMNEFT operators + + #RRRR operators + nuBeta["nd"] = (4/3)*Nc*yd**2*gp**2*my_einsum("prww,st", C["nd"], I3)\ + +(4/3)*Nc*yd*yu*gp**2*my_einsum("prww,st", C["nu"], I3) \ + +(4/3)*yd*ye*gp**2*my_einsum("prww,st", C["ne"], I3) \ + +(8/3)*Nc*yd*yq*gp**2*my_einsum("wwpr,st", C["qn"], I3) \ + +(8/3)*yd*yl*gp**2*my_einsum("wwpr,st", C["ln"], I3) \ + +(4/3)*yd*yh*gp**2*my_einsum("pr,st", C["phin"], I3) \ + - 2*my_einsum("pr,st", Gn.conj().T @ Gn, C["phid"]) \ + + 2*my_einsum("st,pr", Gd.conj().T @ Gd, C["phin"]) \ + - 2*my_einsum("vp,wr,vwst", np.conj(Gn), Gn, C["ld"]) \ + - 2*my_einsum("vs,wt,vwpr", np.conj(Gd), Gd, C["qn"]) \ + - 1*my_einsum("vp,ws,vrwt", np.conj(Gn), np.conj(Gd), C["lnqd1"]) \ + - 1*my_einsum("vr,wt,vpws", Gn, Gd, np.conj(C["lnqd1"])) \ + + 12*my_einsum("vp,ws,vrwt", np.conj(Gn), np.conj(Gd), C["lnqd3"]) \ + + 12*my_einsum("vr,wt,vpws", Gn, Gd, np.conj(C["lnqd3"])) \ + + my_einsum("pv,vrst", Gamman, C["nd"]) \ + + my_einsum("sv,prvt", Gammad, C["nd"]) \ + + my_einsum("pvst,vr", C["nd"], Gamman) \ + + my_einsum("prsv,vt", C["nd"], Gammad) + + + nuBeta["nu"] = (4/3)*Nc*yu*yd*gp**2*my_einsum("prww,st", C["nd"], I3) \ + +(4/3)*Nc*yu**2*gp**2*my_einsum("prww,st", C["nu"], I3) \ + +(4/3)*yu*ye*gp**2*my_einsum("prww,st", C["ne"], I3) \ + +(8/3)*Nc*yu*yq*gp**2*my_einsum("wwpr,st", C["qn"], I3) \ + +(8/3)*yu*yl*gp**2*my_einsum("wwpr,st", C["ln"], I3) \ + +(4/3)*yu*yh*gp**2*my_einsum("pr,st", C["phin"], I3) \ + - 2*my_einsum("pr,st", Gn.conj().T @ Gn, C["phiu"]) \ + - 2*my_einsum("st,pr", Gu.conj().T @ Gu, C["phie"]) \ + - 2*my_einsum("vp,wr,vwst", np.conj(Gn), Gn, C["lu"]) \ + - 2*my_einsum("vs,wt,vwpr", np.conj(Gu), Gu, C["qn"]) \ + + my_einsum("vp,wt,vrsw", np.conj(Gn), Gu, C["lnuq"]) \ + + my_einsum("vr,ws,vptw", Gn, np.conj(Gu), np.conj(C["lnuq"])) \ + + my_einsum("pv,vrst", Gamman, C["nu"]) \ + + my_einsum("sv,prvt", Gammau, C["nu"]) \ + + my_einsum("pvst,vr", C["nu"], Gamman) \ + + my_einsum("prsv,vt", C["nu"], Gammau) + + nuBeta["ne"] = (4/3)*Nc*ye*yd*gp**2*my_einsum("prww,st", C["nd"], I3) \ + +(4/3)*Nc*ye*yu*gp**2*my_einsum("prww,st", C["nu"], I3) \ + +(4/3)*ye**2*gp**2*my_einsum("prww,st", C["ne"], I3) \ + +(8/3)*Nc*ye*yq*gp**2*my_einsum("wwpr,st", C["qn"], I3) \ + +(8/3)*ye*yl*gp**2*my_einsum("wwpr,st", C["ln"], I3) \ + +(4/3)*ye*yh*gp**2*my_einsum("pr,st", C["phin"], I3) \ + -2*my_einsum("pr,st",Gn.conj().T @ Gn, C["phie"]) \ + +2*my_einsum("st,pr",Ge.conj().T @ Ge, C["phin"]) \ + +2*my_einsum("sr,pt",Ge.conj().T @ Gn, C["phine"]) \ + +2*my_einsum("pt,rs",Gn.conj().T @ Ge, np.conj(C["phine"])) \ + -2*my_einsum("vp,wr,vwst",np.conj(Gn), Gn, C["le"]) \ + -1*my_einsum("vp,ws,vrwt",np.conj(Gn), np.conj(Ge), C["lnle"]) \ + -1*my_einsum("vr,wt,vpws",Gn, Ge, np.conj(C["lnle"]) ) \ + +1*my_einsum("wp,vs,vrwt",np.conj(Gn), np.conj(Ge), C["lnle"]) \ + +1*my_einsum("wr,vt,vpws",Gn, Ge, np.conj(C["lnle"]) ) \ + -2*my_einsum("vs,wt,vwpr",np.conj(Ge), Ge, C["ln"]) \ + + my_einsum("pv,vrst", Gamman, C["ne"]) \ + + my_einsum("sv,prvt", Gammae, C["ne"]) \ + + my_einsum("pvst,vr", C["ne"], Gamman) \ + + my_einsum("prsv,vt", C["ne"], Gammae) + + nuBeta["nn"] = - my_einsum("pr,st", Gn.conj().T @ Gn, C["phin"]) \ + - my_einsum("st,pr", Gn.conj().T @ Gn, C["phin"]) \ + - my_einsum("vp,wr,vwst", np.conj(Gn), Gn, C["ln"]) \ + - my_einsum("vs,wt,vwpr", np.conj(Gn), Gn, C["ln"]) \ + + my_einsum("pv,vrst", Gamman, C["nn"]) \ + + my_einsum("sv,prvt", Gamman, C["nn"]) \ + + my_einsum("pvst,vr", C["nn"], Gamman) \ + + my_einsum("pvsv,vt", C["nn"], Gamman) + + nuBeta["nedu"] =((yd-yu)**2 + ye*(ye + 8*yu - 2*yd))*gp**2*my_einsum("prst",C["nedu"])\ + + 2*my_einsum("st,pr", Gd.conj().T @ Gu, C["phine"]) \ + + 2*my_einsum("pr,ts", Gn.conj().T @ Ge, C["phiud"]) \ + - 1*my_einsum("vp,ws,vrwt", np.conj(Gn), np.conj(Gd), C["lequ1"]) \ + + 12*my_einsum("vp,ws,vrwt", np.conj(Gn), np.conj(Gd), C["lequ3"]) \ + + 1*my_einsum("vp,wt,vrsw", np.conj(Gn), Gu, C["ledq"]) \ + + 1*my_einsum("vr,wt,vpws", Ge, Gu, np.conj(C["lnqd1"]) ) \ + - 12*my_einsum("vr,wt,vpws", Ge, Gu, np.conj(C["lnqd3"]) ) \ + + 1*my_einsum("vr,ws,vptw", Ge, np.conj(Gd), C["lnuq"]) \ + + my_einsum("pv,vrst", Gamman, C["nedu"]) \ + + my_einsum("sv,prvt", Gammad, C["nedu"]) \ + + my_einsum("pvst,vr", C["nedu"], Gammae) \ + + my_einsum("prsv,vt", C["nedu"], Gammau) + + #LLRR Operators + + nuBeta["qn"] =(4/3)*Nc*yq*yd*gp**2*my_einsum("stww,pr", C["nd"], I3) \ + +(4/3)*Nc*yq*yu*gp**2*my_einsum("stww,pr", C["nu"], I3) \ + +(4/3)*yq*ye*gp**2*my_einsum("stww,pr", C["ne"], I3) \ + +(8/3)*Nc*yq**2*gp**2*my_einsum("wwst,pr", C["qn"], I3) \ + +(8/3)*yq*yl*gp**2*my_einsum("wwst,pr", C["ln"], I3) \ + +(4/3)*yq*yh*gp**2*my_einsum("st,pr", C["phin"], I3) \ + + my_einsum("pv,vrst", Gammaq, C["qn"]) \ + + my_einsum("sv,prvt", Gamman, C["qn"]) \ + + my_einsum("pvst,vr", C["qn"], Gammaq) \ + + my_einsum("prsv,vt", C["qn"], Gamman) + + nuBeta["ln"] = (4/3)*Nc*yl*yd*gp**2*my_einsum("stww,pr", C["nd"], I3) \ + +(4/3)*Nc*yl*yu*gp**2*my_einsum("stww,pr", C["nu"], I3) \ + +(4/3)*yl*ye*gp**2*my_einsum("stww,pr", C["ne"], I3) \ + +(8/3)*Nc*yl*yq*gp**2*my_einsum("wwst,pr", C["qn"], I3) \ + +(8/3)*yl**2*gp**2*my_einsum("wwst,pr", C["ln"], I3) \ + +(4/3)*yl*yh*gp**2*my_einsum("st,pr", C["phin"], I3) \ + + my_einsum("pv,vrst", Gammal, C["ln"]) \ + + my_einsum("sv,prvt", Gamman, C["ln"]) \ + + my_einsum("pvst,vr", C["ln"], Gammal) \ + + my_einsum("prsv,vt", C["ln"], Gamman) + + #LRRL and LRLR operators + + nuBeta["lnle"] = ((ye**2 - 8*ye*yl + 6*yl**2)*gp**2 - (3/2)*g**2)*my_einsum("prst", C["lnle"]) \ + -(4*yl*(ye+yl)*gp**2 - 3*g**2)*my_einsum("srpt",C["lnle"])\ + -4*my_einsum("vr,wt,pvsw", Gn, Ge, C["ll"])\ + +4*my_einsum("vr,wt,svpw", Gn, Ge, C["ll"])\ + +4*my_einsum("wr,vt,pvsw", Gn, Ge, C["ll"])\ + -4*my_einsum("wt,vt,svpw", Gn, Ge, C["ll"])\ + -4*my_einsum("pv,sw,vrwt", Gn, Ge, C["ne"])\ + +4*my_einsum("sv,pw,vrwt", Gn, Ge, C["ne"])\ + +4*my_einsum("sw,vt,pvwr", Gn, Ge, C["ln"])\ + -4*my_einsum("vr,pw,svwt", Gn, Ge, C["le"])\ + +4*gp*(ye+yl)*my_einsum("pr,st", C["nB"], Ge)\ + -8*gp*(ye+yl)*my_einsum("sr,pt", C["nB"], Ge)\ + -6*g*my_einsum("pr,st", C["nB"], Ge)\ + +12*g*my_einsum("sr,pt", C["nW"], Ge)\ + +4*gp*(yn+yl)*my_einsum("st,pr", C["eB"], Gn)\ + -8*gp*(yn+yl)*my_einsum("pt,sr", C["eB"], Gn)\ + -6*g*my_einsum("st,pr", C["eW"], Gn)\ + +12*g*my_einsum("pt,sr", C["eW"], Gn)\ + -2*my_einsum("pr,st", Xin, Ge) \ + -2*my_einsum("st,pr", Xie, Gn) \ + + my_einsum("pv,vrst", Gammal, C["lnle"]) \ + + my_einsum("sv,prvt", Gammal, C["lnle"]) \ + + my_einsum("pvst,vr", C["lnle"], Gamman) \ + + my_einsum("prsv,vt", C["lnle"], Gammae) + + nuBeta["lnqd1"] = ((yd**2 - 2*yd*(yl + 4*yq) + (yl + yq)**2)*gp**2 - 8*gs**2)*my_einsum("prst", C["lnqd1"]) \ + + (-24*yl*(yd + yq)*gp**2 + 18*g**2)*my_einsum("prst", C["lnqd3"])\ + -2*my_einsum("vr,pw,vwts", Gn, Ge, np.conj(C["ledq"]))\ + +2*my_einsum("pw,vt,svwr", Gn, Gd, C["qn"])\ + +2*my_einsum("pw,sv,rwtv", Ge, Gu, np.conj(C["nedu"]))\ + +2*my_einsum("vr,sw,pvwt", Gn, Gd, C["ld"])\ + -2*my_einsum("wt,sv,prvw", Gd, Gu, np.conj(C["lnuq"]))\ + -2*my_einsum("pw,sv,wrvt", Gn, Gd, C["nd"])\ + -2*my_einsum("vr,wt,pvsw", Gn, Gd, np.conj(C["lq1"]))\ + +6*my_einsum("vr,wt,pvsw", Gn, Gd, C["lq3"])\ + -2*my_einsum("pr,st", Xin, Gd) \ + -2*my_einsum("st,pr", Xid, Gn) \ + + my_einsum("pv,vrst", Gammal, C["lnqd1"]) \ + + my_einsum("sv,prvt", Gammaq, C["lnqd1"]) \ + + my_einsum("pvst,vr", C["lnqd1"], Gamman) \ + + my_einsum("prsv,vt", C["lnqd1"], Gammad) + + + nuBeta["lnqd3"] = (-(1/2)*yl*(yd+yq)*gp**2 + (3/8)*g**2)*my_einsum("prst",C["lnqd1"])\ + +((yd**2 - 6*yd*yl + yl**2 + 6*yl*yq + yq**2)*gp**2 - 3*g**2 + (8/3)*gs**2)*my_einsum("prst",C["lnqd3"])\ + -(1/2)*my_einsum("pw,sv,rwtv", Ge, Gu, np.conj(C["nedu"]))\ + +(1/2)*my_einsum("vr,wt,pvsw", Gn, Gd, C["lq1"])\ + -(3/2)*my_einsum("vr,wt,pvsw", Gn, Gd, C["lq3"])\ + +(1/2)*my_einsum("vr,sw,pvwt", Gn, Gd, C["ld"])\ + +(1/2)*my_einsum("pw,vt,svwr", Gn, Gd, C["qn"])\ + +(1/2)*my_einsum("pw,sv,wrvt", Gn, Gd, C["nd"])\ + -gp*(yd+yq)*my_einsum("pr,st", C["nB"], Gd)\ + -gp*(yn+yl)*my_einsum("st,pr", C["dB"], Gn)\ + +(3/2)*g*my_einsum("pr,st", C["nW"], Gd)\ + +(3/2)*g*my_einsum("st,pr", C["dW"], Gn)\ + + my_einsum("pv,vrst", Gammal, C["lnqd3"]) \ + + my_einsum("sv,prvt", Gammaq, C["lnqd3"]) \ + + my_einsum("pvst,vr", C["lnqd3"], Gamman) \ + + my_einsum("prsv,vt", C["lnqd3"], Gammad) + + + + nuBeta["lnuq"] = (((yl+yu)**2 + yq*(yq - 2*yl -8*yu))*gp**2 - 8*gs**2)*my_einsum("prst",C["lnuq"]) \ + +2*my_einsum("wr,pv,wvts", Gn, Ge, np.conj(C["lequ1"])) \ + -2*my_einsum("pw,vs,vtwr", Gn, np.conj(Gu), C["qn"])\ + +2*my_einsum("pw,tv,wrsv", Gn, np.conj(Gu), C["nu"])\ + +2*my_einsum("pw,tv,rwvs", Ge, np.conj(Gd), np.conj(C["nedu"]))\ + +2*my_einsum("vr,tw,pvwt", Gn, np.conj(Gu), C["lq1"])\ + +6*my_einsum("vr,tw,pvwt", Gn, np.conj(Gu), C["lq3"])\ + -2*my_einsum("vr,tw,pvsw", Gn, np.conj(Gu), C["lu"])\ + -2*my_einsum("sv,wt,prvw", Gu, np.conj(Gd), C["lnqd1"])\ + -2*my_einsum("pr,ts", Xin, np.conj(Gu)) \ + -2*my_einsum("st,pr", np.conj(Xiu), Gn) \ + + my_einsum("pv,vrst", Gammal, C["lnuq"]) \ + + my_einsum("sv,prvt", Gammau, C["lnuq"]) \ + + my_einsum("pvst,vr", C["lnuq"], Gamman) \ + + my_einsum("prsv,vt", C["lnuq"], Gammaq) + + + + #psi2phi3 tems + nuBeta["nphi"] = -(9*yl**2*gp**2 + (27/4)*g**2)*my_einsum("pr",C["nphi"]) \ + -6*(4*yh**2*yl*gp**3 - yh*gp*g**2)*my_einsum("pr",C["nB"]) \ + +3*(4*yh*yl*gp**2*g - 3*g**3)*my_einsum("pr",C["nW"]) + + + #psi2phi2D terms + + nuBeta["phin"] = (4/3)*yh**2*gp**2*my_einsum("pr",C["phin"])\ + +(4/3)*Nc*yd*yh*gp**2*my_einsum("prww",C["nd"])\ + +(4/3)*Nc*yu*yh*gp**2*my_einsum("prww", C["nu"]) \ + +(4/3)*ye*yh*gp**2*my_einsum("prww", C["ne"]) \ + +(8/3)*Nc*yq*yh*gp**2*my_einsum("wwpr", C["qn"]) \ + +(8/3)*yl*yh*gp**2*my_einsum("wwpr", C["ln"]) + + #psi2Xphi + nuBeta["nW"] = ((3*CF2 - b02)*g**2 - 3*yl**2*gp**2 )*my_einsum("pr",C["nW"])\ + +3*yl*gp*g*my_einsum("pr",C["nB"]) + (g/4)*my_einsum("srpt,ts", C["lnle"], Ge) + 2*g*Nc*my_einsum("prst,st", C["lnqd3"], Gd) + + nuBeta["nB"] = (-3*CF2*g**2 + (3*yl**2 - b01)*gp**2 )*my_einsum("pr",C["nB"]) \ + +12*CF2*yl*gp*g*my_einsum("pr",C["nW"]) + + + nuBeta["phine"] = -3*ye**2*my_einsum("pr", C["phine"])*gp**2 + + + return nuBeta diff --git a/wilson/run/nusmeft/classes.py b/wilson/run/nusmeft/classes.py new file mode 100644 index 0000000..f9ea687 --- /dev/null +++ b/wilson/run/nusmeft/classes.py @@ -0,0 +1,318 @@ +"""Defines the EFTevolve class that provides the main API to smeft and nusmeft.""" +from wilson.run.smeft import smpar +from wilson.run.nusmeft import definitions +import numpy as np +import ckmutil.phases, ckmutil.diag +import wilson +from wilson.util import smeftutil, nusmeftutil, wetutil +from wilson import wcxf +from copy import deepcopy +from math import pi, log, sqrt +from scipy.integrate import solve_ivp +from collections import OrderedDict + +eftutil = {'SMEFT': smeftutil, 'nuSMEFT': nusmeftutil, 'WET': wetutil} + +class EFTevolve: + """?? Class representing a parameter point in the Standard Model Effective + Field Theory and other EFTs such as nuSMEFT and allowing the evolution + of the Wilson Coefficients. + + Methods: + + - __init__: Initialize, given a wcxf.WC instance + - run: solve the RGE and return a wcxf.WC instance + """ + + def __init__(self, wc, beta, get_smpar=True): + """Initialize the SMEFT instance. + + Parameters: + + - `wc`: the Wilson coefficients as `wcxf.WC` instance. + `beta`: is a function of WCs which returns beta-fucntion dictionary + for SMEFT (beta.beta) or nuSMEFT (beta_numseft.nubeta) + """ + self.wc = wc +# self.scale_in = wc.scale + self.beta = beta + + self.scale_in = None + self.C_in = None + + if wc is not None: + self.basis = wc.basis + self.eft = wc.eft + self._set_initial_wcxf(wc, get_smpar=get_smpar) + + def ext_par_scale_in(self,par): + if par is not None: + self.C_in.update(par) + + def _set_initial(self, C_in, scale_in): + r"""Set the initial values for parameters and Wilson coefficients at + the scale `scale_in`.""" + self.C_in = C_in + self.scale_in = scale_in + + def _set_initial_wcxf(self, wc, get_smpar=True): + """Load the initial values for Wilson coefficients from a + wcxf.WC instance. + + Parameters: + + - `get_smpar`: boolean, optional, defaults to True. If True, an attempt + is made to determine the SM parameters from the requirement of + reproducing the correct SM masses and mixings at the electroweak + scale. As approximations are involved, the result might or might not + be reliable, depending on the size of the Wilson coefficients + affecting the SM masses and mixings. If False, Standard Model + parameters have to be provided separately and are assumed to be in + the weak basis used for the Warsaw basis as defined in WCxf, + i.e. in the basis where the down-type and charged lepton mass + matrices are diagonal. + """ + if wc.eft not in ('SMEFT','nuSMEFT'): + raise ValueError("Wilson coefficients use wrong or unknown EFT.") + if wc.basis not in ('Warsaw'): + raise ValueError("Wilson coefficients use wrong or unknown basis.") + self.scale_in = wc.scale + C = eftutil[self.eft].wcxf2arrays_symmetrized(wc.dict) + # fill in zeros for missing WCs + for k, s in eftutil[self.eft].C_keys_shape.items(): + if k not in C and k not in eftutil[self.eft].dim4_keys: + if s == 1: + C[k] = 0 + else: + C[k] = np.zeros(s) + if self.C_in is None: + self.C_in = C + else: + self.C_in.update(C) + if get_smpar: + self.C_in.update(self._get_sm_scale_in()) + + + def _to_wcxf(self, C_out, scale_out): + """Return the Wilson coefficients `C_out` as a wcxf.WC instance. + + Note that the Wilson coefficients are rotated into the Warsaw basis + as defined in WCxf, i.e. to the basis where the down-type and charged + lepton mass matrices are diagonal.""" + C = self._rotate_defaultbasis(C_out) + d = eftutil[self.eft].arrays2wcxf_nonred(C) + d = wcxf.WC.dict2values(d) + wc = wcxf.WC(self.eft, self.basis, scale_out, d) + return wc + + def beta_reordered(self, C, HIGHSCALE=1, *args, **kwargs): + """Returns a beta dictionary but changes its order to match that is obtained by + EFTutil class. It is essential, because when flattening the array the + order is lost and it is assumed that keys from beta and keys from EFTutil are in order. + Subsequently, they are used in evolution of WCs, where a mismatch of flattened arrays + could lead to incorrect solutions, if not reordered""" + ordered_beta_dict = {} + _keys = eftutil[self.eft].C_keys + _beta = self.beta(C, HIGHSCALE, *args, **kwargs) + ordered_beta_dict = OrderedDict({k: _beta[k] for k in _keys}) + return ordered_beta_dict + + def beta_array_reordered(self, C, HIGHSCALE=1, *args, **kwargs): + """Return the beta functions of all SM parameters and SMEFT/nuSMEFT Wilson + coefficients as a 1D numpy array.""" + return np.hstack([np.asarray(b).ravel() for b in self.beta_reordered(C, HIGHSCALE=1, *args, **kwargs).values()]) + + def eft_evolve_leadinglog(self, scale_out, newphys=True): + """Solve the EFT (SMEFT and nuSMEFT) RGEs in the leading log approximation. + Input C_in and output C_out are dictionaries of arrays.""" + self._check_initial() + C_out = deepcopy(self.C_in) + b = self.beta_reordered(self.C_in, newphys=newphys) + for k, C in C_out.items(): + C_out[k] = C + b[k] / (16 * pi**2) * log(scale_out / self.scale_in) + return C_out + + def _eft_evolve(self, scale_out, newphys=True, **kwargs): + """Axuliary function used in `eft_evolve` and `eft_evolve_continuous`""" + def fun(t0, y): + return self.beta_array_reordered(C=eftutil[self.eft].C_array2dict(y.view(complex)), + newphys=newphys).view(float) / (16 * pi**2) + y0 = eftutil[self.eft].C_dict2array(self.C_in).view(float) + sol = solve_ivp(fun=fun, + t_span=(log(self.scale_in), log(scale_out)), + y0=y0, **kwargs) + return sol + + def eft_evolve(self, scale_out, newphys=True, **kwargs): + """Solve the SMEFT/nuSMEFT RGEs by numeric integration. + + Input C_in and output C_out are dictionaries of arrays.""" + self._check_initial() + sol = self._eft_evolve(scale_out, newphys=newphys, **kwargs) + C_out= eftutil[self.eft].C_array2dict(sol.y[:, -1].view(complex)) + return C_out + + def run(self, scale, accuracy='integrate', **kwargs): + """Return the Wilson coefficients (as wcxf.WC instance) evolved to the + scale `scale`. + + Parameters: + + - `scale`: scale in GeV + - accuracy: whether to use the numerical solution to the RGE + ('integrate', the default, slow but precise) or the leading logarithmic + approximation ('leadinglog', approximate but much faster). + """ + if accuracy == 'integrate': + C_out = self.eft_evolve(scale, **kwargs) + elif accuracy == 'leadinglog': + C_out = self.eft_evolve_leadinglog(scale) + else: + raise ValueError(f"'{accuracy}' is not a valid value of 'accuracy' (must be either 'integrate' or 'leadinglog').") + return self._to_wcxf(C_out, scale) + + + def _run_sm_scale_in(self, C_out, scale_sm=91.1876): + """Get the dim4 parameters at the EW scale according to EFT, using an estimate `C_out` + of the Wilson coefficients at that scale, and run them to the input scale.""" + """New: Before, this function ran only SM parameters. But now it makes more sense to + call them dim4 parameters, as for nusmeft,it will contain Gn, which technically is not + a SM parameter""" + + # set the SM parameters to the values obtained from smpar.eftpar + eftpar = {'SMEFT': smpar.smeftpar, 'nuSMEFT': smpar.nusmeftpar} #FIXME Correct Gn in smpar file + C_SM = eftpar[self.eft](scale_sm, C_out, basis='Warsaw') + SM_keys = set(eftutil[self.eft].dim4_keys) # to speed up lookup + C_SM = {k: v for k, v in C_SM.items() if k in SM_keys} + + # set the Wilson coefficients at the EW scale to C_out + ##First defining C_in_sm - a dictionary with all SMpar and WCs with zeros + C_in_sm = eftutil[self.eft].C_array2dict(np.zeros(9999)) + C_in_sm.update(C_out) + ##Then we update the value of SM parameters which we extracted above + C_in_sm.update(C_SM) + + # initialize an empty EFT instance + eft_sm = EFTevolve(None, self.beta) + # update C_in in eft_sm + eft_sm.C_in = C_in_sm + eft_sm.scale_in = scale_sm + eft_sm.eft = self.eft + + # run up (with 1% relative precision, ignore running of Wilson coefficients) + # putting newphys false makes beta function of new physics WC zero, but it doesnt + # change the beta functions of SM. Therefore, the SM function evolves as usual - + # and it evolves using the C_out. newphys = False doesn't make C_out zero, but it + # makes its corresponding beta function zero. + C_SM_high = eft_sm.eft_evolve(self.scale_in, newphys=False, rtol=0.001, atol=1) + C_SM_high = self._rotate_defaultbasis(C_SM_high) + return {k: v for k, v in C_SM_high.items() if k in SM_keys} + + def _get_sm_scale_in(self, scale_sm=91.1876): + """Get an estimate of the SM parameters at the input scale by running + them from the EW scale using constant values for the Wilson coefficients + (corresponding to their leading log approximated values at the EW + scale). + + Note that this is not guaranteed to work and will fail if some of the + Wilson coefficients (the ones affecting the extraction of SM parameters) + are large.""" + # intialize a copy of ourselves without SM parameters + _smeft = EFTevolve(self.wc, self.beta, get_smpar=False) + # Step 1: run the SM up, using the WCs at scale_input as (constant) estimate + _smeft.C_in.update(self._run_sm_scale_in(self.C_in, scale_sm=scale_sm)) + # Step 2: run the WCs down in LL approximation + C_out = _smeft.eft_evolve_leadinglog(scale_sm) + # Step 3: run the SM up again, this time using the WCs at scale_sm as (constant) estimate + return self._run_sm_scale_in(C_out, scale_sm=scale_sm) + + + + + +# @staticmethod + def _rotate_defaultbasis(self, C): + """Rotate all parameters to the basis where the running down-type quark + and charged lepton mass matrices are diagonal and where the running + up-type quark mass matrix has the form V.S, with V unitary and S real + diagonal, and where the CKM and PMNS matrices have the standard + phase convention.""" + "Equation 3.27 1704.03888" + "Need to implement neutrino sector" #FIXME + v = 246.22 + Mep = v/sqrt(2) * (C['Ge'] - C['ephi'] * v**2/2) + Mup = v/sqrt(2) * (C['Gu'] - C['uphi'] * v**2/2) + Mdp = v/sqrt(2) * (C['Gd'] - C['dphi'] * v**2/2) + Mnup = -v**2 * C['llphiphi'] + UeL, Me, UeR = ckmutil.diag.msvd(Mep) + UuL, Mu, UuR = ckmutil.diag.msvd(Mup) + UdL, Md, UdR = ckmutil.diag.msvd(Mdp) + Unu, Mnu = ckmutil.diag.mtakfac(Mnup) + UuL, UdL, UuR, UdR = ckmutil.phases.rephase_standard(UuL, UdL, UuR, UdR) + Unu, UeL, UeR = ckmutil.phases.rephase_pmns_standard(Unu, UeL, UeR) + if self.eft == 'SMEFT': + C_out = definitions.flavor_rotation_smeft(C, Uq=UdL, Uu=UuR, Ud=UdR, Ul=UeL, Ue=UeR) + elif self.eft == 'nuSMEFT': + C_out = definitions.flavor_rotation_nusmeft(C, Uq=UdL, Uu=UuR, Ud=UdR, Ul=UeL, Ue=UeR, Un=np.eye(3)) + return C_out + + # If we don't input scale_in or input wilson coefficient, the wcxf.WC function will automatically + # not run, so this _check_initial function seems useles. + # Sill will see later if its used. I am just thinking, if the original author put it, + # so mt must have some use + + def get_smpar(self, accuracy='integrate', scale_sm=91.1876): + """Compute the SM MS-bar parameters at the electroweak scale. + + This method can be used to validate the accuracy of the iterative + extraction of SM parameters. If successful, the values returned by this + method should agree with the values in the dictionary + `wilson.run.smeft.smpar.p`.""" + if accuracy == 'integrate': + C_out = self.eft_evolve(scale_sm) + elif accuracy == 'leadinglog': + C_out = self.eft_evolve_leadinglog(scale_sm) + else: + raise ValueError(f"'{accuracy}' is not a valid value of 'accuracy' (must be either 'integrate' or 'leadinglog').") + return smpar.smpar(C_out) + + def eft_evolve_continuous(self, C_in, scale_in, scale_out, newphys=True, **kwargs): + """Solve the SMEFT RGEs by numeric integration, returning a function that + allows to compute an interpolated solution at arbitrary intermediate + scales.""" + sol = self._eft_evolve(scale_out, newphys=newphys, dense_output=True, **kwargs) + @np.vectorize + def _rge_solution(scale): + t = log(scale) + y = sol.sol(t).view(complex) + yd = eftutil[self.eft].C_array2dict(y) + yw = eftutil[self.eft].arrays2wcxf_nonred(yd) + return yw + def rge_solution(scale): + # this is to return a scalar if the input is scalar + return _rge_solution(scale)[()] + return rge_solution + + def run_continuous(self, scale): + """Return a continuous solution to the RGE as `RGsolution` instance.""" + if scale == self.scale_in: + raise ValueError("The scale must be different from the input scale") + elif scale < self.scale_in: + scale_min = scale + scale_max = self.scale_in + elif scale > self.scale_in: + scale_max = scale + scale_min = self.scale_in + fun = self.eft_evolve_continuous(C_in=self.C_in, + scale_in=self.scale_in, + scale_out=scale) + return wilson.classes.RGsolution(fun, scale_min, scale_max) + + def _check_initial(self): + """Check if initial values and scale as well as the new physics scale + have been set.""" + if self.C_in is None: + raise Exception("You have to specify the initial conditions first.") + if self.scale_in is None: + raise Exception("You have to specify the initial scale first.") + diff --git a/wilson/run/nusmeft/definitions.py b/wilson/run/nusmeft/definitions.py new file mode 100644 index 0000000..9f641d5 --- /dev/null +++ b/wilson/run/nusmeft/definitions.py @@ -0,0 +1,109 @@ +"""Definitions of auxiliary objects and operator properties.""" +import numpy as np +from wilson.util import smeftutil +from wilson.util import nusmeftutil + + +def flavor_rotation_smeft(C_in, Uq, Uu, Ud, Ul, Ue, sm_parameters=True): + """Gauge-invariant $U(3)^5$ flavor rotation of all Wilson coefficients and + SM parameters.""" + C = {} + if sm_parameters: + # nothing to do for scalar SM parameters + for k in ['g', 'gp', 'gs', 'Lambda', 'm2']: + C[k] = C_in[k] + C['Ge'] = Ul.conj().T @ C_in['Ge'] @ Ue + C['Gu'] = Uq.conj().T @ C_in['Gu'] @ Uu + C['Gd'] = Uq.conj().T @ C_in['Gd'] @ Ud + # nothing to do for purely bosonic operators + for k in smeftutil.WC_keys_0f: + C[k] = C_in[k] + # see 1704.03888 table 4 (but staying SU(2) invariant here) + # LR + for k in ['ephi', 'eW', 'eB']: + C[k] = Ul.conj().T @ C_in[k] @ Ue + for k in ['uphi', 'uW', 'uB', 'uG']: + C[k] = Uq.conj().T @ C_in[k] @ Uu + for k in ['dphi', 'dW', 'dB', 'dG']: + C[k] = Uq.conj().T @ C_in[k] @ Ud + # LL + for k in ['phil1', 'phil3']: + C[k] = Ul.conj().T @ C_in[k] @ Ul + for k in ['phiq1', 'phiq3']: + C[k] = Uq.conj().T @ C_in[k] @ Uq + C['llphiphi'] = Ul.T @ C_in['llphiphi'] @ Ul + # RR + C['phie'] = Ue.conj().T @ C_in['phie'] @ Ue + C['phiu'] = Uu.conj().T @ C_in['phiu'] @ Uu + C['phid'] = Ud.conj().T @ C_in['phid'] @ Ud + C['phiud'] = Uu.conj().T @ C_in['phiud'] @ Ud + # 4-fermion + C['ll'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Ul, Ul.conj(), Ul.conj(), C_in['ll']) + C['ee'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Ue, Ue.conj(), Ue.conj(), C_in['ee']) + C['le'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Ue, Ul.conj(), Ue.conj(), C_in['le']) + C['qq1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uq, Uq.conj(), Uq.conj(), C_in['qq1']) + C['qq3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uq, Uq.conj(), Uq.conj(), C_in['qq3']) + C['dd'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ud, Ud, Ud.conj(), Ud.conj(), C_in['dd']) + C['uu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Uu, Uu.conj(), Uu.conj(), C_in['uu']) + C['ud1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uu.conj(), Ud.conj(), C_in['ud1']) + C['ud8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uu.conj(), Ud.conj(), C_in['ud8']) + C['qu1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uu, Uq.conj(), Uu.conj(), C_in['qu1']) + C['qu8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Uu, Uq.conj(), Uu.conj(), C_in['qu8']) + C['qd1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ud, Uq.conj(), Ud.conj(), C_in['qd1']) + C['qd8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ud, Uq.conj(), Ud.conj(), C_in['qd8']) + C['quqd1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uq.conj(), Uq.conj(), C_in['quqd1']) + C['quqd8'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ud, Uq.conj(), Uq.conj(), C_in['quqd8']) + C['lq1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Uq, Ul.conj(), Uq.conj(), C_in['lq1']) + C['lq3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Uq, Ul.conj(), Uq.conj(), C_in['lq3']) + C['ld'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Ud, Ul.conj(), Ud.conj(), C_in['ld']) + C['lu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Uu, Ul.conj(), Uu.conj(), C_in['lu']) + C['qe'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ue, Uq.conj(), Ue.conj(), C_in['qe']) + C['ed'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Ud, Ue.conj(), Ud.conj(), C_in['ed']) + C['eu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Ue.conj(), Uu.conj(), C_in['eu']) + C['ledq'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uq, Ul.conj(), Ud.conj(), C_in['ledq']) + C['lequ1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Ul.conj(), Uq.conj(), C_in['lequ1']) + C['lequ3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Ul.conj(), Uq.conj(), C_in['lequ3']) + C['duql'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ul, Ud, Uq, C_in['duql']) + C['qque'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ue, Uq, Uu, C_in['qque']) + C['qqql'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Ul, Uq, Uq, C_in['qqql']) + C['duue'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uu, Ue, Ud, Uu, C_in['duue']) + return C + +#assumption - Un is unitary matrix (mass and flavour basis are same). +def flavor_rotation_nusmeft(C_in, Uq, Uu, Ud, Ul, Ue, Un, sm_parameters=True): + """Gauge-invariant $U(3)^6$ flavor rotation of all Wilson coefficients and + SM parameters.""" + C = {} + + # smeft + C = flavor_rotation_smeft(C_in, Uq, Uu, Ud, Ul, Ue, sm_parameters) + + if sm_parameters: + C['Gn'] = Ul.conj().T @ C_in['Gn'] @ Un #TODO + + # 2-fermion nusmeft + for k in ['nphi', 'nW', 'nB']: # new terms added 20 oct + C[k] = Ul.conj().T @ C_in[k] @ Un + + C['phin'] = Un.conj().T @ C_in['phin'] @ Un # + C['phine'] = Un.conj().T @ C_in['phine'] @ Ue # + + # 4-fermion nusmeft + C['nd'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Ud, Un.conj(), Ud.conj(), C_in['nd']) # + C['nu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Uu, Un.conj(), Uu.conj(), C_in['nu']) # + C['ln'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ul, Un, Ul.conj(), Un.conj(), C_in['ln']) # + C['qn'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Uq, Un, Uq.conj(), Un.conj(), C_in['qn'])# + C['ne'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Ue, Un.conj(), Ue.conj(), C_in['ne']) # + C['nn'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Un, Un.conj(), Un.conj(), C_in['nn']) # + C['nedu'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Ue, Uu, Un.conj(), Ud.conj(), C_in['nedu']) # + C['lnle'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Ue, Ul.conj(), Ul.conj(), C_in['lnle'])# + C['lnqd1'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Ud, Ul.conj(), Uq.conj(), C_in['lnqd1'])# + C['lnqd3'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Ud, Ul.conj(), Uq.conj(), C_in['lnqd3'])# + C['lnuq'] = np.einsum('jb,ld,ia,kc,ijkl->abcd', Un, Uq, Ul.conj(), Uu.conj(), C_in['lnuq'])# + return C + + + + + + diff --git a/wilson/run/smeft/smpar.py b/wilson/run/smeft/smpar.py index fe387f1..d6833f0 100644 --- a/wilson/run/smeft/smpar.py +++ b/wilson/run/smeft/smpar.py @@ -188,3 +188,11 @@ def smpar(C): sm['m_mu'] = Me[1] sm['m_tau'] = Me[2] return {k: v.real for k, v in sm.items()} + +def nusmeftpar(scale, C, basis): + c=smeftpar(scale, C, basis) + GF = p['GF'] + vb = sqrt(1 / sqrt(2) / GF) + v = vb + c['Gn'] = np.zeros((3,3))#C['nphi'] * (v**2) / 2 + return c diff --git a/wilson/util/__init__.py b/wilson/util/__init__.py index 888b496..385f9fc 100644 --- a/wilson/util/__init__.py +++ b/wilson/util/__init__.py @@ -20,6 +20,27 @@ 1: ['Gu', 'Gd', 'Ge'], }, ) + +nusmeftutil = EFTutil( + 'nuSMEFT', + 'Warsaw', + dim4_keys_shape = { + 'g': 1, + 'gp': 1, + 'gs': 1, + 'Lambda': 1, + 'm2': 1, + 'Gu': (3, 3), + 'Gd': (3, 3), + 'Ge': (3, 3), + 'Gn': (3, 3), + }, + dim4_symm_keys = { + 0: ['g', 'gp', 'gs', 'Lambda', 'm2'], + 1: ['Gu', 'Gd', 'Ge','Gn'], + }, +) + wetutil = EFTutil( 'WET', 'JMS',