Changes made on the private repo for uncertainties paper

CHANGELOG.rst

@@ -2,6 +2,20 @@
 Changelog
 ==================
 
+1.1.0 (2025-08-20)
+==================
+
+* Phase tracing updates to improve algorithm stability.
+* User may now pass the free energy of a phase as arrays.
+* Spectral truncation options added to avoid aliasing, with an automatic method used by default.
+* Error estimate for the wall speed updated to make it more theoretically sound.
+* Linearization criterion updated so that `linearizationCriterion2` estimates the second order correction.
+* Error estimate due to the Tanh ansatz now computed and added to `WallGoResults`.
+* Accuracy of energy-momentum conservation computed, and printed for logging level `DEBUG`.
+* New tests added for the Standard Model with light Higgs.
+* `PTTools <https://github.com/CFT-HY/pttools>`_ example file added.
+
+
 1.0.0 (2024-11-07)
 ==================
 

Models/InertDoubletModel/inertDoubletModel.py

@@ -378,7 +378,7 @@ def jCW(
 
         # Note that we are taking the absolute value of the mass in the log here,
         # instead of using EImaginaryOption = ABS_ARGUMENT, because we do not 
-        # want the absolute value in the product of massSq ans rgScale
+        # want the absolute value in the product of massSq and rgScale
         return degreesOfFreedom*np.array( 
             massSq * massSq * (np.log(np.abs(massSq / rgScale**2)) - c)
             + 2 * massSq * rgScale**2

Models/InertDoubletModel/inertDoubletModelConfig.ini

@@ -83,6 +83,9 @@ phaseTracerTol = 1e-8
 # First step size in units of the maximum step size. Use None for default algorithm.
 phaseTracerFirstStep = None
 
+# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives.
+interpolationDegree = 1
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility
@@ -96,3 +99,6 @@ basisM = Cardinal
 
 # The momentum polynomial basis type, either Cardinal or Chebyshev.
 basisN = Chebyshev
+
+# Whether or not to truncate the spectral expansion
+truncationOption = AUTO

Models/ManySinglets/manySingletsConfig.ini

@@ -83,6 +83,9 @@ phaseTracerTol = 1e-6
 # First step size in units of the maximum step size. Use None for default algorithm.
 phaseTracerFirstStep = None
 
+# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives.
+interpolationDegree = 1
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility
@@ -96,3 +99,6 @@ basisM = Cardinal
 
 # The momentum polynomial basis type, either Cardinal or Chebyshev.
 basisN = Chebyshev
+
+# Whether or not to truncate the spectral expansion
+truncationOption = AUTO

Models/SingletStandardModel_Z2/exampleCollisionDefs.py

@@ -1,3 +1,10 @@
+"""
+This Python script, exampleCollisionDefs.py,
+uses the implementation of the minimal Standard Model
+extension in singletStandardModelZ2.py and adds example
+definitions for the Boltzmann collision integrals for
+WallGoCollision.
+"""
 import WallGoCollision
 import pathlib
 import sys

Models/SingletStandardModel_Z2/exampleOutputThermodynamics.py

@@ -0,0 +1,129 @@
+"""
+This Python script, exampleOutputThermodynamics.py,
+uses the implementation of the minimal Standard Model
+extension in singletStandardModelZ2.py and gives methods
+for saving the thermodynamics of the model for later use
+in e.g. PTTools.
+"""
+
+import sys
+import pathlib
+import numpy as np
+import h5py
+
+# WallGo imports
+import WallGo  # Whole package, in particular we get WallGo._initializeInternal()
+
+# Add the SingletStandardModelZ2 folder to the path to import SingletSMZ2
+sys.path.append(str(pathlib.Path(__file__).resolve().parent))
+from singletStandardModelZ2 import SingletSMZ2
+
+
+def main() -> None:
+
+    manager = WallGo.WallGoManager()
+
+    # Model definition is done in the SingletSMZ2 class
+    model = SingletSMZ2()
+    manager.registerModel(model)
+
+    inputParameters = {
+        "RGScale": 125.0,
+        "v0": 246.0,
+        "MW": 80.379,
+        "MZ": 91.1876,
+        "Mt": 173.0,
+        "g3": 1.2279920495357861,
+        "mh1": 125.0,
+        "mh2": 120.0,
+        "lHS": 0.9,
+        "lSS": 1.0,
+    }
+
+    model.updateModel(inputParameters)
+
+    # Creates a Thermodynamics object, containing all thermodynamic functions
+    # by tracing the high- and low-temperature phases.
+    # For temperatures outside of the range of existence of the phases, the
+    # thermodynamic functions are extrapolated using the template model.
+    manager.setupThermodynamicsHydrodynamics(
+        WallGo.PhaseInfo(
+            temperature=100.0,  # nucleation temperature
+            phaseLocation1=WallGo.Fields([0.0, 200.0]),
+            phaseLocation2=WallGo.Fields([246.0, 0.0]),
+        ),
+        WallGo.VeffDerivativeSettings(
+            temperatureVariationScale=10.0, fieldValueVariationScale=[10.0, 10.0]
+        ),
+    )
+
+    Tcrit = manager.thermodynamics.findCriticalTemperature(dT = 0.1)
+
+    # Generate tables for the thermodynamics functions of both phases.
+    # Note that if the minimum or maximum temperature chosen is outside of
+    # the range of existence of the phases, the extrapolation is used.
+
+    # Temperature range and step for high-temperature phase
+    highTPhaseRange = (80.0, 140.0, 0.1)
+    # Temperature range and step for low-temperature phase
+    lowTPhaseRange = (80.0, 110.0, 0.1)
+
+    # Create temperature arrays
+    temp_high = np.arange(highTPhaseRange[0], highTPhaseRange[1], highTPhaseRange[2])
+    temp_low = np.arange(lowTPhaseRange[0], lowTPhaseRange[1], lowTPhaseRange[2])
+
+    # Evaluate thermodynamic functions for high-temperature phase
+    p_high = np.array([manager.thermodynamics.pHighT(T) for T in temp_high])
+    e_high = np.array([manager.thermodynamics.eHighT(T) for T in temp_high])
+    csq_high = np.array([manager.thermodynamics.csqHighT(T) for T in temp_high])
+
+    # Evaluate thermodynamic functions for low-temperature phase
+    p_low = np.array([manager.thermodynamics.pLowT(T) for T in temp_low])
+    e_low = np.array([manager.thermodynamics.eLowT(T) for T in temp_low])
+    csq_low = np.array([manager.thermodynamics.csqLowT(T) for T in temp_low])
+
+    # Get temperature limits; outside of these limits, extrapolation has been used
+    max_temp_high = manager.thermodynamics.freeEnergyHigh.maxPossibleTemperature[0]
+    min_temp_high = manager.thermodynamics.freeEnergyHigh.minPossibleTemperature[0]
+    max_temp_low = manager.thermodynamics.freeEnergyLow.maxPossibleTemperature[0]
+    min_temp_low = manager.thermodynamics.freeEnergyLow.minPossibleTemperature[0]
+
+    # Create HDF5 file
+    filename = pathlib.Path(__file__).resolve().parent/"thermodynamics_data.h5"
+
+    with h5py.File(filename, 'w') as f:
+        # Global attributes
+        f.attrs["model_label"] = "singletSMZ2"
+        f.attrs["critical_temperature"] = Tcrit
+        f.attrs["nucleation_temperature"] = 100.
+
+        # High-temperature phase group
+        high_group = f.create_group("high_temperature_phase")
+        high_group.create_dataset("temperature", data=temp_high)
+        high_group.create_dataset("pressure", data=p_high)
+        high_group.create_dataset("energy_density", data=e_high)
+        high_group.create_dataset("sound_speed_squared", data=csq_high)
+        high_group.attrs["max_possible_temperature"] = max_temp_high
+        high_group.attrs["min_possible_temperature"] = min_temp_high
+
+        # Low-temperature phase group
+        low_group = f.create_group("low_temperature_phase")
+        low_group.create_dataset("temperature", data=temp_low)
+        low_group.create_dataset("pressure", data=p_low)
+        low_group.create_dataset("energy_density", data=e_low)
+        low_group.create_dataset("sound_speed_squared", data=csq_low)
+        low_group.attrs["max_possible_temperature"] = max_temp_low
+        low_group.attrs["min_possible_temperature"] = min_temp_low
+
+    print(f"Thermodynamics data saved to {filename}")
+    print("Model: singletSMZ2")
+    print(f"Critical temperature: {Tcrit:.2f} GeV")
+    print(f"High-T phase: {len(temp_high)} temperature points from {temp_high[0]:.1f} to {temp_high[-1]:.1f} GeV")
+    print(f"Low-T phase: {len(temp_low)} temperature points from {temp_low[0]:.1f} to {temp_low[-1]:.1f} GeV")
+    print(f"High-T phase valid range: {min_temp_high:.1f} to {max_temp_high:.1f} GeV")
+    print(f"Low-T phase valid range: {min_temp_low:.1f} to {max_temp_low:.1f} GeV")
+
+
+## Don't run the main function if imported to another file
+if __name__ == "__main__":
+    main()

Models/SingletStandardModel_Z2/singletStandardModelZ2Config.ini

@@ -83,6 +83,9 @@ phaseTracerTol = 1e-6
 # First step size in units of the maximum step size. Use None for default algorithm.
 phaseTracerFirstStep = None
 
+# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives.
+interpolationDegree = 1
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility
@@ -96,3 +99,6 @@ basisM = Cardinal
 
 # The momentum polynomial basis type, either Cardinal or Chebyshev.
 basisN = Chebyshev
+
+# Whether or not to truncate the spectral expansion
+truncationOption = AUTO

Models/StandardModel/standardModel.py

@@ -25,7 +25,6 @@
 A Study of the electroweak phase transition dynamics, Phys.Rev.D 52 (1995) 7182-7204
 doi:10.1103/PhysRevD.52.7182
 """
-
 import sys
 import pathlib
 import numpy as np
@@ -299,9 +298,9 @@ def evaluate(  # pylint: disable=R0914
         lambdaT = self.modelParameters["lambda"] - 3 / (
             16 * np.pi * np.pi * self.modelParameters["v0"] ** 4
         ) * (
-            2 * mW**4 * np.log(mW**2 / (ab * T**2) + 1e-100)
-            + mZ**4 * np.log(mZ**2 / (ab * T**2) + 1e-100)
-            - 4 * mt**4 * np.log(mt**2 / (af * T**2) + 1e-100)
+            2 * mW**4 * np.log(mW**2 / (ab * T**2) )
+            + mZ**4 * np.log(mZ**2 / (ab * T**2) )
+            - 4 * mt**4 * np.log(mt**2 / (af * T**2) )
         )
 
         cT: float | np.ndarray = self.modelParameters["C0"] + 1 / (
@@ -317,7 +316,7 @@ def evaluate(  # pylint: disable=R0914
 
         potentialT: float | np.ndarray = (
             self.modelParameters["D"] * (T**2 - self.modelParameters["T0sq"]) * v**2
-            - cT * T**2 * pow(v, 2) * np.log(np.abs(v / T))
+            - cT * T**2 * pow(v, 2) * np.log(np.abs(v / T) + 1e-100) # Avoid log(0)
             - eT * T * pow(v, 3)
             + lambdaT / 4 * pow(v, 4)
         )
@@ -529,7 +528,7 @@ def getBenchmarkPoints(self) -> list[ExampleInputPoint]:
         of the Higgs mass.
         """
         valuesMH = [0.0, 34.0, 50.0, 70.0, 81.0]
-        valuesTn = [57.192, 70.579, 83.426, 102.344, 113.575]
+        valuesTn = [57.1958, 70.5793, 83.4251, 102.344, 113.575]
 
         output: list[ExampleInputPoint] = []
 
@@ -555,8 +554,8 @@ def getBenchmarkPoints(self) -> list[ExampleInputPoint]:
                     WallGo.WallSolverSettings(
                         # we actually do both cases in the common example
                         bIncludeOffEquilibrium=True,
-                        meanFreePathScale=100.0, # In units of 1/Tnucl
-                        wallThicknessGuess=20.0, # In units of 1/Tnucl
+                        meanFreePathScale=50.0, # In units of 1/Tnucl
+                        wallThicknessGuess=5.0, # In units of 1/Tnucl
                     ),
                 )
             )

Models/StandardModel/standardModelConfig.ini

@@ -16,7 +16,7 @@ momentumGridSize = 11
 ratioPointsWall = 0.5
 
 # Smoothing factor of the mapping function (the larger the smoother)
-smoothing = 0.1
+smoothing = 0.3
 
 [EquationOfMotion]
 # The absolute error tolerance for the wall-velocity result
@@ -83,6 +83,9 @@ phaseTracerTol = 1e-6
 # First step size in units of the maximum step size. Use None for default algorithm.
 phaseTracerFirstStep = None
 
+# Degree of the splines used in FreeEnergy to interpolate the potential and its derivatives.
+interpolationDegree = 1
+
 [BoltzmannSolver]
 # Factor multiplying the collision term in the Boltzmann equation.
 # Can be used for testing or for studying the solution's sensibility
@@ -96,3 +99,6 @@ basisM = Cardinal
 
 # The momentum polynomial basis type, either Cardinal or Chebyshev.
 basisN = Chebyshev
+
+# Whether or not to truncate the spectral expansion
+truncationOption = AUTO