unitarity constaint

README.rst

@@ -45,6 +45,9 @@ modification, extension and reuse.
   If you use the W boson pole mass prediction in FlexibleSUSY 2.7.0
   (or later), please cite [2204.05285]_.
 
+  If you use unitarity constraints please cite [XXXX.XXXXX]_ and necessarily
+  [1805.07306_].
+
   FlexibleSUSY depends on SARAH_ and contains components from
   SOFTSUSY_. Therefore, please also cite the following publications
   along with FlexibleSUSY:
@@ -1039,6 +1042,7 @@ References
 .. [1708.05720] `Eur. Phys. J. C77 (2017) no. 12, 814 <https://inspirehep.net/record/1617767>`_ [`arxiv:1708.05720 <https://arxiv.org/abs/1708.05720>`_]
 .. [1710.03760] `CPC 230 (2018) 145-217 <https://inspirehep.net/record/1629978>`_ [`arXiv:1710.03760 <https://arxiv.org/abs/1710.03760>`_]
 .. [1804.09410] `Eur. Phys. J. C78 (2018) no. 7, 573 <https://inspirehep.net/record/1670032>`_ [`arxiv:1804.09410 <https://arxiv.org/abs/1804.09410>`_]
+.. [1805.07306] `Eur. Phys. J. C78 (2018) no. 8, 649 <https://inspirehep.net/literature/1673989`_ [`arxiv:1805.07306 <https://arxiv.org/abs/1805.07306>`_]
 .. [1807.03509] `Eur. Phys. J. C78 (2018) no. 10, 874 <https://inspirehep.net/record/1681658>`_ [`arxiv:1807.03509 <https://arxiv.org/abs/1807.03509>`_]
 .. [1910.03595] `Eur. Phys. J. C80 (2020) no. 3, 186 <https://inspirehep.net/record/1758261>`_ [`arxiv:1910.03595 <https://arxiv.org/abs/1910.03595>`_]
 .. [2106.05038] `CPC 283 (2023) 108584 <https://inspirehep.net/literature/1867840>`_ [`arxiv:2106.05038 <http://arxiv.org/abs/2106.05038>`_]

doc/observables.rst

@@ -24,3 +24,4 @@ List of observables computed by ``FlexibleSUSY``
   - LO order (tree-level or loop-induced) decays of Higgs scalars into non-SM state
   - LO order (tree-level or loop-induced) decays of non-Higgs scalars
 
+- Unitarity costraints

doc/observables/unitarity.rst

@@ -0,0 +1,5 @@
+===========
+Unitarity constraints
+===========
+
+Calculation relies on expression provided by SARAH and requires SARAH >= 4.12.2.

meta/FlexibleSUSY.m

@@ -62,7 +62,8 @@
               "WeinbergAngle`",
               "Wrappers`",
               "Himalaya`",
-              "GM2Calc`"
+              "GM2Calc`",
+              "Unitarity`"
 }];
 
 $flexiblesusyMetaDir     = DirectoryName[FindFile[$Input]];
@@ -2500,6 +2501,18 @@ corresponding tadpole is real or imaginary (only in models with CP
       DeleteDuplicates@Join[vertices,npfVertices]
    ];
 
+WriteUnitarityClass[files_List] :=
+    Module[{},
+      matrix = Unitarity`GetScatteringMatrix[];
+      WriteOut`ReplaceInFiles[files,
+        {"@scatteringPairsLength@" -> ToString@matrix[[2]],
+         "@skipZeros@" -> "if (!(" <> StringRiffle["(i==" <> ToString[First@#-1] <> " && j==" <> ToString[Last@#-1] <> ")"& /@ First@matrix, "||"] <> ")) continue;",
+         "@scatteringElements@" -> TextFormatting`IndentText[Last@matrix],
+         "@generationSizes@" -> ToString[{{#}& /@ matrix[[3]]}],
+      Sequence @@ GeneralReplacementRules[]
+        }];
+    ];
+
 (* Write the AMM c++ files *)
 WriteAMMClass[fields_List, files_List] :=
     Module[{calculation, getMSUSY,
@@ -5201,6 +5214,13 @@ corresponding tadpole is real or imaginary (only in models with CP
                              FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_FFV_form_factors.cpp"}]}}
                ];
 
+           Print["Creating unitarity class..."];
+           WriteUnitarityClass[{{FileNameJoin[{$flexiblesusyTemplateDir, "unitarity.hpp.in"}],
+                               FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_unitarity.hpp"}]},
+                           {FileNameJoin[{$flexiblesusyTemplateDir, "unitarity.cpp.in"}],
+                               FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_unitarity.cpp"}]}}
+           ];
+
            Print["Creating C++ QFT class..."];
            cxxQFTTemplateDir = FileNameJoin[{$flexiblesusyTemplateDir, "cxx_qft"}];
            cxxQFTOutputDir = FileNameJoin[{FSOutputDir, "cxx_qft"}];

meta/Unitarity.m

@@ -0,0 +1,91 @@
+(* :Copyright:
+
+   ====================================================================
+   This file is part of FlexibleSUSY.
+
+   FlexibleSUSY is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published
+   by the Free Software Foundation, either version 3 of the License,
+   or (at your option) any later version.
+
+   FlexibleSUSY is distributed in the hope that it will be useful, but
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with FlexibleSUSY.  If not, see
+   <http://www.gnu.org/licenses/>.
+   ====================================================================
+
+*)
+
+BeginPackage["Unitarity`", {"SARAH`", "Parameters`", "CConversion`", "TreeMasses`"}];
+
+GetScatteringMatrix::usage = "";
+
+Begin["`Private`"];
+
+(* return a C++ lambda computing a given scattering eigenvalue in function of sqrt(s) *)
+ExpressionToCPPLambda[expr__, istates_List, fstates_List] := Module[{params = Parameters`FindAllParametersClassified[expr], paramsCPP, mixingCPP},
+
+   (* CPP definitions of parameters present in the expression *)
+   paramsCPP =
+      StringJoin[
+         ("const auto " <> ToString@CConversion`ToValidCSymbol[#] <>
+            " = model_.get_" <> ToString@CConversion`ToValidCSymbol[#] <> "();\n")& /@ (Parameters`FSModelParameters /. params)
+      ];
+   (* definition of mixing matrices *)
+   mixingCPP = ("auto " <> ToString[#] <> " = [&model_] (int i, int j) { return model_.get_" <> ToString@CConversion`ToValidCSymbol[#] <> "(i-1,j-1); };\n")& /@ (Parameters`FSOutputParameters /. params);
+
+   (* replace input parameters with their FS names *)
+   newExpr = Simplify[expr, Assumptions->s>0] /. Thread[(Parameters`FSModelParameters /. params) -> CConversion`ToValidCSymbol /@ (Parameters`FSModelParameters /. params)];
+   newExpr = newExpr /. Susyno`LieGroups`conj -> Conj;
+   newExpr = newExpr //. SARAH`sum[idx_, start_, stop_, exp_] :> FlexibleSUSY`SUM[idx, start, stop, exp];
+   newExpr = ToString@CForm[newExpr];
+   (* in expressions from SARAH masses are denoted as pmass(X)
+      this converts them to proper C++ form *)
+   newExpr = StringReplace[newExpr, "pmass(" ~~ f:Except[")"].. ~~ "(" ~~ i:DigitCharacter.. ~~ "))" :> "context.mass<" <> f <> ">({" <> ToString[ToExpression[i]-1] <> "})"];
+   newExpr = StringReplace[newExpr, "pmass(" ~~ f:Except[")"].. ~~ "(" ~~ i:Except[")"].. ~~ "))" :> "context.mass<" <> f <> ">({" <> ToString[i] <> "-1})"];
+   newExpr = StringReplace[newExpr, "pmass(" ~~ f:Except[")"]..  ~~ ")" :> "context.mass<" <> f <> ">({})"];
+
+If[expr === 0,
+"[](double sqrtS, int in1, int in2, int out1, int out2) { return 0.; };\n\n",
+"[&model" <> If[!FreeQ[expr, sChan], ",sChan", ""] <> If[!FreeQ[expr, tChan], ",tChan", ""] <> If[!FreeQ[expr, uChan], ",uChan", ""] <> If[!FreeQ[expr, qChan], ",qChan", ""] <> "](double sqrtS, int in1, int in2, int out1, int out2) {
+      auto model_ = model;
+      // couplings should be evaluated at the renormalization scale sqrt(s)
+      // see comment around eq. 20 of 1805.07306
+      model_.run_to(sqrtS);
+      model_.solve_ewsb();
+      const double s = Sqr(sqrtS);
+      const " <> FlexibleSUSY`FSModelName <> "_cxx_diagrams::context_base context {model_};\n" <>
+      "if (sqrtS < context.mass<" <> ToString[fstates[[1]] /. Susyno`LieGroups`conj->Identity] <> ">(" <> If[TreeMasses`GetDimension[fstates[[1]]]>1, "{out1-1}", "{}"] <> ") + context.mass<" <> ToString[fstates[[2]] /. Susyno`LieGroups`conj->Identity] <> ">(" <> If[TreeMasses`GetDimension[fstates[[2]]]>1, "{out2-1}", "{}"] <> ")) return 0.;\n" <>
+      paramsCPP <>
+      mixingCPP <>
+      (* TODO: can this really be complex? *)
+"     return " <> If[FreeQ[expr, Susyno`LieGroups`conj], newExpr, "std::real(" <> newExpr <> ")"] <> ";
+   };\n\n"
+]
+];
+
+GetScatteringMatrix[] := Module[{result, generationSizes},
+   InitUnitarity[];
+   (*RemoveParticlesFromScattering={Se , Sv, Sd, Su};*)
+
+   a0 = Outer[GetScatteringDiagrams[#1 -> #2]&, scatteringPairs, scatteringPairs, 1];
+   generationSizes = Table[{i, j}, {i,1, Length[scatteringPairs]}, {j,1, Length[scatteringPairs]}];
+   generationSizes = Apply[Join[TreeMasses`GetDimension /@ scatteringPairs[[#1]], TreeMasses`GetDimension /@ scatteringPairs[[#2]]]&, generationSizes, {2}];
+
+   result = "";
+   For[i=1, i<=Length[a0], i++,
+      For[j=i, j<=Length[a0[[i]]], j++,
+         result = result <> "// " <> ToString[scatteringPairs[[i]]] <> "->" <> ToString[scatteringPairs[[j]]] <> "\n" <>
+                  "matrix[" <> ToString[i-1] <> "][" <> ToString[j-1] <> "] = " <> ExpressionToCPPLambda[a0[[i,j]], scatteringPairs[[i]], scatteringPairs[[j]]]
+      ]
+   ];
+   {SparseArray[a0]["NonzeroPositions"], Dimensions[a0][[1]], generationSizes, result}
+];
+
+End[];
+EndPackage[];
+

meta/module.mk

@@ -112,6 +112,7 @@ META_SRC     := \
 		$(DIR)/TwoLoopMSSM.m \
 		$(DIR)/TwoLoopQCD.m \
 		$(DIR)/TwoLoopSM.m \
+		$(DIR)/Unitarity.m \
 		$(DIR)/Utils.m \
 		$(DIR)/Vertices.m \
 		$(DIR)/WeinbergAngle.m \

templates/module.mk

@@ -70,6 +70,8 @@ BASE_TEMPLATES := \
 		$(DIR)/susy_parameters.hpp.in \
 		$(DIR)/susy_parameters.cpp.in \
 		$(DIR)/susy_scale_constraint.hpp.in \
+		$(DIR)/unitarity.hpp.in \
+		$(DIR)/unitarity.cpp.in \
 		$(DIR)/utilities.hpp.in \
 		$(DIR)/utilities.cpp.in
 

templates/module.mk.in

@@ -71,6 +71,7 @@ LIB@CLASSNAME@_SRC := \
 		$(DIR)/@CLASSNAME@_slha_io.cpp \
 		$(DIR)/@CLASSNAME@_soft_parameters.cpp \
 		$(DIR)/@CLASSNAME@_susy_parameters.cpp \
+		$(DIR)/@CLASSNAME@_unitarity.cpp \
 		$(DIR)/@CLASSNAME@_utilities.cpp \
 		$(DIR)/@CLASSNAME@_weinberg_angle.cpp
 
@@ -118,6 +119,7 @@ LIB@CLASSNAME@_HDR := \
 		$(DIR)/@CLASSNAME@_soft_parameters.hpp \
 		$(DIR)/@CLASSNAME@_susy_parameters.hpp \
 		$(DIR)/@CLASSNAME@_susy_scale_constraint.hpp \
+		$(DIR)/@CLASSNAME@_unitarity.hpp \
 		$(DIR)/@CLASSNAME@_utilities.hpp \
 		$(DIR)/@CLASSNAME@_weinberg_angle.hpp
 

templates/run.cpp.in

@@ -25,6 +25,7 @@
 #include "@ModelName@_spectrum_generator.hpp"
 #include "@ModelName@_utilities.hpp"
 @decaysIncludes@
+#include "@ModelName@_unitarity.hpp"
 
 @solverIncludes@
 #include "physical_input.hpp"
@@ -97,6 +98,7 @@ int run_solver(flexiblesusy::@ModelName@_slha_io& slha_io,
 
 @calculateDecaysForModel@
 
+   @ModelName@_unitarity::max_scattering_eigenvalue(std::get<0>(models));
    const bool show_result = !problems.have_problem() ||
       spectrum_generator_settings.get(Spectrum_generator_settings::force_output);
    // SLHA output

templates/unitarity.cpp.in

@@ -0,0 +1,157 @@
+// ====================================================================
+// This file is part of FlexibleSUSY.
+//
+// FlexibleSUSY is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License,
+// or (at your option) any later version.
+//
+// FlexibleSUSY is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with FlexibleSUSY.  If not, see
+// <http://www.gnu.org/licenses/>.
+// ====================================================================
+
+
+/**
+ * @file @ModelName@_unitarity.cpp
+ *
+ * This file was generated with FlexibleSUSY @FlexibleSUSYVersion@ and SARAH @SARAHVersion@ .
+ */
+
+#include "@ModelName@_unitarity.hpp"
+#include "@ModelName@_mass_eigenstates.hpp"
+#include "@ModelName@_context_base.hpp"
+
+#include "minimizer.hpp"
+#include "wrappers.hpp"
+#include "sum.hpp"
+
+#include <gsl/gsl_min.h>
+#include <algorithm>
+#include <array>
+#include <functional>
+
+namespace flexiblesusy {
+namespace @ModelName@_unitarity {
+
+inline double Sqrt2(int i, int j) {
+   // 1/Sqrt[2] or 1
+   return i==j ? 0.707106781186547524 : 1;
+}
+
+double max_scattering_eigenvalue(const @ModelName@_mass_eigenstates& model) {
+   constexpr size_t size = @scatteringPairsLength@;
+   std::array<std::array<std::function<double(double, int, int, int, int)>, size>, size> matrix = {};
+   std::array<std::function<double(double)>, size> row_temp = {};
+   const double sChan = 1.;
+   const double tChan = 1.;
+   const double uChan = 1.;
+   const double qChan = 1.;
+
+   using namespace @ModelName@_cxx_diagrams::fields;
+
+@scatteringElements@
+   std::array<std::array<double, size>, size> max_eigenvalues = {};
+
+   constexpr std::array<std::array<std::array<int, 4>, size>, size> generationSizes =
+      @generationSizes@;
+
+   for (int i=0; i<size; i++) {
+      for (int j=i; j<size; j++) {
+         @skipZeros@
+         std::array<int, 4>  externalGenerationSizes = generationSizes[i][j];
+         // suply indices in mathematica counting (starting from 1)
+         for (int iIn1 = 1; iIn1 <= externalGenerationSizes[0]; ++iIn1) {
+         for (int iIn2 = 1; iIn2 <= externalGenerationSizes[1]; ++iIn2) {
+         for (int iOut1 = 1; iOut1 <= externalGenerationSizes[2]; ++iOut1) {
+         for (int iOut2 = 1; iOut2 <= externalGenerationSizes[3]; ++iOut2) {
+         std::function<double(double)> f = std::bind(matrix[i][j], std::placeholders::_1, iIn1, iIn2, iOut1, iOut2);
+
+         int status;
+         int iter = 0, max_iter = 100;
+         const gsl_min_fminimizer_type *T;
+         gsl_min_fminimizer *s;
+         double a = 100.0, b = 5000.0;
+
+         static constexpr int grid_size = 10;
+         std::array<double, grid_size> grid;
+         for (int idx = 0; idx < grid_size; ++idx) { grid[idx] = a + (idx+1)*(b-a)/(grid.size()+1); }
+         std::array<double, grid_size> f_grid;
+         for (int idx = 0; idx < grid_size; ++idx) { f_grid[idx] = -std::abs(f(grid[idx])); }
+         double m = grid[std::distance(f_grid.begin(), std::min_element(f_grid.begin(), f_grid.end()))];
+         // the minimum is between [a,b]
+         if (-std::abs(f(m)) < std::min(-std::abs(f(a)), -std::abs(f(b)))) {
+         gsl_function F = {
+            [](double d, void* vf) -> double {
+               auto& f = *static_cast<std::function<double(double)>*>(vf);
+               return -std::abs(f(d));
+            },
+            &f
+         };
+         T = gsl_min_fminimizer_brent;
+         s = gsl_min_fminimizer_alloc (T);
+         gsl_min_fminimizer_set (s, &F, m, a, b);
+
+         printf ("using %s method\n",
+            gsl_min_fminimizer_name (s));
+
+         printf ("%5s [%9s, %9s] %9s %10s %9s\n",
+            "iter", "lower", "upper", "min",
+            "err", "err(est)");
+
+         printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
+            iter, a, b,
+            m, b - a);
+
+         do
+         {
+           iter++;
+           status = gsl_min_fminimizer_iterate (s);
+
+           m = gsl_min_fminimizer_x_minimum (s);
+           a = gsl_min_fminimizer_x_lower (s);
+           b = gsl_min_fminimizer_x_upper (s);
+
+          // determine sqrt(s) with a relative error < 0.1%
+          status
+             = gsl_min_test_interval (a, b, 0.0, 1e-3);
+
+          if (status == GSL_SUCCESS)
+             printf ("Converged:\n");
+
+             printf ("%5d [%.7f, %.7f] "
+                "%.7f %.7f\n",
+                iter, a, b,
+                m, b - a);
+         }
+         while (status == GSL_CONTINUE && iter < max_iter);
+         max_eigenvalues[i][j] = std::abs(f(m));
+
+         gsl_min_fminimizer_free (s);
+         }
+         // the minimum is at a or b
+         else {
+            max_eigenvalues[i][j] = std::max(std::abs(f(m)), max_eigenvalues[i][j]);
+         }
+      }
+   }
+         }}}}
+   std::cout << "=======================================\n";
+   double max_elem = 0.;
+   for (int i=0; i<size; i++) {
+      for (int j=i; j<size; j++) {
+         std::cout << max_eigenvalues[i][j] << std::endl;
+         max_elem = std::max(max_elem, max_eigenvalues[i][j]);
+      }
+   }
+   std::cout << max_elem << std::endl;
+   return max_elem;
+}
+
+} // namespace @ModelName@_unitarity
+} // namespace flexiblesusy