Fix tail length (#351)

* Fixed the tail lengths in grid3Scales

* Changing input meanFreePathScales to half of the original

* EOM minimum tails to scale with new code

* One more place for scaling

* Fixed a small error in the doc of grid3Scales

---------

Co-authored-by: Benoit <benoit.laurent96@gmail.com>
Co-authored-by: Joonas Hirvonen <joonas.hirvonen@nottingham.ac.uk>

Author: 3 people

Models/InertDoubletModel/inertDoubletModel.py

@@ -861,7 +861,7 @@ 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
+                    meanFreePathScale=50.0, # In units of 1/Tnucl
                     wallThicknessGuess=5.0, # In units of 1/Tnucl
                 ),
             )

Models/ManySinglets/manySinglets.py

@@ -742,7 +742,7 @@ 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
+                    meanFreePathScale=50.0,  # In units of 1/Tnucl
                     wallThicknessGuess=5.0,  # In units of 1/Tnucl
                 ),
             )

Models/SingletStandardModel_Z2/singletStandardModelZ2.py

@@ -711,7 +711,7 @@ 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
+                    meanFreePathScale=50.0, # In units of 1/Tnucl
                     wallThicknessGuess=5.0, # In units of 1/Tnucl
                 ),
             )

Models/StandardModel/standardModel.py

@@ -564,7 +564,7 @@ def getBenchmarkPoints(self) -> list[ExampleInputPoint]:
                     WallGo.WallSolverSettings(
                         # we actually do both cases in the common example
                         bIncludeOffEquilibrium=True,
-                        meanFreePathScale=200.0, # In units of 1/Tnucl
+                        meanFreePathScale=100.0, # In units of 1/Tnucl
                         wallThicknessGuess=20.0, # In units of 1/Tnucl
                     ),
                 )

Models/Yukawa/yukawa.py

@@ -205,7 +205,7 @@ def main() -> None:
         # meanFreePathScale is determined here by the annihilation channels,
         # and scales inversely with y^4 or lam^2. This is why
         # meanFreePathScale has to be so large.
-        meanFreePathScale=10000.0,  # In units of 1/Tnucl
+        meanFreePathScale=5000.0,  # In units of 1/Tnucl
         wallThicknessGuess=10.0,  # In units of 1/Tnucl
     )
 

Models/Yukawa/yukawaWithCollisionGeneration.py

@@ -202,7 +202,7 @@ def getBenchmarkPoints(self) -> list[ExampleInputPoint]:
                     # meanFreePathScale is determined here by the annihilation channels,
                     # and scales inversely with y^4 or lam^2. This is why
                     # meanFreePathScale has to be so large.
-                    meanFreePathScale=10000.0,  # In units of 1/Tnucl
+                    meanFreePathScale=5000.0,  # In units of 1/Tnucl
                     wallThicknessGuess=10.0,  # In units of 1/Tnucl
                 ),
             )

src/WallGo/equationOfMotion.py

@@ -1251,17 +1251,17 @@ def _updateGrid(self, wallParams: WallParams, velocityMid: float) -> None:
         gammaWall = 1 / np.sqrt(1 - velocityMid**2)
         """ The length of the tail inside typically scales like gamma, while the one
         outside like 1/gamma. We take the max because the tail lengths must be larger
-        than wallThicknessGrid*(1+2*smoothing)/ratioPointsWall """
+        than wallThicknessGrid*(1/2+smoothing)/ratioPointsWall """
         tailInside = max(
             self.meanFreePathScale * gammaWall * self.includeOffEq,
             wallThicknessGrid
-            * (1 + 2.1 * self.grid.smoothing)
+            * (0.5 + 1.05 * self.grid.smoothing)
             / self.grid.ratioPointsWall,
         )
         tailOutside = max(
             self.meanFreePathScale / gammaWall * self.includeOffEq,
             wallThicknessGrid
-            * (1 + 2.1 * self.grid.smoothing)
+            * (0.5 + 1.05 * self.grid.smoothing)
             / self.grid.ratioPointsWall,
         )
         self.grid.changePositionFalloffScale(

src/WallGo/grid3Scales.py

@@ -24,7 +24,7 @@ class Grid3Scales(Grid):
     thickness, so that
 
     .. math::
-        z'(\chi) \approx \frac{L}{r}\chi, \quad \chi \in [-r, r].
+        z'(\chi) \approx \frac{L}{r}, \quad \chi \in [-r, r].
 
     It is easier to find the derivative of a function that has these properties,
     and then integrate it. We choose here
@@ -36,9 +36,9 @@ class Grid3Scales(Grid):
     
     .. math::
         f(\chi) \approx \begin{cases}
-            \lambda_-,& \chi<-r,\\
+            2\lambda_-,& \chi<-r,\\
             L/r,& -r<\chi<r,\\
-            \lambda_+,& \chi>r.
+            2\lambda_+,& \chi>r.
         \end{cases}
     
     We choose :math:`f(\chi)` to be a sum of functions like
@@ -78,11 +78,11 @@ def __init__(
             (and :math:`\rho_z` and :math:`\rho_\Vert`) directions.
         tailLengthInside : float
             Decay length of the solution's tail inside the wall. Should be larger
-            than wallThickness*(1+2*smoothing)/ratioPointsWall. Should be
+            than wallThickness*(1/2+smoothing)/ratioPointsWall. Should be
             expressed in physical units (the units used in EffectivePotential).
         tailLengthOutside : float
             Decay length of the solution's tail outside the wall. Should be larger
-            than wallThickness*(1+2*smoothing)/ratioPointsWall. Should be
+            than wallThickness*(1/2+smoothing)/ratioPointsWall. Should be
             expressed in physical units (the units used in EffectivePotential).
         wallThickness : float
             Thickness of the wall. Should be expressed in physical units
@@ -154,11 +154,11 @@ def _updateParameters(
         assert wallThickness > 0, "Grid3Scales error: wallThickness must be positive."
         assert smoothing > 0, "Grid3Scales error: smoothness must be positive."
         assert (
-            tailLengthInside > wallThickness * (1 + 2 * smoothing) / ratioPointsWall
+            tailLengthInside > wallThickness * (1 / 2 + smoothing) / ratioPointsWall
         ), """Grid3Scales error: tailLengthInside must be greater than
         wallThickness*(1+2*smoothness)/ratioPointsWall."""
         assert (
-            tailLengthOutside > wallThickness * (1 + 2 * smoothing) / ratioPointsWall
+            tailLengthOutside > wallThickness * (1 / 2 + smoothing) / ratioPointsWall
         ), """Grid3Scales error: tailLengthOutside must be greater than
         wallThickness*(1+2*smoothness)/ratioPointsWall."""
         assert (
@@ -181,18 +181,18 @@ def _updateParameters(
             * smoothing
             * wallThickness
             * ratioPointsWall**2
-            * (ratioPointsWall * tailLengthInside - wallThickness * (1 + smoothing))
+            * (2 * ratioPointsWall * tailLengthInside - wallThickness * (1 + smoothing))
         ) / abs(
-            ratioPointsWall * tailLengthInside - wallThickness * (1 + 2 * smoothing)
+            2 * ratioPointsWall * tailLengthInside - wallThickness * (1 + 2 * smoothing)
         )
         self.aOut = np.sqrt(
             4
             * smoothing
             * wallThickness
             * ratioPointsWall**2
-            * (ratioPointsWall * tailLengthOutside - wallThickness * (1 + smoothing))
+            * (2 * ratioPointsWall * tailLengthOutside - wallThickness * (1 + smoothing))
         ) / abs(
-            ratioPointsWall * tailLengthOutside - wallThickness * (1 + 2 * smoothing)
+            2 * ratioPointsWall * tailLengthOutside - wallThickness * (1 + 2 * smoothing)
         )
 
     def decompactify(
@@ -213,7 +213,7 @@ def decompactify(
 
         def term1(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
             return np.array((1 - r)
-                * (r * tailOut - L)
+                * (2 * r * tailOut - L)
                 * np.arctanh(
                     (1 - x + np.sqrt(aOut**2 + (x - r) ** 2))
                     / np.sqrt(aOut**2 + (1 - r) ** 2)
@@ -225,7 +225,7 @@ def term1(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
 
         def term2(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
             return np.array(-(1 + r)
-                * (r * tailOut - L)
+                * (2 * r * tailOut - L)
                 * np.arctanh(
                     (1 + x - np.sqrt(aOut**2 + (x - r) ** 2))
                     / np.sqrt(aOut**2 + (1 + r) ** 2)
@@ -236,7 +236,7 @@ def term2(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
             )
         def term3(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
             return np.array((1 - r)
-                * (r * tailIn - L)
+                * (2 * r * tailIn - L)
                 * np.arctanh(
                     (1 + x - np.sqrt(aIn**2 + (x + r) ** 2))
                     / np.sqrt(aIn**2 + (1 - r) ** 2)
@@ -247,7 +247,7 @@ def term3(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
             )
         def term4(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
             return np.array(-(1 + r)
-                * (r * tailIn - L)
+                * (2 * r * tailIn - L)
                 * np.arctanh(
                     (1 - x + np.sqrt(aIn**2 + (x + r) ** 2))
                     / np.sqrt(aIn**2 + (1 + r) ** 2)
@@ -257,7 +257,7 @@ def term4(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
                 / r
             )
         def term5(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
-            return np.array((tailIn + tailOut - 4 * self.smoothing * L / r)
+            return np.array((2 * tailIn + 2 * tailOut - 4 * self.smoothing * L / r)
                             * np.arctanh(x)
             )
         def totalMapping(x: np.ndarray) -> np.ndarray: # pylint: disable=invalid-name
@@ -289,12 +289,12 @@ def compactificationDerivatives(
         aOut = self.aOut
 
         dzdzCompact = (
-            (tailIn - L / r)
+            (2 * tailIn - L / r)
             * (1 - (zCompact + r) / np.sqrt(aIn**2 + (zCompact + r) ** 2))
             / 2
         )
         dzdzCompact += (
-            (tailOut - L / r)
+            (2 * tailOut - L / r)
             * (1 + (zCompact - r) / np.sqrt(aOut**2 + (zCompact - r) ** 2))
             / 2
         )

src/WallGo/manager.py

@@ -33,7 +33,7 @@ class WallSolverSettings:
     """If False, will ignore all out-of-equilibrium effects (no Boltzmann solving).
     """
 
-    meanFreePathScale: float = 100.0
+    meanFreePathScale: float = 50.0
     """Estimate of the mean free path of the plasma in units of 1/Tn. This will be used
     to set the tail lengths in the Grid object. Default is 100.
     """
@@ -584,7 +584,7 @@ def buildGrid(
 
         # We divide by Tnucl to get it in physical units of length
         tailLength = max(
-            meanFreePathScale, wallThicknessIni * (1.0 + 3.0 * smoothing) / ratioPointsWall
+            meanFreePathScale, 0.5 * wallThicknessIni * (1.0 + 3.0 * smoothing) / ratioPointsWall
         ) / Tnucl
 
         if gridN % 2 == 0: