""" Module providing a class to go from matrix elements squared,
statistics and thermal masses to the numerical evaluation of the
rate with necessary HTL resummations """
# This software is covered by the GNU General Public License 3.
# Copyright (C) 2024-2025 Killian Bouzoud, Jacopo Ghiglieri, Greg Jackson
#use this for typing
from typing import Tuple, Union, Optional
#handle warnings
import warnings
#start by importing sympy and other needed modules
import sympy
import numpy
from scipy import LowLevelCallable
from scipy.integrate import nquad, IntegrationWarning
import integrator # type: ignore
#create low-level scipy functions from the cython functions
#these are slightly faster for numerics
t4fast=LowLevelCallable.from_cython(integrator,'dim4tvec')
s4fast=LowLevelCallable.from_cython(integrator,'dim4svec')
t5fast=LowLevelCallable.from_cython(integrator,'dim5tvec')
s5fast=LowLevelCallable.from_cython(integrator,'dim5svec')
[docs]
class AutothermWarning(Warning):
"""Empty placeholder class to be able to customize warnings
occurring during numerical integration"""
#we want these warning to be always displayed (can of course be turned off by the user)
warnings.filterwarnings("always", category=AutothermWarning)
[docs]
def decompose(msq: str,stats: Tuple[int],couplings: Tuple[str]) -> dict:
r"""
Decomposes a string matrix element squared ``msq`` with statistics
``stats``, couplings ``couplings`` into the coefficients :math:`a_1`,
:math:`a_2`, :math:`b_1` (and :math:`b_2` for dimension 5). In the
case of higher-dimensional couplings, the function returns :math:`a_t`
and :math:`b_s`, where the former contains all pieces of the original matrix
element squared that feature positive powers of :math:`t` at the denominator
and the latter is the remainder.
``stats`` should be a tuple (t1,s1,s2) for the statistics of the other outgoing
particle t1 and the two incoming ones s1 and s2. +1 for a boson, -1 for a fermion.
``couplings`` should be a tuple of strings. Interactions of dimension different
from 4 or 5 are supported only if IR finite.
:param str msq: A string for a matrix element squared
:param tuple stats: A tuple (t1,s1,s2) of integers
:param tuple couplings: A tuple of strings for the couplings
:return: A dict with stats as keys and a list of [a1,a2,b1], [a1,a2,b1,b2] or
[at,bs] as values, depending on the dimensionality of the interaction
:rtype: dict
Examples
========
.. code-block:: python
from numerical.manipulate import *
decompose("kappa^2*g1**2*(s*t/u + s*u/t + t*u/s)",(1,1,1),("kappa","g1"))
"""
#warning: input is not fully sanitized:
# errors might not be thrown if msq and couplings are not str objects
#first of all add mandelstam variables and couplings to the list of known symbols
nslist={}
for coupling in couplings:
nslist[coupling] = sympy.Symbol(coupling)
s = sympy.Symbol("s")
t = sympy.Symbol("t")
u = sympy.Symbol("u")
for mandelstam in ["s","t","u"]:
nslist[mandelstam] = sympy.Symbol(mandelstam)
#create a list of the Mandelstam variables and of the type of denominators we expect
mandvars=[s,t,u]
wanteddenons=[sympy.S("1"),s,t,u,s**2,t**2,u**2]
#create a symbolic version of the matrix element squared with sympy
#expects a matrix element squared of the form \sum_i e^a g_i^b f_i(s,t,u)
# with g_i the couplings of the bath particles and e the bath-relic
#coupling. I guess it would be easy to generalize to multiple such
# couplings at once. And also to account for some
# masses for the bath particles??
msq_symb_start=sympy.sympify(msq,strict=False,locals=nslist)
#handle directly the case in which the expression trivially simplifies
# to zero (vanishing contribution in a flavor index projection, for instance)
if msq_symb_start == 0 :
return {}
#extract the dimensionality of the matrix elements squared
#to do so, rescale s, t and u by a dimensionful constant and see how the expression scales in it
energyscale=sympy.S("energyscale")
msq_symb_scale=msq_symb_start.subs({s: energyscale*s,t: energyscale*t,u:energyscale*u})
dim = 4+sympy.simplify(sympy.Derivative(msq_symb_scale,energyscale)*energyscale/msq_symb_scale)
if not isinstance(dim,sympy.core.numbers.Integer):
raise RuntimeError(f"could not work out the dimensionality {dim} "
f"of the interaction for the expression {msq}")
else:
dim=int(dim)
#fail if unsupported dimensionality of the interaction operator
fallback=False
if dim !=4 and dim != 5:
warnings.warn("Interaction of unsupported dimensionality, "
"using C fallback mode",AutothermWarning,stacklevel=2)
fallback=True
#now prepare to do some manipulations in case the matrix element is not
# in a clean form such as (s u/t+tu/s+st/u) but in some convoluted one like
#(s^2+st+t^2)^2/(stu)
#to do that, introduce some control variables for a while loop
fixed_denominators=False
niters=0
#we want to avoid an infinite loop
nitersmax=20
#expand out the matrix element squared
msq_symb=msq_symb_start.expand()
while (not fixed_denominators and niters<nitersmax):
newexpr=0
fixed_denominators=True
#the way to proceed is to decompose the matrix element into its basis elements
#with this we separate the matrix element squared into each component of a sum
for term in msq_symb.as_coeff_add()[1]:
#and for each term we look at the numerator and denominator
num, denom = sympy.cancel(term).as_numer_denom()
#if each denominator is in the form of a constant times
# one of the accepted forms (1,s,t,u,s^2,t^2,u^2) then
#we are good. Otherwise we need to do extra work
denomtest=True
for den in wanteddenons:
ratioden = denom/den
#we want this ratio to be independent of the mandelstam variables
if (sympy.simplify(sympy.Derivative(ratioden,s))==0 and\
sympy.simplify(sympy.Derivative(ratioden,t))==0 and\
sympy.simplify(sympy.Derivative(ratioden,u))==0):
denomtest=False
if denomtest:
fixed_denominators=False
#now decompose each numerator in a dictionary of powers, e.g. 2*s*t - > 2:1, s:1, t:1
ndict=num.as_powers_dict()
for var in mandvars:
if var in ndict:
#if a variable is at the numerator, exclude it from the denominator
newdenom=denom.subs(var,var-s-t-u)
#rewrite the numerator by lowering one power of that variable
newnum=num/var*sympy.simplify(var-s-t-u)
#add this new term to the new expression
newexpr=newexpr+newnum/newdenom
#do this only once
break
#add an escape route after a certain number of iterations
niters = niters+1
if not fixed_denominators:
msq_symb=sympy.simplify(newexpr).expand()
if niters>= nitersmax:
raise RuntimeError(f"failed to reduce the matrix element squared {msq} "
f"to a compatible form after {niters} iterations. Obtained {msq_symb}")
#the way to proceed is to decompose the matrix element into its basis elements
#with this we separate the matrix element squared into each component of a sum
if msq_symb.as_coeff_add()[0] != 0:
raise RuntimeError(f"The matrix element squared {msq_symb} "
"does not appear to be in the form of a sum")
newmsq=sympy.S("0")
#define the overall sign factor that accounts for n_F=-n_{-1}
#BETA: account for double production by adding the possibility that stats[0]==0. In
#case of double production the two incoming states [1] and [2] have the same statistics
if stats[0] == 0:
overall =1
else:
#if not dealing with pair production, the sign is given by tau_1, i.e. stats[0]
overall = stats[0]
#create a dict for the result
result={}
for term in msq_symb.as_coeff_add()[1]:
finalterm = term
#find the u dependence of each term
ucoeff = term.as_coeff_exponent(u)
#ucoeff is ( a,b) for each a u^b
#if u is at the denominator: do a u<->t switch
if ucoeff[1]<0:
finalterm = term.subs({t: -s-t, u:t },simultaneous=True)
if ucoeff[1] ==-1 and stats[1]!= stats[2]:
#create a new tuple with exchanged s1 s2 statistics
statterm = stats[0], stats[2], stats[1]
#destroy the variable, so that it does not end up in the sum later
strterm=str(finalterm)
finalterm=sympy.S("0")
#recursively call the function on the new matrix element
newresult = decompose(strterm,statterm,couplings)
#if such a new statistics has already been called, we must add this new term to it
if statterm in result:
tempresult = result[statterm]
for i, tempres in enumerate(tempresult):
tempresult[i]=sympy.sympify(tempres)+sympy.sympify(newresult[statterm][i])
result[statterm] =tempresult
else:
#otherwise we just add it to the dict of the results
result.update(newresult)
else:
#if the denominator is more IR sensitive than 1/u we throw an exception
if ucoeff[1]<-1:
raise RuntimeError(f"This code cannot handle a term like {term} "
f"for a dimension-{dim} coupling")
#if u is at the numerator: do a u = -s -t substitution
if ucoeff[1]>0:
finalterm = term.subs(u,-s-t)
newmsq = newmsq+finalterm*overall
#create lists for the terms that go like s/t, 1 (in t-channel parametrisation)
# and u/s and (s^2+u^2)/t, t, t^2/s and s in dim 5
a1=sympy.S("0")
a2=sympy.S("0")
b1=sympy.S("0")
if dim==5:
b2=sympy.S("0")
#again throw an exception if the result does not appear in the form of a sum
if newmsq.simplify().expand().as_coeff_add()[0] != 0:
raise RuntimeError(f"The matrix element squared {newmsq} "
"does not appear to be in the form of a sum")
for term in newmsq.simplify().expand().as_coeff_add()[1]:
#assume that the result term by term is a power-law in t,
# and look separately at each possible power
#we later check that no piece of the original matrix element squared is left behind
tsplit = term.as_coeff_exponent(t)
if dim ==4:
if tsplit[1]==-1:
#then this is a s/t
a1 =a1 +tsplit[0]/s
if tsplit[1]==1:
#then this is t/s
b1= b1 -tsplit[0]*s
a2 = a2 -tsplit[0]*s
if tsplit[1]==0:
#then this is a constant
a2 = a2 +tsplit[0]
if dim ==5:
if tsplit[1]==-1:
#then this is a s^2/t
a1 =a1 +tsplit[0]/s/s/2
a2 = a2- tsplit[0]/s/s/2
b2 = b2 - tsplit[0]/s/s
if tsplit[1]==1:
#then this is t
a2 = a2 +tsplit[0]
if tsplit[1]==2:
#then this is t^2/s
b1 = b1 +tsplit[0]*s
if tsplit[1]==0:
#then this is s
b2 = b2 +tsplit[0]/s
#if the dimension is not 4 or 5 and we are falling back to 4D numerics,
# then we just need to separate t at the numerator (s-parametrisation)
# from t at the denominator (t-parametrisation). These will be divergent, so
#their evaluation will not be possible at the moment
#we use a1 to store t denominators and b1 to store the rest
if fallback:
if tsplit[1]<0:
a1 = a1 + term
else:
b1 = b1 + term
#test the new expression
checkmsq=sympy.sympify('0')
if dim ==4:
checkmsq = a1*s/t+b1*(-t/s-1)+a2
if dim ==5:
checkmsq = a1*(2*s*s/t+2*s+t)+a2*t+b1*t*t/s+b2*s
if fallback:
checkmsq= a1 + b1
#check that the new expression agrees with the old one,
# i.e. that the basis decomposition went well
if sympy.simplify(checkmsq-newmsq) != 0:
raise RuntimeError(f"The dimension-{dim} decomposition "
f"for expression {sympy.simplify(newmsq)} has failed")
#return the dict of the coefficients
if dim==4:
result.update({stats: [a1,a2,b1] })
if dim==5:
#exploit \int s+2t=0 when s_1=s_2
if stats[1]==stats[2]:
result.update({stats: [a1,a2-2*b2,b1,sympy.S("0")] })
else:
result.update({stats: [a1,a2,b1,b2] })
if fallback:
result.update({stats: [a1,b1]})
return result
[docs]
class NumRate:
r"""
This class maps the 2<->2 phase-space integration for a dict
``statmsqdict`` of matrix elements squared and statistics to
all the coefficients :math:`a_i` and :math:`b_i`. Its method
:py:func:`~NumRate.rate` returns a function performing the
two-dimensional numerical integration.
In cases where IR divergences are present, HTL resummation is
performed in the strict LO, subtracted and tuned schemes.
In cases of matrix elements squared arising from higher-dimension
operators, a C fallback routine is called for
their 4D integrations if no t-channel divergence is present
:param dict statsmsqdict: uses tuples of (t1,s1,s2) statistics as keys
and strings for the corresponding matrix elements squared. t1,s1 and s2 are +1
for bosons and -1 for fermions. They correspond to the statistics of the
outgoing and two incoming bath particles.
:param dict couplingdict: is a dict with three keys. ``noneq`` is the
non-equilibrium coupling, ``gauge`` are the gauge couplings (a
tuple) and ``others`` all other couplings
:param tuple masses: is needed for IR-divergent processes. In case of
a dimension-4 operator a string with the analytical expression of
:math:`m_f^2/T^2` is expected, with :math:`m_f^2` the asymptotic mass
of the fermion coupling to the non-equilibrium particle. In case of
dimension 5 it is a tuple of strings containing :math:`m_D^2/T^2`
for each gauge boson, with :math:`m_D` the Debye mass
:param int deg: is the degeneracy of the produced particle. The returned
rate is then the production rate per degree of freedom. If this is
not needed, ``deg`` can be set to 1
Examples
========
Axion production in the KSVZ model
.. code-block:: python
from numerical.manipulate import *
axiondict = {(1,1,1):"32*kappa^2*8*g3^6*3*(s*u/t+t*u/s+s*t/u)",\
(-1,-1,1):"32*kappa^2*(s^2+u^2)/t*(-1/2*Nf*8*g3^6)",\
(-1,1,-1):"32*kappa^2*(s^2+t^2)/u*(-1/2*Nf*8*g3^6)",\
(1,-1,-1):"32*kappa^2*(u^2+t^2)/s*(1/2*Nf*8*g3^6)"}
axionclass=NumRate(axiondict, {"noneq":("kappa",),\
"gauge": ("g3",), "others": ("Nf",)},("g3^2*(1+Nf/6)",),1)
axionclass.rate(numpy.array([1.,3.,5.]),1,(numpy.sqrt(4.*numpy.pi*0.3),3),0)
"""
def __init__(self,statmsqdict: dict, couplingdict: dict,\
masses: Union[str,Tuple[str]],deg: int):
self.deg=deg
"""Stores the degeneracy of the produced particle"""
#this constructor takes as input the R.H.S of a Boltzmann eq.
# in the form of eq:prodtot in the autotherm note summed over all matrix
# elements and statistics and iterates over each combination of
# statistics x matrix elements to yield a function to compute
# the rate decomposed in all the basis elements, as a function of
# the momentum k/T and of the numerical value of all coupling
#The factor of 1/(4 deg k) is pulled out, where deg is the
# degeneracy of the particle being produced, e.g. 1 for axion, 2 for photon
#it takes as input a dict of {stats : "msq"} and another dict of the couplings
#msq is to be intended as summed over the degeneracies of all species
#dim is the dimensionality of the coupling to the bath
# couplingdict has to have the following structure, with keys 'noneq',
# 'gauge' (needed for gauge-boson masses), 'others' (which may also
# include stuff like Nf)
#for each key we have a tuple of string variable
#masses is the expression for the mass. if dim 4 it is the fermion's mass squared over T^2.
#if dim 5 it is a tuple of the masses for each gauge boson, in the same order
# or the couplings in the 'gauge' dict entry
#get the symbolic names of the couplings
noneq = couplingdict["noneq"]
gauge = couplingdict["gauge"]
others = couplingdict["others"]
couplings=noneq+gauge+others
#first of all add mandelstam variables and couplings to the list of known symbols
nslist={}
for coupling in couplings:
nslist[coupling] = sympy.Symbol(coupling)
#save the non-equilibrium coupling
self.ynoneq=nslist[noneq[0]]
"""Stores the :external+sympy:py:class:`sympy symbol <sympy.core.symbol.Symbol>`
for the non-equilibrium coupling"""
s = sympy.Symbol("s")
t = sympy.Symbol("t")
for mandelstam in ["s","t","u"]:
nslist[mandelstam] = sympy.Symbol(mandelstam)
#create a list of symbols for the bath and gauge couplings
self.bathcouplings=[]
"""Stores the :external+sympy:py:class:`sympy symbol <sympy.core.symbol.Symbol>`
for the equilibrium couplings"""
self.gaugecouplings=[]
"""Stores the :external+sympy:py:class:`sympy symbol <sympy.core.symbol.Symbol>`
for the gauge couplings"""
for coupling in gauge+others:
if coupling in gauge:
self.gaugecouplings.append(sympy.Symbol(coupling))
self.bathcouplings.append(sympy.Symbol(coupling))
#TODO: fix when there are no bathcouplings
#create an empty dict for the results
self.results={}
"""Stores a ``stats``, ``msq`` dict for the decomposition into the basis functions.
This follows the definition given in :py:func:`decompose`"""
#run the decomposition on each stats*msq pair
for stats, msq in statmsqdict.items():
#run the decomposition
temp = decompose(msq,stats,couplings)
#if there is already an entry in the dictionary for this statistics, update it
for resstat in list(temp):
if resstat in self.results:
oldresults = self.results[resstat]
for i,oldres in enumerate(oldresults):
oldresults[i]=sympy.simplify(sympy.sympify(oldres)\
+sympy.sympify(temp[resstat][i]))
self.results[resstat] = oldresults
#if not, just add it
else:
self.results[resstat]=temp[resstat]
#reconstruct the dimension from the output of the decomposition
self.dim = len(list(self.results.values())[0])+1
"""Stores the dimension of the operator coupling the produced particle to the bath"""
#prepare to use the fall-back method if D=!4 and D!=5
self.fallback=False
"""If ``True`` the C fallback mode is active"""
if self.dim != 4 and self.dim != 5:
self.fallback=True
#prepare to ignore dict entries having vanishing coefficients
emptyset = [0,0,0]
if self.dim ==5:
emptyset.append(0)
#switch to find out if we need HTL resummation
self.irdiv = False
"""If ``True`` the naive rate is IR divergent"""
#dict in which the divergent pieces are stored
ircoeffs={}
#dict in which the clean results are stored in the form of lambdified
# functions of the bath couplings
self.cleanresults={}
#create a set to store the power of the non-equilibrium coupling
powerlaws=set()
for stats, msq in self.results.items():
#start by filtering out entries having vanishing coefficients
if msq != emptyset:
#create a vector to store the cleaned-up matrix elements
msqsimple=emptyset.copy()
#deal with the 0th coefficient (a1) separately, as it is the
#one which is divergent
msq[0]= msq[0].simplify()
#If it has a non-zero component for any statistics, we have an IR divergence
#and we need resummation
if msq[0] !=0:
self.irdiv = True
#if the dimensionality is unsupported, we can't treat this case at the moment.
# For instance, Gs->Gs (graviton-scalar scattering) would go like
# \kappa^4 s^2u^2/t^2 and we would need the graviton HTL to
# address its IR behavior
if self.fallback:
warnings.warn(f"unsupported divergent expression {msq[0]} for "
f"statistics {stats}. Setting it to zero",AutothermWarning,stacklevel=2)
self.irdiv = False
#in dimension four, simply store this coefficient
if self.dim==4:
ircoeffs.update({stats:msq[0]})
#each msq[i] will be of the form a noneq^b. extract a and b
coeffexp=msq[0].as_coeff_exponent(self.ynoneq)
#create a lamdbdified function for this coefficient
if self.dim ==4:
msqsimple[0] = sympy.lambdify(self.bathcouplings,\
coeffexp[0],modules="numpy")
if self.dim ==5:
#in this case we need to go gauge coupling by gauge
# coupling and create a vector function
#start by converting this expression to a poly object
polya1 = coeffexp[0].as_poly(self.gaugecouplings)
#create a dict of the coefficient of the different powers
# in the gauge couplings
polya1dict= polya1.as_dict()
#create an array to store the output split by gauge group
a1gauge=len(self.gaugecouplings)*[0]
#go through them and create the coefficient for each coupling
# such as g_i^2 g_j^2 (s^2+u^2)/t
for powers, a1coeff in polya1dict.items():
for i, gaugec in enumerate(self.gaugecouplings):
if powers[i] !=0:
# a1gauge[i] ==0
# else:
a1gauge[i]=a1coeff*gaugec**powers[i]
#now create a lambdified vector function (will return a numpy array)
msqsimple[0] = sympy.lambdify(self.bathcouplings,a1gauge,modules="numpy")
#and store this function in the dict
ircoeffs.update({stats:a1gauge})
#add the power law exponent of the non-eq coupling to the list
powerlaws.add(coeffexp[1])
else:
#we still need a function to give 0 if there is no divergence
if self.dim==4:
msqsimple[0] = sympy.lambdify(self.bathcouplings,0,modules="numpy")
if self.dim==5:
msqsimple[0] = sympy.lambdify(self.bathcouplings,\
len(self.gaugecouplings)*[0],modules="numpy")
#now go through the coefficients of the non-divergent elements
for i in range(1,len(msq)):
msq[i]= msq[i].simplify()
#each msq[i] will be of the for a noneq^b. extract a and b
coeffexp=msq[i].as_coeff_exponent(self.ynoneq)
if self.fallback and i==1:
#in this case we need to create the type of structure
# handled by the fallback C integration assuming the
# non-divergent contribution is of the form
# \sum_i c_i s^{m_i} t^{l_i}, with m_i+l_i=dim -4
# and l_i>=0, we need the (c_i, m_i, l_i) elements
#TODO: add a check for this!
#the polynomial method fails if there are negative powers,
# so we keep multiplying by s until it works
tempexpr=coeffexp[0].simplify().expand()
spower=0
while tempexpr.as_poly(s,t) is None:
tempexpr = s*tempexpr
spower = spower+1
#now we can create a dict
polydict = tempexpr.as_poly(s,t).as_dict(s,t)
# initialise an empty array to store the results
fallbackvec=[]
for stpowers, coeffs in polydict.items():
#create a lambidified function for the c_i
fallbackvec.append(sympy.lambdify(self.bathcouplings,\
coeffs.as_expr(),modules="numpy"))
#add the values of m_i and l_i. For the former fix the extra
# powers of s we added before
fallbackvec.append(stpowers[0]-spower)
fallbackvec.append(stpowers[1])
msqsimple[i]= fallbackvec
else:
#create a lambdified function of the bath couplings for each coefficient
msqsimple[i] = sympy.lambdify(self.bathcouplings,\
coeffexp[0],modules="numpy")
#do not consider cases where the matrix element coefficient is zero
if msq[i] !=0:
powerlaws.add(coeffexp[1])
#add this cleaned-up matrix elements to the new dict
self.cleanresults.update({stats:msqsimple})
#check that we have no more than a single power law
if len(powerlaws) >1:
raise NotImplementedError("Non-homogeneous matrix element in terms of the "
f"non-equilibrium coupling {noneq[0]}\n"
"This is currently unsupported, please split your matrix "
"element squared into homogeneous subsets.")
else:
#save the exponent in a variable
self.powerlawexp=powerlaws.pop()
#if we have an IR divergence, we need to make sure the thermal masses
# match with what comes out of the IR limit of the matrix elements.
if self.irdiv:
#start with the 4-dimensional case
if self.dim==4:
#here I can only compute the IR limit of the integrations,
#first a sanity check: the coeff of the (-1,1,whatever) piece must be opposite to
#that of the (1,-1,whatever) coefficient. We also check that there is a single
#value of "whatever" contributing and fail otherwise (it would mean that the
# non-equilbrium particle is both a boson and a fermion)
#start by introducing a control variable
foundcoeff = False
poscoeff=sympy.sympify('0')
for i in (-1,1):
if (-1,1,i) in ircoeffs:
coeff = ircoeffs[(-1,1,i)]
if coeff !=0:
#store for later use the value of s2
self.s2div=i
if sympy.simplify(coeff + ircoeffs[(1,-1,i)]) != 0 :
raise RuntimeError("Lack of crossing symmetry for dimension-4 IR "
"coefficients "+str(coeff)+" and "\
+str(-ircoeffs[(1,-1,i)]))
else:
poscoeff=coeff
if foundcoeff:
raise RuntimeError("Multiple nonzero dimension-4 IR "
"coefficients "+str(coeff)+" "\
+str(-ircoeffs[(1,-1,i)]))
foundcoeff = True
#let us start by stripping the dependence on the noneq
# coupling constant and multiply by pi^2/4,
#which is the typical T^2 \pi^2/4 accompanying a fermion mass
#ratsimp does a rational simplification
poscoeffstrip=sympy.ratsimp(poscoeff.subs(noneq[0],1)*sympy.pi*sympy.pi/4)
self.mass=(sympy.sympify(masses)).ratsimp()
#check that the ratio between the mass and the IR-divergence coefficient is real
#meaning that it is a pure number and does not depend on any constant, which
#would mean that there is a mismatch between the two.
#NOTE: This could be problematic for some models, say if an overall constant
#considered as a "coupling" multiplies the divergent matrix element
#for this reason we print out a warning but do not fail
#TODO: consider using the is_constant() method
self.massratio_expr=sympy.ratsimp(poscoeffstrip/self.mass)
if not self.massratio_expr.is_real :
warnings.warn("The ratio between the fermion mass "\
+str(self.mass)+" and the coefficient of the IR "
"divergence "+str(poscoeffstrip)+" depends on"
" the coupling constants entering the mass.\n"
"This is most likely not problematic, continuing with "
"evaluation.",AutothermWarning,stacklevel=2)
#lambdify this coefficient
self.massratio=sympy.lambdify(self.bathcouplings,self.massratio_expr)
#the value of \xi is
self.xi=0.5*numpy.exp(1.)
#create a lambdified function for the mass
self.massfun=sympy.lambdify(self.bathcouplings,self.mass)
#continue with the gauge-boson exchange in the dim-5 case
if self.dim==5:
#here I can only compute the IR limit of the integrations, along the lines
#of what is done starting around eq. eq:irlmitstart in the autotherm note
#first we collect the bosonic and fermionic parts of the divergence, i.e. the
#coefficients of the (1,1,whatever) and (-1,-1,whatever) statistics.
#note that whatever has to correspond to the statistics of the tracked
# particle of momentum k, as we are looking at a t-channel gauge-boson exchange
#these are tracked by the ircoeffs vector
#we now use this to add a safety check on the masses
massescoeffs=len(self.gaugecouplings)*[0]
for i in range(0,len(self.gaugecouplings)):
for stats, a1coeffs in ircoeffs.items():
if stats[0]==1 and stats[1]==1:
#we have a boson. The corresponding factor is \pi^2/6
massescoeffs[i]=massescoeffs[i]-a1coeffs[i]/6*sympy.pi*sympy.pi
if stats[0]==-1 and stats[1]==-1:
#we have a fermion. The corresponding factor is -\pi^2/12,
# with a factor of -1/2 from statistics and \tau_1 normalisation
# FIXME refer to equations in the manual
massescoeffs[i]=massescoeffs[i]+a1coeffs[i]/12*sympy.pi*sympy.pi
#if the two statistics are different, something is very wrong
if stats[0]!= stats[1]:
raise RuntimeError("Different statistics for the two particles "
f"connected to a t-channel gauge boson {a1coeffs}")
massescoeffs[i]=sympy.ratsimp(massescoeffs[i])
#compare the IR coefficients with the masses we obtained
self.massvec=[]
if isinstance(masses,tuple):
for item in masses:
masstmp=(sympy.sympify(item)).ratsimp()
self.massvec.append(masstmp)
#in this case we have a single mass
else:
masstmp=(sympy.sympify(masses)).ratsimp()
self.massvec.append(masstmp)
self.massratio=[]
self.massratio_expr=[]
for i in range(0,len(self.gaugecouplings)):
#compute the ratio between the two and store it in a function,
# it will be needed later to compute the coefficient of the HTL log
massratiotemp=sympy.ratsimp(massescoeffs[i]/self.massvec[i])
self.massratio_expr.append(massratiotemp)
self.massratio.append(sympy.lambdify(self.bathcouplings,\
massratiotemp,modules="numpy"))
#do again the check that the ratio is a pure number. This fails in the axion
# case if the explicit powers of the gauge coupling in the
# effective axion-gauge coupling are not included in non-eq
if not massratiotemp.is_real :
warnings.warn("The ratio between the gauge boson "
"mass "+str(self.massvec[i])+" and "
"the coefficient of the IR divergence "+str(massescoeffs[i])+" depends on"
" the coupling constants entering the mass.\n"
"This is most likely not problematic, continuing with evaluation.",\
AutothermWarning,stacklevel=2)
#the value of \xi is
self.xi=0.5*numpy.exp(1./3.)
#create a lambdified function for the mass
self.massfun=sympy.lambdify(self.bathcouplings,self.massvec)
[docs]
def get_coeffs(self):
"""This method returns the :external+sympy:py:class:`symbolic form <sympy.core.expr.Expr>`
of the coefficients"""
return self.results
[docs]
def get_leadlog(self):
"""This method returns the :external+sympy:py:class:`symbolic form <sympy.core.expr.Expr>`
of the leading-logarithmic part"""
#in this case there is no log
if not self.irdiv :
return sympy.sympify('0')
else:
k=sympy.S('k')
temp=sympy.S('T')
prefac = self.ynoneq**self.powerlawexp/(128*sympy.pi**3*self.deg*k)
if self.dim ==4:
msquared = self.mass*temp**2
logcoeff = self.massratio_expr*msquared*(sympy.exp(k/temp)+self.s2div)
logcoeff = logcoeff*sympy.log(4*k*k/msquared)/(sympy.exp(k/temp)-self.s2div)
if self.dim==5:
msquared= self.massvec
logcoeff = 0
#in this case we need to sum over the gauge groups
for i, msqi in enumerate(msquared):
if msqi != 0:
logcoeff += 8*k*self.massratio_expr[i]*msqi*temp**3*\
sympy.log(4*k*k/(msqi*temp**2))
return sympy.simplify(prefac*logcoeff)
[docs]
def rate(self,k: Union[float,numpy.ndarray], y: float, numcouplings: Tuple[float],
mode: int, intoptions : Optional[list]=None) ->numpy.ndarray:
r"""
This method returns the rate :math:`\Gamma/T^{1+2dim-8}` as a function of an
array ``k`` of :math:`k/T` values, ``y`` value for the non-equilibrium coupling, a
``numcoupling`` tuple with the numerical values of the other couplings.
:param ndarray k: an array ``k`` of :math:`k/T` values
:param float y: the non-equilibrium coupling
:param tuple numcoupling: the numerical values of the other couplings
:param int mode: sets the desired output. 0 returns everything, 1 returns
the strict rate only, 2 the subtr rate only, 3 the tuned rate only.
A negative value uses the C fallback routines
:param list intoptions: can be passed optionally as a dict or lists of
dicts for the two nested integrations in :math:`q_+` and :math:`q_-`.
See the :external+scipy:py:func:`scipy documentation<scipy.integrate.quad>`.
The default options are ``{'limit':1000,'epsrel':1e-5}`` for both integrations.
:return: a list of ndarrays() with `k` in the first position and the
requested rates according to `mode`
:rtype: ndarray
"""
#first of all, fail if k is not a numpy array or a single float or integer
if not isinstance(k,numpy.ndarray):
if isinstance(k,float) or isinstance(k,int) or isinstance(k,numpy.float64):
k=numpy.array([k])
else:
raise TypeError("unexpected type of variable k")
klen=len(k)
#mode works as follows
# <0 - use fallback C routines
# 0 - returns k/T, strict, subtr, tuned
# 1 - returns k/T, strict
# 2 - returns k/T, subtr,
# 3 - returns k/T, tuned
# for positive modes, ignore this distinction if the expression is not
# infrared divergent, as all methods are bound to agree
cfallback=False
if mode <0:
cfallback = True
if not self.irdiv :
mode = 1
mass = 0.
else:
#compute m/T
msquared= self.massfun(*numcouplings)
mass = numpy.sqrt(msquared)
if mode >3 or not isinstance(mode, int):
raise ValueError(f"Invalid value of mode: {mode}")
if mode ==1 or mode ==2:
tuned = False
else:
tuned = True
if mode !=3:
subtract = True
else:
subtract = False
#compute the overall prefactor
prefac= y**self.powerlawexp/(4*self.deg*k)
if not self.fallback and not cfallback :
prefac=prefac/((4.*numpy.pi)**3*k)
#create the variables for the t and s parts of the result in the subtr
# and strict schemes, as well as in the C fallback
if subtract:
tpartsubtr=numpy.zeros(klen)
if tuned:
tparttuned=numpy.zeros(klen)
spart = numpy.zeros(klen)
#in the dim-5 case it is convenient to split out the finite
# part (a2-proportional) of the t-channel
tpartfinite=numpy.zeros(klen)
# Save the original warning handler
original_showwarning = warnings.showwarning
#prepare a custom function for the integration warnings thrown by scipy.integrate
def create_custom_warning_handler(k,stat,intmessage,aibi):
def custom_warning_handler(message, category, filename, lineno, file=None, line=None):
if category is IntegrationWarning:
print(f"Warning raised by the integration routine for the {intmessage} for "
f"momentum {k}, statistics {stat} and coefficients {aibi}.\nThis is "
f"possibly not problematic. The warning message is:\n\n{message}")
else:
original_showwarning(message, category, filename, lineno, file, line)
return custom_warning_handler
# Restore the original warning handler
warnings.showwarning = original_showwarning
#set up the default integration options
if intoptions is None:
intoptions=[{'limit':1000,'epsrel':1e-5},{'limit':1000,'epsrel':1e-5}]
#now start to cycle through the statistics
if self.fallback:
for stats, msqfun in self.cleanresults.items():
coeffslist=msqfun[1]
n=int(len(coeffslist)/3)
#create a new list with evaluated lambdified functions
coeffslistnew=coeffslist.copy()
for i in range(0,n):
coeffslistnew[3*i]=coeffslist[3*i](*numcouplings)
for j, kj in enumerate(k):
spart[j] = spart[j]+integrator.fallback_s_poly_vec(kj,*stats, n,\
numpy.asarray(coeffslistnew,dtype=float))
else:
for stats, msqfun in self.cleanresults.items():
#evaluate the lambdified functions of the couplings
a1 = msqfun[0](*numcouplings)
a2 = msqfun[1](*numcouplings)
b1 = msqfun[2](*numcouplings)
if self.dim==5:
b2 =msqfun[3](*numcouplings)
#now prepare for the integrations. First we handle the case where we
# are using the C fallback routines
if cfallback:
#work through the different cases
if self.dim ==4:
if a1 !=0 or a2!=0:
paramvec = numpy.asarray([a1,a2,self.xi*mass])
for j, kj in enumerate(k):
tpartfinite[j] += integrator.fallback_t_poly_vec(kj,\
*stats,self.dim,paramvec)
if b1 !=0:
paramvec = numpy.array([-b1,0.,0.,-b1,-1.,1.])
for j, kj in enumerate(k):
spart[j] += integrator.fallback_s_poly_vec(kj,*stats,2,paramvec)
if self.dim ==5:
#store separately the finite part
if a2!=0:
paramvec=numpy.array([0.,a2,0.])
for j, kj in enumerate(k):
tpartfinite[j] += integrator.fallback_t_poly_vec(kj,*stats,\
self.dim,paramvec)
#now go gauge group by gauge group
for i, a1i in enumerate(a1):
if a1i!=0. or a1i!=0:
paramvec=numpy.array([a1i,0.,self.xi*mass[i]])
for j, kj in enumerate(k):
tpartfinite[j] += integrator.fallback_t_poly_vec(kj,\
*stats,self.dim,paramvec)
#now the s-channel
if b1 !=0 or b2!=0:
paramvec= numpy.array([b1,-1.,2.,b2,1.,0.])
for j, kj in enumerate(k):
spart[j] += integrator.fallback_s_poly_vec(kj,*stats,2,paramvec)
else:
if self.dim==4:
if a1 !=0 or a2!=0:
#deal with the different cases
if subtract:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,stats,\
"dim-4 t-channel subtracted",[a1,a2])
tpartsubtr[j] += nquad(t4fast,[[0.,kj],[-numpy.inf,0.]],\
args=(kj,a1,a2,*stats,-1),opts=intoptions)[0]
if tuned:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,stats,\
"dim-4 t-channel tuned",[a1,a2])
#recall that t4fast and t5fast take care automatically of the
# a_1,a_2 manipulations needed to deal with the tuned scheme
tparttuned[j] += nquad(t4fast,[[0.,kj],[-numpy.inf,0.]],\
args=(kj,a1,a2,*stats,self.xi*mass),opts=intoptions)[0]
if b1 !=0:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,stats,\
"dim-4 s channel",1)
spart[j] += b1*nquad(s4fast,[[kj,numpy.inf],[0.,kj]],\
args=(kj,*stats),opts=intoptions)[0]
if self.dim==5:
#store separately the finite part
if a2!=0:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,\
stats,"dim-5 t-channel a_2",1)
tpartfinite[j] += a2*nquad(t5fast,[[0.,kj],[-numpy.inf,0.]],\
args=(kj,0.,1.,*stats,-1),opts=intoptions)[0]
if subtract:
#in this case we do not need to go gauge group by gauge group,
# we can sum the a1 vector
a1tot=numpy.sum(a1)
if a1tot !=0:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,\
stats,"dim-5 t-channel subtracted",1)
tpartsubtr[j] += a1tot*nquad(t5fast,[[0.,kj],[-numpy.inf,0.]]\
,args=(kj,1.,0.,*stats,-1),opts=intoptions)[0]
if tuned:
for i,a1i in enumerate(a1):
if a1i!=0. or a1i!=0:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,\
stats,"dim-5 t-channel tuned",1)
tparttuned[j] += a1i*nquad(t5fast,[[0.,kj],[-numpy.inf,0.]]\
,args=(kj,1.,0.,*stats,self.xi*mass[i]),opts=intoptions)[0]
if b1 !=0 or b2!=0:
for j, kj in enumerate(k):
warnings.showwarning = create_custom_warning_handler(kj,stats,\
"dim-5 s channel",[b1,b2])
spart[j] += nquad(s5fast,[[kj,numpy.inf],[0.,kj]],\
args=(kj,b1,b2,*stats),opts=intoptions)[0]
#add tpartfinite to spart to simplify the expressions below
spart = spart +tpartfinite
#in the subtract case, find the logarithmic parts
if subtract:
if not self.irdiv :
logcoeffstrict=numpy.zeros(klen)
logcoeffsubtr=numpy.zeros(klen)
else:
if self.dim ==4:
logcoeffstrict=self.massratio(*numcouplings)*msquared*2*k*\
(1+self.s2div*numpy.exp(-k))/(1-self.s2div*numpy.exp(-k))
logcoeffsubtr=logcoeffstrict*numpy.log1p(4*k*k/msquared)
logcoeffstrict=logcoeffstrict*numpy.log(4*k*k/msquared)
if self.dim==5:
logcoeffstrict=numpy.zeros(klen)
logcoeffsubtr=numpy.zeros(klen)
#in this case we need to sum over the gauge groups
for i, msqi in enumerate(msquared):
if msqi != 0:
logcoeffstrict += 16*k*k*((self.massratio[i])(*numcouplings))\
*msqi*numpy.log(4*k*k/msqi)
logcoeffsubtr += 16*k*k*((self.massratio[i])(*numcouplings))\
*msqi*numpy.log1p(4*k*k/msqi)
#prepare to return the output
#NOTE: the dimensionful factor of T is scaled out
if mode ==0:
return [k,prefac*(spart + tpartsubtr+logcoeffstrict),\
prefac*(spart + tpartsubtr+logcoeffsubtr),\
prefac*(spart + tparttuned)]
if mode ==1:
return [k,prefac*(spart + tpartsubtr+logcoeffstrict)]
if mode ==2:
return [k,prefac*(spart + tpartsubtr+logcoeffsubtr)]
if mode ==3:
return [k,prefac*(spart + tparttuned)]
if mode < 0:
return [k,prefac*spart]
