Skip to content

Commit 233a199

Browse files
jorindevandevisjorindevandevis
authored
Merge pull request #136 from Wall-Go/checkThermo
Check if phases are not identical and if the sound speeds are real
2 parents f00076e + 0ce1849 commit 233a199

File tree

3 files changed

+277
-92
lines changed

3 files changed

+277
-92
lines changed

‎Models/StandardModel/StandardModel.py

+237-89
Original file line numberDiff line numberDiff line change
@@ -1,23 +1,156 @@
11
import numpy as np
2+
import numpy.typing as npt
3+
import os
4+
import pathlib
25

36
## WallGo imports
7+
import WallGo ## Whole package, in particular we get WallGo.initialize()
48
from WallGo import GenericModel
59
from WallGo import Particle
6-
from WallGo import EffectivePotential
710
from WallGo import WallGoManager
11+
## For Benoit benchmarks we need the unresummed, non-high-T potential:
12+
from WallGo import EffectivePotential
13+
from WallGo import Fields
814

915
### LN: This file is very WIP, test with SingletStandardModel_Z2.py instead
1016

17+
## V = msq |phi|^2 + lambda (|phi|^2)^2
18+
class StandardModel(GenericModel):
19+
20+
particles = []
21+
outOfEquilibriumParticles = []
22+
modelParameters = {}
23+
24+
## Specifying this is REQUIRED
25+
fieldCount = 1
26+
27+
28+
def __init__(self, initialInputParameters: dict[str, float]):
29+
30+
self.modelParameters = self.calculateModelParameters(initialInputParameters)
31+
32+
# Initialize internal Veff with our params dict. @todo will it be annoying to keep these in sync if our params change?
33+
self.effectivePotential = EffectivePotentialSM(self.modelParameters, self.fieldCount)
34+
35+
## Define particles. this is a lot of clutter, especially if the mass expressions are long,
36+
## so @todo define these in a separate file?
37+
38+
# NB: particle multiplicity is pretty confusing because some internal DOF counting is handled internally already.
39+
# Eg. for SU3 gluons the multiplicity should be 1, NOT Nc^2 - 1.
40+
# But we nevertheless need something like this to avoid having to separately define up, down, charm, strange, bottom
41+
42+
## === Top quark ===
43+
topMsqVacuum = lambda fields: 0.5 * self.modelParameters["yt"]**2 * fields[0]**2
44+
topMsqThermal = lambda T: self.modelParameters["g3"]**2 * T**2 / 6.0
45+
46+
topQuark = Particle("top",
47+
msqVacuum = topMsqVacuum,
48+
msqThermal = topMsqThermal,
49+
statistics = "Fermion",
50+
inEquilibrium = False,
51+
ultrarelativistic = True,
52+
multiplicity = 1
53+
)
54+
self.addParticle(topQuark)
55+
56+
## === Light quarks, 5 of them ===
57+
lightQuarkMsqThermal = lambda T: self.modelParameters["g3"]**2 * T**2 / 6.0
58+
59+
lightQuark = Particle("lightQuark",
60+
msqVacuum = 0.0,
61+
msqThermal = lightQuarkMsqThermal,
62+
statistics = "Fermion",
63+
inEquilibrium = True,
64+
ultrarelativistic = True,
65+
multiplicity = 5
66+
)
67+
self.addParticle(lightQuark)
68+
69+
## === SU(3) gluon ===
70+
gluonMsqThermal = lambda T: self.modelParameters["g3"]**2 * T**2 * 2.0
71+
72+
gluon = Particle("gluon",
73+
msqVacuum = 0.0,
74+
msqThermal = gluonMsqThermal,
75+
statistics = "Boson",
76+
inEquilibrium = True,
77+
ultrarelativistic = True,
78+
multiplicity = 1
79+
)
80+
self.addParticle(gluon)
81+
82+
## Go from whatever input params --> action params
83+
## This function was just copied from SingletStandardModel_Z2
84+
def calculateModelParameters(self, inputParameters: dict[str, float]) -> dict[str, float]:
85+
super().calculateModelParameters(inputParameters)
86+
87+
modelParameters = {}
88+
89+
v0 = inputParameters["v0"]
90+
# Higgs zero-temperature mass
91+
mH = inputParameters["mH"]
92+
93+
## this are direct input:
94+
modelParameters["RGScale"] = inputParameters["RGScale"]
95+
96+
modelParameters["lambda"] = 0.5 * mH**2 / v0**2
97+
#modelParameters["msq"] = -mh1**2 / 2. # should be same as the following:
98+
modelParameters["msq"] = -modelParameters["lambda"] * v0**2
99+
100+
## Then the gauge/Yukawa sector
101+
102+
Mt = inputParameters["Mt"]
103+
MW = inputParameters["MW"]
104+
MZ = inputParameters["MZ"]
105+
106+
# helper
107+
g0 = 2.*MW / v0
108+
modelParameters["g1"] = g0*np.sqrt((MZ/MW)**2 - 1)
109+
modelParameters["g2"] = g0
110+
# Just take QCD coupling as input
111+
modelParameters["g3"] = inputParameters["g3"]
112+
113+
modelParameters["yt"] = np.sqrt(1./2.)*g0 * Mt/MW
114+
115+
return modelParameters
116+
117+
118+
119+
# ## Define parameter dict here. @todo correct values, these are from top of my head.
120+
# ## In practice the user would define a function that computes these from whatever their input is (mW, mt, mH etc).
121+
# ## But this still needs to be initialized here as required by the interface. modelParameters = None is fine.
122+
# modelParameters = {
123+
# "RGScale" : 91, # MS-bar scale. Units: GeV
124+
# "yt" : 1.0, # Top Yukawa
125+
# "g1" : 0.3, # U1 coupling
126+
# "g2" : 0.6, # SU2 coupling
127+
# "g3" : 1.4, # SU3 coupling
128+
# "lambda" : 0.13,
129+
# "msq" : -7000 # Units: GeV^2
130+
# }
131+
132+
## @todo kinda would want to have these as individual member variables for easier access.
133+
## But that alone is not good enough as we need to pass the params to other things like the collision module,
134+
## and for that we want some common routine that does not involve hardcoding parameter names. So I anticipate that
135+
## some combination of these approaches would be good.
136+
137+
# end StandardModel
138+
11139

12140
class EffectivePotentialSM(EffectivePotential):
13141

14-
def __init__(self, modelParameters: dict[str, float]):
15-
super().__init__(modelParameters)
142+
def __init__(self, modelParameters: dict[str, float], fieldCount: int):
143+
super().__init__(modelParameters, fieldCount)
16144
## ... do SM specific initialization here. The super call already gave us the model params
17145

18-
19-
def evaluate(self, fields: np.ndarray[float], temperature: float) -> complex:
20-
v = fields[0] # phi ~ 1/sqrt(2) (0, v)
146+
## Count particle degrees-of-freedom to facilitate inclusion of light particle contributions to ideal gas pressure
147+
self.num_boson_dof = 28 #check if correct
148+
self.num_fermion_dof = 90
149+
150+
def evaluate(self, fields: Fields, temperature: float) -> complex:
151+
# phi ~ 1/sqrt(2) (0, v)
152+
v = fields.GetField(0)
153+
21154
T = temperature
22155

23156
# 4D units
@@ -41,13 +174,31 @@ def evaluate(self, fields: np.ndarray[float], temperature: float) -> complex:
41174
# NLO 1-loop correction in Landau gauge, so g^3. Debyes are integrated out by getThermalParameters
42175
V1 = 2*(3-1) * J3(mWsq) + (3-1) * J3(mZsq) + J3(mHsq) + 3.*J3(mGsq)
43176

44-
VTotal = V0 + V1
177+
VTotal = V0 + V1 + self.constantTerms(T)
178+
45179
return VTotal
46-
180+
181+
182+
def constantTerms(self, temperature: npt.ArrayLike) -> npt.ArrayLike:
183+
"""Need to explicitly compute field-independent but T-dependent parts
184+
that we don't already get from field-dependent loops. At leading order in high-T expansion these are just
185+
(minus) the ideal gas pressure of light particles that were not integrated over in the one-loop part.
186+
"""
187+
188+
## See Eq. (39) in hep-ph/0510375 for general LO formula
189+
190+
## How many degrees of freedom we have left. I'm hardcoding the number of DOFs that were done in evaluate(), could be better to pass it from there though
191+
dofsBoson = self.num_boson_dof - 13 # 13 = 6 + 3 + 4 (W + Z + Higgs)
192+
dofsFermion = self.num_fermion_dof - 12 ## we only did top quark loops
193+
194+
## Fermions contribute with a magic 7/8 prefactor as usual. Overall minus sign since Veff(min) = -pressure
195+
return -(dofsBoson + 7./8. * dofsFermion) * np.pi**2 * temperature**4 / 90.
196+
47197

48198
## Calculates thermally corrected parameters to use in Veff. So basically 3D effective params but keeping 4D units
49199
def getThermalParameters(self, temperature: float) -> dict[str, float]:
50200
T = temperature
201+
51202
msq = self.modelParameters["msq"]
52203
lam = self.modelParameters["lambda"]
53204
yt = self.modelParameters["yt"]
@@ -81,112 +232,109 @@ def getThermalParameters(self, temperature: float) -> dict[str, float]:
81232
return thermalParameters
82233

83234

235+
def main():
84236

237+
WallGo.initialize()
85238

86-
## V = msq |phi|^2 + lambda (|phi|^2)^2
87-
class StandardModel(GenericModel):
239+
## Create WallGo control object
240+
manager = WallGoManager()
88241

89-
particles = np.array([], dtype=Particle)
90-
outOfEquilibriumParticles = np.array([], dtype=Particle)
242+
"""Initialize your GenericModel instance.
243+
The constructor currently requires an initial parameter input, but this is likely to change in the future
244+
"""
91245

246+
## QFT model input. Some of these are probably not intended to change, like gauge masses. Could hardcode those directly in the class.
247+
inputParameters = {
248+
"RGScale" : 91.1876,
249+
"v0" : 246.0,
250+
"MW" : 80.379,
251+
"MZ" : 91.1876,
252+
"Mt" : 173.0,
253+
"g3" : 1.2279920495357861,
254+
"mH" : 50.0,
255+
}
92256

93-
def __init__(self):
94257

95-
# Initialize internal Veff with our params dict. @todo will it be annoying to keep these in sync if our params change?
96-
self.effectivePotential = EffectivePotentialSM(self.modelParameters)
258+
model = StandardModel(inputParameters)
97259

98-
## Define particles. this is a lot of clutter, especially if the mass expressions are long,
99-
## so @todo define these in a separate file?
100-
101-
# NB: particle multiplicity is pretty confusing because some internal DOF counting is handled internally already.
102-
# Eg. for SU3 gluons the multiplicity should be 1, NOT Nc^2 - 1.
103-
# But we nevertheless need something like this to avoid having to separately define up, down, charm, strange, bottom
260+
""" Register the model with WallGo. This needs to be done only once.
261+
If you need to use multiple models during a single run, we recommend creating a separate WallGoManager instance for each model.
262+
"""
263+
manager.registerModel(model)
264+
265+
## ---- File name for collisions integrals. Currently we just load this
266+
collisionFileName = pathlib.Path(__file__).parent.resolve() / "Collisions/collisions_top_top_N11.hdf5"
267+
manager.loadCollisionFile(collisionFileName)
268+
269+
## ---- This is where you'd start an input parameter loop if doing parameter-space scans ----
270+
271+
""" Example mass loop that just does one value of mH. Note that the WallGoManager class is NOT thread safe internally,
272+
so it is NOT safe to parallelize this loop eg. with OpenMP. We recommend ``embarrassingly parallel`` runs for large-scale parameter scans.
273+
"""
274+
values_mH = [ 50.0 ]
275+
276+
for mH in values_mH:
277+
278+
inputParameters["mH"] = mH
279+
280+
modelParameters = model.calculateModelParameters(inputParameters)
281+
282+
"""In addition to model parameters, WallGo needs info about the phases at nucleation temperature.
283+
Use the WallGo.PhaseInfo dataclass for this purpose. Transition goes from phase1 to phase2.
284+
"""
285+
286+
Tn = 63.1 ## nucleation temperature
287+
288+
phaseInfo = WallGo.PhaseInfo(temperature = Tn,
289+
phaseLocation1 = WallGo.Fields( [0.0] ),
290+
phaseLocation2 = WallGo.Fields( [246.0] ))
104291

105-
## === Top quark ===
106-
topMsqVacuum = lambda fields: 0.5 * self.modelParameters["yt"]**2 * fields[0]**2
107-
topMsqThermal = lambda T: self.modelParameters["g3"]**2 * T**2 / 6.0
108292

109-
topQuark = Particle("top",
110-
msqVacuum = topMsqVacuum,
111-
msqThermal = topMsqThermal,
112-
statistics = "Fermion",
113-
inEquilibrium = False,
114-
ultrarelativistic = True,
115-
multiplicity = 1
116-
)
117-
self.addParticle(topQuark)
293+
"""Give the input to WallGo. It is NOT enough to change parameters directly in the GenericModel instance because
294+
1) WallGo needs the PhaseInfo
295+
2) WallGoManager.setParameters() does parameter-specific initializations of internal classes
296+
"""
297+
manager.setParameters(modelParameters, phaseInfo)
118298

119-
## === Light quarks, 5 of them ===
120-
lightQuarkMsqThermal = lambda T: self.modelParameters["g3"]**2 * T**2 / 6.0
299+
## TODO initialize collisions. Either do it here or already in registerModel().
300+
## But for now it's just hardcoded in Boltzmann.py and __init__.py
121301

122-
lightQuark = Particle("lightQuark",
123-
msqVacuum = 0.0,
124-
msqThermal = lightQuarkMsqThermal,
125-
statistics = "Fermion",
126-
inEquilibrium = True,
127-
ultrarelativistic = True,
128-
multiplicity = 5
129-
)
130-
self.addParticle(lightQuark)
302+
"""WallGo can now be used to compute wall stuff!"""
131303

132-
## === SU(3) gluon ===
133-
gluonMsqThermal = lambda T: self.modelParameters["g3"]**2 * T**2 * 2.0
304+
## ---- Solve wall speed in Local Thermal Equilibrium approximation
134305

135-
gluon = Particle("gluon",
136-
msqVacuum = 0.0,
137-
msqThermal = gluonMsqThermal,
138-
statistics = "Boson",
139-
inEquilibrium = True,
140-
ultrarelativistic = True,
141-
multiplicity = 1
142-
)
143-
self.addParticle(gluon)
306+
vwLTE = manager.wallSpeedLTE()
144307

308+
print(f"LTE wall speed: {vwLTE}")
145309

146-
## Define parameter dict here. @todo correct values, these are from top of my head.
147-
## In practice the user would define a function that computes these from whatever their input is (mW, mt, mH etc).
148-
## But this still needs to be initialized here as required by the interface. modelParameters = None is fine.
149-
modelParameters = {
150-
"RGScale" : 91, # MS-bar scale. Units: GeV
151-
"yt" : 1.0, # Top Yukawa
152-
"g1" : 0.3, # U1 coupling
153-
"g2" : 0.6, # SU2 coupling
154-
"g3" : 1.4, # SU3 coupling
155-
"lambda" : 0.13,
156-
"msq" : -7000 # Units: GeV^2
157-
}
158-
159-
## @todo kinda would want to have these as individual member variables for easier access.
160-
## But that alone is not good enough as we need to pass the params to other things like the collision module,
161-
## and for that we want some common routine that does not involve hardcoding parameter names. So I anticipate that
162-
## some combination of these approaches would be good.
310+
## ---- Solve field EOM. For illustration, first solve it without any out-of-equilibrium contributions. The resulting wall speed should match the LTE result:
163311

164-
# end StandardModel
312+
## This will contain wall widths and offsets for each classical field. Offsets are relative to the first field, so first offset is always 0
313+
wallParams: WallGo.WallParams
165314

315+
bIncludeOffEq = False
316+
print(f"=== Begin EOM with {bIncludeOffEq=} ===")
166317

318+
wallVelocity, wallParams = manager.solveWall(bIncludeOffEq)
167319

168-
def main():
320+
print(f"{wallVelocity=}")
321+
print(f"{wallParams.widths=}")
322+
print(f"{wallParams.offsets=}")
169323

170-
model = StandardModel()
324+
## Repeat with out-of-equilibrium parts included. This requires solving Boltzmann equations, invoked automatically by solveWall()
325+
bIncludeOffEq = True
326+
print(f"=== Begin EOM with {bIncludeOffEq=} ===")
171327

172-
# test / example
173-
mH = 60
174-
model.modelParameters["msq"] = -mH**2 / 2.
175-
model.modelParameters["lambda"] = 0.5 * mH**2 / 246.**2
328+
wallVelocity, wallParams = manager.solveWall(bIncludeOffEq)
176329

177-
Tn = 200
330+
print(f"{wallVelocity=}")
331+
print(f"{wallParams.widths=}")
332+
print(f"{wallParams.offsets=}")
178333

179-
userInput = {
180-
"Tn" : Tn,
181-
"phaseLocation1" : [ 0.0 ],
182-
"phaseLocation2" : [ 246.0 ]
183-
}
184334

185-
## Create control class
186-
manager = WallGoManager(model, userInput)
335+
# end parameter-space loop
187336

188-
# At this point we should have all required input from the user
189-
# and the manager should have validated it, found phases etc.
337+
# end main()
190338

191339

192340

0 commit comments

Comments
 (0)