Defining and configuring a model#

We consider for illustration two models. The first, henceforth toy model, features two real scalars \(\phi\) and \(\sigma\), with Lagrangian

\begin{align} \mathcal{L}_\mathrm{toy}=&\frac12\partial_\mu\phi\partial^\mu\phi+\frac12\partial_\mu\sigma\partial^\mu\sigma -\frac{\lambda}{4!}\phi^4-\frac{g}{3!}\phi^3\sigma\,, \end{align}

where \(\lambda\) is a quartic coupling that can keep the \(\phi\) scalars in equilibrium and \(g\) is the smaller coupling that governs \(\sigma\) production and equilibration.

Our second model is thermal axion (or ALP) production with the most generic axion-SM effective Lagrangian. It reads — see for instance this paper

\begin{align} \mathcal{L}_\mathrm{ALP}=&\,\mathcal{L}_{\mathrm{SM}}+\frac{1}{2}(\partial_\mu a)^2 -\frac{a}{f_\mathrm{PQ}}\left[c_3\frac{\alpha_3}{8\pi}G^b_{\mu\nu}\tilde{G}^{\mu\nu}_b+c_2\frac{\alpha_2}{8\pi}W^j_{\mu\nu}\tilde{W}^{\mu\nu}_j+c_1\frac{\alpha_1}{8\pi}B_{\mu\nu}\tilde{B}^{\mu\nu}\right] \\ &+ic_t y_t\frac{a}{f_\mathrm{PQ}}\bar{Q}_L\sigma_2\phi^{*}t_R+\text{h.c.}\,, \end{align}

where \(f_\mathrm{PQ}\) is the Peccei-Quinn scale, \(B,\;W,\;G\) are the U(1), SU(2) and SU(3) field-strength tensors, \(Q_L\) is the third-generation quark doublet with \(y_t\) the top Yukawa coupling. \(\phi\) is the Higgs doublet. The \(c_i\) are arbitrary Wilson coefficients.

The model file#

We recommend creating a directory for each model. In the following, we shall call $ModelDir the absolute path to this directory, e.g. /Users/myself/AUTOTHERM/models/axion/ on a mac or /home/myself/AUTOTHERM/models/axion/ on Linux.

To translate the Lagrangian in (2) into a FeynRules model file, we use its support for multiple model files. An implementation of the symmetric-phase SM is provided by AutoTherm, so we only need to add the axion field and its couplings.

We thus create a toymodel.fr/ axionSM.fr in the $ModelDir directory. A first step is the definition of model metadata, i.e.

M$ModelName = "Two-scalar model";

M$Information = {
    Authors -> {"AutoTherm User"},
    Version -> "1.0"
};

M$ModelName = "Axion model";

M$Information = {
    Authors -> {"AutoTherm User"},
    References -> {"1310.6982","2404.06113"},
    Version -> "1.0"
};

In the toy model, both fields have to be introduced in the M$ClassesDescription variable. In the axion case, as we shall import the AutoTherm SM model file, only the axion field needs adding. Explicitly

M$ClassesDescription = {
    S[1] == {
        ClassName -> phi,
        SelfConjugate -> True,
        Mass -> 0,
        Width -> 0,
        PropagatorLabel -> "\\phi",
        PropagatorType -> ScalarDash,
        PropagatorArrow -> None
    },
    S[2] == {
        ClassName -> sigma,
        SelfConjugate -> True,
        Mass -> 0,
        Width -> 0,
        PropagatorLabel -> "\\sigma",
        PropagatorType -> ScalarDashes,
        PropagatorArrow -> None
    }
};

M$ClassesDescription = {
    S[2] == {
        ClassName -> ax,
        SelfConjugate -> True,
        Mass -> 0,
        Width -> 0,
        PropagatorLabel -> "a",
        PropagatorType -> ScalarDash,
        PropagatorArrow -> None
    }
};

See the FeynRules documentation for further information on the meaning of each property.

Before we can define the new terms in the Lagrangian, we need to introduce the new parameters in (1)/ (2). This is done through the M$Sparameters variable, i.e.

M$Parameters = {
    lambda == {
        ParameterType -> Internal
    },
    g == {
        ParameterType -> Internal
    }
};

M$Parameters = {
    c1 == {
        ParameterType -> Internal
    },
    c2 == {
        ParameterType -> Internal
    },
    c3 == {
        ParameterType -> Internal
    },
    ct == {
        ParameterType -> Internal
    },
    fPQ == {
        ParameterType -> Internal
    }
};

Finally, we add the necessary interaction terms to the Lagrangian as in

Ltot := -lambda/(4!)phi^4-g/3! phi^3 sigma

LIntAxion := Block[{mu,nu,rho,sigma,aa,ii},
c3 (g3^2/(32 Pi^2 fPQ)) ax FS[G,mu,nu,aa] 1/2 Eps[mu,nu,rho,sigma] FS[G,rho,sigma,aa]
+ c2 (g2^2/(32 Pi^2 fPQ)) ax FS[Wi,mu,nu,ii] 1/2 Eps[mu,nu,rho,sigma] FS[Wi,rho,sigma,ii]
+ c1 (g1^2/(32 Pi^2 fPQ)) ax FS[B,mu,nu] 1/2 Eps[mu,nu,rho,sigma] FS[B,rho,sigma]];

LYukAxion := Block[{sp,sp1,ii,jj,cc,ff1,ff2,yuk,iii,jjj},
    yuk = -ct/fPQ ax yu[ff1,ff2] QLbar[sp,ii,ff1,cc].ProjP[sp,sp1].uR[sp1,ff2,cc] Phibar[jj] Eps[ii,jj];
    yuk + HC[yuk]  ];

Ltot := LSM  + LIntAxion + LYukAxion;

LSM and all the SM fields are defined in analytical/models/symmetric.fr. Kinetic terms for gauge-singlet fields are not needed by FeynRules, as long as the fields are canonically normalized.

The configuration file#

We create a toymodel.cfg/ axionSM.cfg file in the $ModelDir directory. The first section therein is just used to define the absolute paths of the main tools and libraries. math is the math or WolframKernel executable (the command-line version of Mathematica) and feynrules, feynarts and formcalc the directories where these tools have been installed. On a mac we would then have, for either model

[Tools]
math = /Applications/Mathematica.app/Contents/MacOS/WolframKernel
feynrules = /Users/jacopo/NextCloud/AUTOTHERM/feynrules-current
feynarts = /Users/jacopo/NextCloud/AUTOTHERM/FeynArts-3.12/FeynArts.m
formcalc = /Users/jacopo/NextCloud/AUTOTHERM/FormCalc-9.10/FormCalc.m

We then need the model-specific information. This is achieved through

[Model]
modelpath = $ModelDir/toymodel.fr
# Symbol for the Lagrangian in the model file
lagrangian = Ltot
# "Name" of the particle whose production rate must be computed
produced = S[2]
# List of the particles in the thermal bath. Leave empty for all SM particles
inbath = S[1]
# List of additional assumptions
assumptions = lambda > 0, g > 0
# List of replacements
replacements =
# Whether to include the SM model file while generating the Feynman rules or not
includeSM = no
# Non-equilibrium coupling
noneq = g
# fermions that need to be treated generation by generation because of Yukawa-dependent thermal masses
flavorexpand =

[Model]
modelpath = $ModelDir/axionSM.fr
# Symbol for the Lagrangian in the model file
lagrangian = Ltot
# "Name" of the particle whose production rate must be computed
produced = S[2]
# List of the particles in the thermal bath. Leave empty for all SM particles
inbath =
# List of additional assumptions
assumptions = Element[c1,Reals], Element[ht,Reals], Element[c2,Reals],  Element[c3,Reals],  Element[ct,Reals],  Element[fPQ,Reals]
# List of replacements
replacements =
# Whether to include the SM model file while generating the Feynman rules or not
includeSM = yes
# Non-equilibrium coupling
noneq = fPQ
# fermions that need to be treated generation by generation because of Yukawa-dependent thermal masses
flavorexpand =

We refer to our dedicated reference section and to Appendix D of the manual for an in-depth explanation of each flag. For each model discussed in the provided examples we include this section within the provided jupyter notebooks.

Finally, we need to specify what we want the code to do. We set conf, rules and proc to yes, so that it will generate configuration files for its internal workings, Feynman rules with FeynRules and matrix elements squared and thermal masses for axion production.

[Run]
# Whether to generate the <name>.m file for the Mathematica scripts
conf = yes
# Whether to generate the Feynman rules
rules = yes
# Whether to generate all processes
proc = yes
# Verbosity
verbose = no
# Path to the AutoTherm directory
autothermdir = /path/to/autotherm/

autothermdir should be set to $AutothermDir, introduced in the previous section. In the next section we illustrate how to run the code.