Introduction to CosmoLattice

Here we first present an Overview of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\), where we briefly comment about the basic features of the code. In Section Capabilities and structure we further discuss the capabilities of the code and introduce its file structure. In Section Basic Field Equations implemented (so far) we review the basic field equations that \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is ready to solve.

Overview

The application of numerical methods to studying non-linear field dynamics in the early Universe has grown substantially in recent years. Building increasingly efficient, accurate, and robust numerical tools is essential for placing predictions about the early Universe on a firm theoretical foundation. These developments have led to the emergence of a distinct field—Lattice Cosmology—whose strength lies in its ability to resolve detailed non-linear field dynamics and provide reliable predictions for the resulting observables. Lattice Cosmology Techniques (LCT) represent nowadays a well-established route towards a deeper understanding of the early Universe, and we expect them to play an increasingly prominent role in shaping observational strategies designed to probe its physics. It was in this context that \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) was originally developed with a clear aim: to provide a powerful, flexible and robust package, publicly available, and dedicated to solving the non-linear dynamics of field theories in an expanding Universe, thus allowing to investigate the phenomenology and observational implications of early Universe scenarios. At the time of writing (July 2026), the latest version the code, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) \({\tt v2.0}\), can simulate:

Singlet scalar field interactions Abelian U(1) scalar-gauge interactions Non-Abelian SU(2) scalar-gauge interactions Axion like particle - Abelian U(1) gauge field interactions Non-minimally coupled-to-gravity scalar field dynamics Global cosmic defects (domain walls, strings, monopoles, textures) Local topological defects (local strings) Gravitational waves from scalar field theories Gravitational waves from Abelian U(1) scalar-gauge theories

In other words, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) can handle scenarios with: interacting singlet scalar fields; charged scalar fields under a U(1) and/or SU(2) gauge symmetries, with the corresponding Abelian and/or non-Abelian gauge vector fields; axion-like-particles (ALP) interacting with Abelian gauge sectors, non-minimal gravitationally coupled scalar fields, gravitational wave (GW) emission from singlet scalar fields, GW emission from charged scalars and Abelian gauge fields, global and local cosmic defects, and soon enough (see below Upcoming Physics) fluid dynamics and other interactions. \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) can simulate the dynamics of such fields either in a flat space-time, or in a homogeneous and isotropic (spatially flat) expanding background. In the latter case, the fields can evolve either over a fixed background (e.g. with a power-law scale factor), or with self-consistent expansion of the Universe, i.e. with the fields involved in a simulation determining themselves the expansion rate through the Friedmann equations.

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) provides both symplectic integrators, with accuracy ranging from \(\mathcal{O}(\delta t^2)\) up to \(\mathcal{O}(\delta t^{10})\), and non-symplectic integrators with accuracies from \(\mathcal{O}(\delta t^2)\) up to \(\mathcal{O}(\delta t^{4})\) (See Sections 3.3, 3.4 and 3.5 of \(\mathtt{The~Art-I}\) (Ref. ¹), for a discussion on integrators and their properties). Appropriate observables are also provided for each algorithm, like the energy density components of each field, their relevant spectra, or dynamical constraints. Our algorithms conserve energy up to the accuracy set by the order of the evolution algorithm, reaching even machine precision in the case of the highest order integrators. Notably, our algorithms for scalar-gauge theories, either Abelian or non-Abelian, respect the Gauss constraint to machine precision, independently of the integrator, even in the case of self-consistent expansion. This is thank to the fact that we write the discretization of our scalar-gauge theories with lattice-gauge-invariant techniques (to know more about the latter, see e.g. Section 3 of \(\,\texttt{The}\,\texttt{Art}\)-\(\texttt{I}\) (Ref. ¹), or the Lattice gauge invariant techniques tab in this website).

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is written in C++, and fully exploits the object oriented programming nature of this language, with a modular structure that separates well all the ingredients involved. This allows \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) to have a clear separation between the physics and the technical implementation details. The code is designed so that the user can simulate a given scenario with different parameters, without requiring to re-compile each time that parameter values are changed.

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is fully parallelized using both shared memory parallelization (OpenMP, POSIX Threads or GPU acceleration) and distributed parallelization (Message Passing Interface (MPI)) for use in high-performance clusters, and uses a discrete Fourier Transform parallelized in multiple spatial dimensions. This makes it ideal for probing physical problems with well-separated mass/length scales, running very high resolution simulations, or simply shortening the running time of long simulations. To provide these capabilities, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) relies on TempLat, a general package for lattice simulations that defines field variables and their operations, by introducing its own symbolic language, and managing all aspects of dispatching workload to the available computational resources. Once you become familiar with the basic ‘vocabulary’ of the \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) and TempLat language, editing the code or implementing your own model (resembling how you would write it in the continuum), should become a simple task.

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) can be used at multiple levels of complexity:

Basic level

A basic level user with no experience at all in parallelization techniques and with little to no programming experience will be able to run fully parallelized simulations of their favourite models on hundreds of cores in a cluster, while being completely oblivious to the technical details of the algorithm implementation or parallelization. CosmoLattice automatically prints a collection of relevant observables, such as volume averages, field spectra, or dynamical constraints, which can be used to monitor the evolution of the system.

Intermediate level

An intermediate level user with a certain programming experience, will want to modify the type of output CosmoLattice generates, change the initial condition routine (e.g. turn it into a Monte-Carlo generator for thermal configurations), or even design their own evolution algorithms and add these to the family of integrators available in CosmoLattice.

Advanced level

An advanced level user may want to play as well with the hard core inner parts of CosmoLattice and TempLat, in an attempt to understand or to improve the most advanced technical aspects of the code. This could include e.g. the handling of the parallelization, which is something that typically will remain as a black box for the majority of users. This means contributing to TempLat, the lattice framework at the core of CosmoLattice.

What type of user you want to be is entirely up to you: it simply depends on your programming expertise (or lack of it), and on your will to learn how the code works internally.

Upcoming Physics. We have developed \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) with the intention of providing an evolving and up-to-date relevant numerical tool for the scientific community working on the physics of the early universe. While presently (at the time of writing, July 2026) \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is already capable of simulating canonical scalar-gauge field theories, non-minimal scalar theories, Abelian gauge-ALP interactions, global and local cosmic defects, gravitational waves from scalar and/or gauge fields, etc, we have conceived \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) as an evolving package that we plan to upgrade regularly with new features or specialized tasks, or extra modules dedicated to new physics. For example, to mention just a few of the latter, we are currently working on enhancing the capabilities of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) to simulate as well:

Fluid dynamics (fully relativistic and non-relativistic) Fluid - Abelian U(1) gauge field interactions Fluid - scalar field interactions Non-minimal kinetic scalar theories Gravitational waves from Non-Abelian SU(2) scalar-gauge theories Axion like particle - Non-Abelian SU(2) gauge field interactions [...]

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is freely available to anyone who wants to use or modify it, as long as you give us credit for its creation, see How to cite \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\). If you have any questions or comments about \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\), please contact us, as we would love to hear how the program performs for you, and we will be happy to assist you with questions you might have, bug reports, suggestions for future improvements, etc.

Note

If you would like to help developing some aspect of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\), or even implement your own modules with some new functionality we have not envisaged, please contact us and let us know about your idea(s). TempLat, which is developed by the \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) team, introduces a natural language describing fields and operations between them, so it is a natural platform to implement new libraries (related or not to cosmology).

Capabilities and structure

As explained in the previous Overview section, rather than being a conventional code designed to solve a predetermined set of equations and compute a fixed collection of observables, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) constitutes a platform for the implementation of field-theoretical systems governed by partial differential equations that can be discretized on a lattice. Written in C++, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) follows a modular architecture that cleanly separates the underlying physics from its technical implementation. It also introduces a dedicated symbolic language in which field variables and the operations acting on them are expressed in a form closely resembling their continuum counterparts. To be more specific, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is characterized by a number of aspects that we list below, which constitute all together a clear advantage for using \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) to investigate early Universe scenarios, over writing your own code from scratch:

Beyond singlet scalars

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is ready to simulate not only the evolution of singlet scalar theories, but also Abelian U(1) and non-Abelian SU(2) gauge theories. Since \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) \(\tt v2.0\) was released (on July 2026), one can also simulate scalars with a non-minimal coupling to gravity, \(\phi^2R\), and axion-like particle (ALP) interactions with gauge fields, \(\phi F_{\mu\nu}\tilde F^{\mu\nu}\); it incorporates methods to set up special field configurations, like cosmic defect networks in scaling (e.g. cosmic strings and domain walls), or arbitrary field power spectra; it incorporates non-symplectic evolution algorithms suitable e.g. for non-minimal scalar kinetic theories of the type \(\mathcal{G}_{ab}(\lbrace\phi_c\rbrace)\partial_\mu\phi^a\partial^\mu\phi^b\); it also allows for optimized gravitational wave dynamics on the lattice sourced by scalar and gauge fields, as well as for scalar field dynamics in \(d + 1\) dimensions with \(d = 1, 2\). All simulations can be done either in a flat space-time background, or in an expanding background.

TempLat, the lattice engine

The lattice engine of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is TempLat, which has been designed as a platform to implement any system of dynamical equations suitable for discretization on a lattice. \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is therefore not just a code for one type of simulation, but rather a more evolved concept. Via TempLat, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) uses a symbolic language that incorporates real, complex and SU(2) algebras, allowing to use vectorial and matrix notations without sacrificing performance. Once the user becomes familiar with the basic 'vocabulary' of the new language, they can modify \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) for their own purpose: let it be for the time evolution of new field variables in a given relevant model of interest, or for some other operation, like e.g. a Monte-Carlo generator for thermal configurations.

CosmoInterface, the physics library

The heart of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is CosmoInterface, which handles all relevant aspects of the physical theories, such as the field content, initialization, evolution equations, or relevant output observables. This makes the physics part of the code easy to understand and well separated from technical details from TempLat. This separation significantly simplifies the process of writing new operations for your own purposes. CosmoInterface represents essentially a library where \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) includes all the necessary ingredients that make possible to simulate a rich variety of scenarios. At the time of writing (July 2026), \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) \({\tt v2.0}\) is ready to simulate scenarios that include interacting singlet scalar fields; charged scalar fields under a U(1) and/or SU(2) gauge symmetries, and their corresponding Abelian and/or non-Abelian gauge vector fields; axion-like-particles (ALP) interacting with gauge sectors, non-minimal scalar field dynamics, and soon enough fluid dynamics and other interactions (see above Upcoming Physics). Simulations can be done in a flat space-time background, or in a spatially flat expanding FLRW background. In the latter case the fields can evolve either over an external background (e.g. with a power-law scale factor), or self-consistently, i.e. 'dictating' themselves the expansion of the universe as sourced by their volume averaged energy and pressure densities.

Lattice loops written for you

One of the main advantages of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is that it clearly separates the physics (e.g. definition of the field content, operations between fields, evolution equations, etc) in CosmoInterface, from the more technical details (e.g. parallelization aspects, Fourier transforms, etc) in TempLat [which is actually not stored inside the \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) source tree, as it is an external dependency fetched upon compilation of the code]. TempLat automates tecnical tasks fully, and many user will not never need to edit it. For example, if we have two fields f and g on a lattice and we want to sum them. Typically, one would explicitly write a loop that sums the amplitudes of both fields at each node of the lattice. Instead, in \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) we can just write f + g, and TempLat handles the whole operation of summing their values everywhere on the lattice. Essentially, all the technicalities such as memory handling or parallelization tasks, remain mostly hidden from a typical user.

From beginner to core developer

As a consequence of the above separation between physics and technicalities, a beginner user, say with little experience in programming and no experience at all in parallelization techniques, will be able to run a fully parallelized simulation of their favourite model (perhaps using thousands of processors in a cluster), while being completely oblivious to the technical details. They will just need to write a basic model file in the language of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\), containing the details of the model being simulated. If, on the contrary, the user is rather an experienced programmer and wants to look inside the core routines of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) and TempLat, and modify, say the MPI-implementation, they can always do so, and perhaps even contribute to their improvement.

High-order integrators, parallel FFTs

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) provides symplectic integrators in CosmoInterface, with accuracy ranging from \(\mathcal{O}(\delta t^2)\) up to \(\mathcal{O}(\delta t^{10})\), to simulate the non-linear dynamics of the appropriate fields in comoving two- or three-dimensional lattices. Appropriate observables are also provided for each case, like the energy density components of each field, their relevant spectra, or dynamical constraints. Our algorithms conserve energy up to the accuracy set by the order of the evolution algorithm, reaching down to machine precision in the case of the highest order integrators. Algorithms for gauge theories, either Abelian or non-Abelian, respect (independently of the integrator) the Gauss constraint to machine precision, even when considering self-consistent expansion. Moreover, through TempLat \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) uses a discrete Fourier Transform parallelized in multiple spatial dimensions ², this enables the code to run simulations with very high resolution, very long time scales, or well-separated scales.

The basic folder tree structure of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) is the following:

Folder treeWhat each folder contains

Click any node to open it on GitHub; dashed nodes are generated locally at build time.

with the content of each folder summarized as:

folder name(s)	brief description of each folder
`CMakeLists.txt`	top-level build configuration; fetches `TempLat`, selects the model, creates the executable, and links `TempLat::TempLat` with `CosmoInterface`
`cmake`	CosmoLattice-specific CMake helper files
`source`	executable entry point, currently `cosmolattice.cpp`
`include/CosmoInterface`	CosmoLattice physics interface: model abstraction, field definitions, evolution algorithms, initial conditions, measurements, and simulation management
`abstractmodel`	base classes and shared data structures for scalar, gauge, scale-factor, gravitational-wave, and non-minimal-coupling model sectors
`definitions`	reusable physical definitions such as potentials, gauge derivatives, Gauss laws, power-spectrum helpers, and anisotropic-stress tensors
`definitions/defectsmodule`	cosmic-defect machinery: defect observables, diffusion, resolution-preserving evolution, and the associated lattice weights
`evolvers`	evolution algorithms and evolution kernels
`initializers`	initialization for fields, scale factor and model state
`measurements`	observables, spectra, energy outputs, gravitational-wave outputs, and measurement I/O helpers
`measurements/measurementsIO`	output back-ends: `std` writes the plain-text measurement and spectra files, `hdf5` writes the same data as HDF5
`models`	user-facing model files selected with `-DMODEL=...`
`models/parameter-files`	example input parameter files for the implemented models
`tests`	CosmoInterface tests and benchmarks, enabled with `-DCOSMOINTERFACE_TEST=ON`
`profile`	profiling helpers: `gprof`/`massif` scripts, a sample input file, and a standalone thread-scaling executable
`build`	local build directory generated by CMake; contains object files and fetched dependency sources under `_deps/`. This directory is created by the user.
`build/_deps/templat-src`	generated checkout of the external `TempLat` repository, which provides the lattice language and technical infrastructure
`build/_deps/kokkos-src`	generated Kokkos source tree used by `TempLat` as the default device backend
`build/_deps/parafaft-src`	generated optional ParaFaFT source tree for MPI parallel Fourier transforms
`build/_deps/...`	generated dependency trees for optional packages such as KokkosFFT, FFTW, or HDF5, depending on the build options

The CMake configuration first selects the model to compile from models/, then builds the executable from source/cosmolattice.cpp. The executable includes CosmoInterface from this repository and links against the external TempLat::TempLat target. TempLat is fetched at configure time and is responsible for the lattice-expression machinery, memory layout, parameter parsing, Fourier-transform infrastructure, and device backend setup. Depending on the selected CMake options, TempLat may in turn fetch or find external packages such as Kokkos, KokkosFFT, ParaFaFT, FFTW, HDF5, MPI, OpenMP, or C++ threads.

Note

A remarkable feature of TempLat is that its operations and parallelization can work in an arbitrary number of spatial dimensions \(d\), including \(d < 3\) and \(d > 3\). This makes it a natural basis for future interfaces dealing with field dynamics on lower- or higher-dimensional lattices. In the latest version of the code, \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) \(\tt v2.0\), we have exploited this feature to run scalar field simulations in \(d = 1\) or \(2\) spatial dimensions. Keep an eye into the \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) News tab in case we eventually consider other cases.

Basic Field Equations implemented (so far)

We introduce now the equations of motion (EOM) that the latest version of \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) can solve. First we review the FLRW solution describing an expanding Universe, and then introduce the EOM of canonical scalar-gauge theories, various non-canonical theories, and gravitational waves, all embedded in a spatially-flat FLRW background. The derivation of the EOM can be found in our monographic reviews: \(\mathtt{The~Art-I}\) (Ref. ¹), and \(\mathtt{The~Art-II}\) (Ref. ³), see Lattice-Cosmology Reviews.

FLRW Expanding Background

In \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\), we consider all fields living in a spatially-flat homogeneous and isotropic spacetime background. This is described by the Friedmann–Lemaître–Robertson–Walker (FLRW) solution, characterized by the line element

\[ \begin{equation} d s^2 \equiv g_{\mu\nu}d x^\mu d x^\nu = - a^{2\alpha}d \eta^2 + a^2(\eta) \delta_{ij} d x^i d x^j \ , \label{eq_FLRWlineElem} \end{equation} \]

where \(a(\eta)\) is the scale factor, and \(\eta\) the so-called \(\alpha\)-time. The parameter \(\alpha\) is a constant that can be fixed conveniently depending on the scenario. The choice \(\alpha = 0\) identifies \(\eta\) with the coordinate time \(t\), whereas \(\alpha = 1\) makes \(\eta\) the conformal time \(\tau \equiv \int {dt'\over a(t')}\). In general, we will keep \(\alpha\) as an unspecified real number, writing all equations in \(\alpha\)-time. Only when a physical problem is chosen, then one needs to make a concrete choice for \(\alpha\).

The evolution of the scale factor is dictated by the first and second Friedmann equations (Ref. ¹)

\[ \begin{eqnarray} \mathcal{H}^2 \equiv \left({a'\over a}\right)^2 = a^{2 \alpha} \frac{\bar {\rho}}{3 m_p^2} \,,~~~~~~~ {a''\over a} = \frac{a^{2 \alpha}}{6 m_p^2}\Big[ (2 \alpha - 1) \bar{\rho} - 3 \bar{ p} \Big]\,, \label{eq_Friedmann-full} \end{eqnarray} \]

where \(\bar{p}\) and \(\bar{\rho}\) are the background values of the total pressure and energy densities of the field sectors, we use the notation \({f}' \equiv \text{d} f/ \text{d} \eta\), and \(m_p \simeq 2.435\cdot10^{18}\) GeV is the reduced Planck mass.

In \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\), all field sectors are considered to live by default in the above spatially-flat FLRW background. The user can either choose between an external expanding background, characterized by a constant equation of state \(w\), or self-consistent expansion of the Universe, solving the Friedmann Eqs. (\ref{eq_Friedmann-full}) as sourced with the volume-avaraged pressure and energy densities of the field sectors, \(\bar{p}\) and \(\bar{\rho}\). Alternatively, the user can choose instead to enforce the matter field sectors to live in flat-spacetime, i.e. using the Minkowski metric \(\eta_{\mu\nu} = (-,+,+,+)\).

Canonical Scalar-Gauge Theories

In canonical scalar-gauge theories we consider three types of scalar fields: a real scalar singlet \(\phi\), a U(1)-charged complex scalar field \(\varphi\), and a \([SU(2) \times U(1)]\)-charged doublet scalar field \(\Phi\). In the last two cases, there are also present, respectively, Abelian and non-Abelian gauge fields, \(A_{\mu}\) and \(C_{\mu} \equiv C_{\mu}^a T_a\), with \(T_a \equiv \sigma_a / 2\) the SU(2) group generators, where \(a=1,2,3\) and \(\sigma_a\) are the Pauli matrices. Like any \(SU(N)\) Lie algebra, the generators satisfy that \([T_a, T_b] = i f_{abc} T_c\) , \({\rm Tr}(T_a) = 0\) , \({\rm Tr}(T_a T_b) = \frac{1}{2} \delta_{ab}\), and \(T^{\dagger}_a = T_a\), with \([\cdot,\cdot]\) the usual matrix commutator, and \(f_{abc}\) the totally anti-symmetric structure constants, which in SU(2) are the totally anti-symmetric tensor, \(f_{abc} = \epsilon_{abc}\). For further details see Chapter 2 of \(\mathtt{The~Art-I}\) (Ref. ¹). The scalar fields can be explicitly written in terms of real components as follows,

\[ \begin{eqnarray} \label{eq_ChargedScalars} \begin{array}{ccccc} \phi \in \mathcal{R}e & , & \varphi \equiv {1\over\sqrt{2}}(\varphi_1 +i\varphi_2) & , & \Phi \equiv \begin{pmatrix} \varphi^{(1)} \\[1mm] \varphi^{(2)} \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} \varphi_1+i\varphi_2 \\[1mm] \varphi_3+i\varphi_4 \end{pmatrix}\,. \end{array} \end{eqnarray} \]

More specifically, for canonically-normalized scalar-gauge theories, we consider an action \(S = \int d^4x \sqrt{-g}\, \mathcal{L}\), with \(g \equiv {\rm det} (g_{\mu \nu})\) and Lagrangian

\[ \begin{align} -\mathcal{L} = \frac{1}{2}\partial_{\mu} \phi \partial ^{\mu}\phi + (D_{\mu}^A \varphi)^{*}(D_A^{\mu} \varphi) + (D_{\mu}\Phi )^{\dagger} (D^{\mu} \Phi) + \frac{1}{4} F_{\mu \nu} F^{\mu \nu} + \frac{1}{2}{\rm Tr}\{G_{\mu \nu}G^{\mu \nu}\} + V \ , \label{eq_lagrangian} \end{align} \]

where \(V \equiv V(\phi,|\varphi|, |\Phi|)\) is the potential describing the interactions between the scalar fields. The covariant derivatives and field strength tensors associated to the gauge fields, are defined as

\[ \begin{eqnarray} D_{\mu}^{\rm A} &\equiv & \partial _{\mu} - i g_A Q_AA_\mu \ , \hspace{4.5cm} F_{\mu \nu}\equiv \partial_{\mu} A_{\nu} - \partial_{\nu} A_{\mu} \ , \label{eq_AbCovDerivCont} \\ D_{\mu} & \equiv & \mathcal{I}D^{\rm A}_\mu - i g_C Q_C C_{\mu}^a \,T_a \ , \hspace{2.76cm} G_{\mu \nu} \equiv \partial_{\mu} C_{\nu} - \partial_{\nu} C_{\mu} - i[C_\mu,C_\nu]\,, \label{eq_CovDerivCont} \end{eqnarray} \]

with \(g_{A}\) and \(g_C\) the Abelian and non-Abelian gauge couplings, \(Q_{A}\) and \(Q_C\) the Abelian and non-Abelian charges of the scalar fields, \(\mathcal{I}\) the \(N\times N\) identity matrix. The gauge-invariant electric and magnetic fields associated to the Abelian and non-Abelian fields can be written as

\[ \begin{equation}\label{eq_ElectricMagneticDefs} E_i \equiv F_{0i} , \,\,\,\,\,\,\,\, B_i \equiv \frac{1}{2} \epsilon_{i j k} F^{j k} , \,\,\,\,\,\,\,\, E_i^a \equiv G_{0i}^a , \,\,\,\,\,\,\,\, B_i^a \equiv \frac{1}{2} \epsilon_{i j k} G^{j k}_a \ , \end{equation} \]

where \(\epsilon_{ijk}\) is the Levi-Civita symbol in three dimensions with normalization \(\epsilon_{123}=+1\), and \(G_{\mu \nu}^a \equiv {\rm Tr}(2G_{\mu \nu} T_a) = {\rm Tr}(G_{\mu \nu} \sigma_a)\). We note that the electric field definitions above depend on the \(\alpha\)-time \(\eta\), as \(F_{0i}\) and \(G_{0i}\) are defined with respect to \(\eta\), not \(t\).

The equations of motion for the matter fields and the scale factor have been derived in more detail in \(\mathtt{The~Art-I}\) (Ref. ¹). Here we simply quote their resulting form, which read

\[ \begin{eqnarray} \phi'' - a^{-2(1 - \alpha)} {\vec\nabla}^{\,2} \hspace{-1mm}\phi + (3 - \alpha)\mathcal{H} {\phi'} &=& - a^{2 \alpha} V_{,\phi} \ , \label{eq_singlet-eomCONT} \\ \varphi'' - a^{-2(1 - \alpha)} {\vec D}_{\hspace{-0.5mm}A}^{\,2}\varphi + (3 - \alpha) \mathcal{H} {\varphi'} &=& - \frac{a^{2 \alpha}V_{,|\varphi|} }{2} \frac{\varphi}{|\varphi |} \ , \label{eq_higgsU1-eom}\\ \Phi'' - a^{-2(1 - \alpha)} {\vec D}^{\,2}\Phi + (3 - \alpha) \mathcal{H} {\Phi'} &=& - \frac{a^{2 \alpha} V_{,|\Phi|}}{2} \frac{\Phi}{|\Phi |} \ , \label{eq_higgsSU2-eom} \\ \partial_0 F_{0i} - a^{-2(1 - \alpha )}\partial_j F_{ji} + (1 - \alpha) \mathcal{H} F_{0i} &=& a^{2 \alpha}J^A_i \ , \label{eq_U1eom} \\[2em] (\mathcal{D}_0 )_{a b} (G_{0i})^b - a^{-2(1 - \alpha )} ( \mathcal{D}_j )_{a b} (G_{ji} )^b + (1 - \alpha) \mathcal{H} (G_{0i} )^b &=& a^{2 \alpha}(J_i)_a \ , \label{eq_SU2eom} \end{eqnarray} \]

where we have introduced the derivative operator \((\mathcal{D}_{\nu}O)_a = (\mathcal{D}_{\nu})_{a b}O_b \equiv ( \delta_{a b} \partial_{\nu} - f_{abc} C_{\nu}^c ) O_b\), and defined the matter currents

\[ \begin{eqnarray} \label{eq_AbelianCurrent} J_A^\mu & \equiv & 2g_A Q_A^{(\varphi)} \mathcal{I}m [ \varphi^{*} ( D_A^{\mu} \varphi )] + 2g_A Q_A^{(\Phi)} \mathcal{I}m [ \Phi^\dagger (D^{\mu} \Phi )]\,,\\ \label{eq_NonAbelianCurrent} J_a^\mu & \equiv & 2g_C Q_C\mathcal{I}m [ \Phi^{\dagger} T_a( D^{\mu} \Phi )]\,. \end{eqnarray} \]

We note that the Gauss constraints of the Abelian and non-Abelian sectors, respectively, must be preserved at all times during the evolution:

\[ \begin{eqnarray} \partial_i F_{0i} &=& a^2J^A_0 \ , \label{eq_GaussU1-eom}\\ (\mathcal{D}_i )_{a b} (G_{0i})^b &=& a^2(J_0)_a \ . \label{eq_GaussSU2-eom} \end{eqnarray} \]

If the fields dominate the energy budget of the Universe, the expansion rate can be determined through the Friedmann equations, which in our case can be written as

\[ \begin{align} \label{eq_FriedmannHub} \mathcal{H}^2 \equiv \frac{a'^{ 2}}{a^2} &= \frac{a^{2 \alpha}}{3} \left( \frac{f_*}{m_p}\right)^2 \left[ E_K^{\phi} + E_K^{\varphi} + E_K^{\Phi} + E_G^{\phi} + E_G^{\varphi} + E_G^{\Phi} + E_K^A + E_K^B + E_G^A + E_G^B + E_V \right] , \\ \label{eq_FriedmannD2a} {a''\over a} &= \frac{a^{2 \alpha}}{3} \left( \frac{f_*}{m_p}\right)^2 \left[ (\alpha-2)(E_K^{\phi} + E_K^{\varphi} + E_K^{\Phi}) + \alpha(E_G^{\phi} + E_G^{\varphi} + E_G^{\Phi}) + (\alpha + 1)E_V \right.\\ & \hspace{2cm} \left. + (\alpha-1)(E_K^A + E_K^B + E_G^A + E_G^B) \right] ,\nonumber \end{align} \]

where \(\langle \dots \rangle\) denotes an average over sufficiently large volumes that encompass all relevant wavelengths of the fields, and we have defined the volume-averaged energy contributions as \(E_{K}^{f} = \langle K_{f} \rangle\) and \(E_{G}^{f} = \langle G_{f} \rangle\) for the scalar fields \(f=\phi,\varphi,\Phi\), \(E_{K}^{A} = \langle K_{U(1)} \rangle\), \(E_{G}^{A} = \langle G_{U(1)} \rangle\), \(E_{K}^{B} = \langle K_{SU(2)} \rangle\), and \(E_{G}^{B} = \langle G_{SU(2)} \rangle\) for the gauge fields, and \({E}_V = \langle {V} \rangle\) for the potential energy. The different energy density contributions are given by

\[ \begin{align} \label{eq_energy-contributions} \hspace{-2cm}\left\lbrace \begin{array}{rcl} {K}_{\phi} &=& \frac{1}{2 a^{2\alpha} } \phi'^2 \\ {K}_{\varphi} &=& \frac{1}{a^{2\alpha} } (D_0^A \varphi)^*(D_0^A \varphi) \\ {K}_{\Phi} &=& \frac{1}{a^{2\alpha} } (D_0 \Phi )^\dagger(D_0 \Phi) \\ \end{array}\right. \hspace{0.1cm};\hspace{0.75cm} \left\lbrace \begin{array}{rcl} {G}_{\phi} &=& \frac{1}{2 a^2} \sum_i (\partial_i \phi)^2 \\ {G}_{\varphi} &=& \frac{1}{a^2} \sum_i (D_i^A \varphi)^*(D_i^A \varphi) \\ {G}_{\Phi} &=& \frac{1}{a^2} \sum_i (D_i\Phi)^\dagger(D_i \Phi) \\ \end{array}\right. \hspace{0.1cm};\hspace{0.75cm} \left\lbrace \begin{array}{rcl} {K}_{U(1)} &=& \frac{1}{2 a^{2 + 2 \alpha}} \sum_{i} F_{0i}^2 \\ {K}_{SU(2)} &=& \frac{1}{2 a^{2 + 2 \alpha}} \sum_{a,i} (G_{0i}^a)^2 \\ {G}_{U(1)} &=& \frac{1}{2 a^4} \sum_{i,j < i} F_{ij}^2 \\ {G}_{SU(2)} &=& \frac{1}{2 a^4} \sum_{a,i,j < i} (G_{ij}^a)^2 . \\ \end{array}\right. \\\nonumber\\ \text{(Kinetic-Scalar)} \hspace{5cm} \text{(Gradient-Scalar)} \hspace{6.5cm} \text{(Electric & Magnetic)} \hspace{3.0cm}\nonumber\\ \end{align} \]

In \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) we use Eq. (\ref{eq_FriedmannD2a}) to solve for the scale factor, while monitoring that the Hubble constraint (\ref{eq_FriedmannHub}) is verified throughout the evolution to some desired accuracy, see Section Evolution Algorithms.

Non-Canonical Field Theories

Among the possible field theories with non-canonical interactions, we consider scalar fields non-minimally coupled (NMC) to gravity through a term of the form \(\propto\phi^2R\), with \(R\) the Ricci scalar, and axion-like particles (ALPs) interacting with gauge fields through a coupling of the form \(\phi F\tilde F\). Further details on these cases can be found in Chapters 3 and 4 of \(\mathtt{The~Art-II}\) (Ref. ³).

Scalar Non-Minimally Coupled to Gravity

We consider a theory containing a curvature interaction proportional to \(\xi R\phi^2\) and a generic potential \(V(\phi,\{\varphi_{\rm m}\})\) for \(\phi\) and the remaining matter fields \(\{\varphi_{\rm m}\}\). In the Jordan-frame, restricting the metric to a spatially flat FLRW background, and hence neglecting gravitational perturbations, the equation of motion for \(\phi\) becomes

\[ \begin{eqnarray} \phi'' +(3-\alpha)\frac{a'}{a}\phi' -a^{-2(1-\alpha)}\nabla^2\phi +a^{2\alpha} \left( \xi\bar R\phi+\frac{\partial V}{\partial\phi} \right) =0\,, \label{eq_eom}\\ {\rm where} ~~~~~~\bar R = \frac{6}{a^{2\alpha}} \left[ \frac{a''}{a} +(1-\alpha)\left(\frac{a'}{a}\right)^2 \right]\,.~~~~~~~~~~~~~~ \label{eq_cosmic_R} \end{eqnarray} \]

where \(R\) is the backgorund Ricci scalar, and \(\xi\) the non-minimal coupling parameter.

The Friedmann equations in \(\alpha\)-time read

\[ \begin{eqnarray} \mathcal{H}^{2} \equiv \left(\frac{a'}{a}\right)^2 = \frac{a^{2\alpha}}{3m_p^2} \left(\bar\rho_\phi+\bar\rho_{\rm m}\right) \label{eq_Hu}\,~~;~~~~~~~~~~ \frac{a''}{a} = -\frac{a^{2\alpha}}{6m_p^2} \left[ (1-2\alpha) \left(\bar\rho_\phi+\bar\rho_{\rm m}\right) +3\left(\bar p_\phi+\bar p_{\rm m}\right) \right]\,, \end{eqnarray} \]

where the energy density and pressure of the non-minimally coupled scalar field are ⁴

\[ \bar\rho_\phi(\eta) = \frac{1}{2a^{2\alpha}}\left\langle\phi'^2\right\rangle +\frac{1}{2a^2}\left\langle(\nabla\phi)^2\right\rangle +\left\langle V(\phi)\right\rangle +\frac{3\xi}{a^{2\alpha}}\mathcal{H}^2 \left\langle\phi^2\right\rangle +\frac{6\xi}{a^{2\alpha}}\mathcal{H} \left\langle\phi\phi'\right\rangle -\frac{\xi}{a^2}\left\langle\nabla^2\phi^2\right\rangle\,. \label{eq_nmcrho}\\ \]

\[ \begin{aligned} \bar p_\phi(\eta) ={}& \frac{1-4\xi}{2a^{2\alpha}}\left\langle\phi'^2\right\rangle -\frac{1-12\xi}{6a^2}\left\langle(\nabla\phi)^2\right\rangle -\left\langle V(\phi)\right\rangle +\frac{2\xi}{a^{2\alpha}}\mathcal{H} \left\langle\phi\phi'\right\rangle -\frac{\xi}{3a^2}\left\langle\nabla^2\phi^2\right\rangle \\ &+ 2\xi\left\langle\phi V_{,\phi}\right\rangle +\frac{\xi}{a^{2\alpha}} \left[ \mathcal{H}^2 +12\left(\xi-\frac{1}{6}\right) \left( \frac{a''}{a} +(1-\alpha)\mathcal{H}^2 \right) \right] \left\langle\phi^2\right\rangle\,, \end{aligned} \label{eq_nmcp} \]

with \(V_{,\phi}\equiv\partial V/\partial\phi\), and \(\langle \dots \rangle\) denoting volume-averaging over sufficiently large scales that encompass all relevant wavelengths of the fields. The scale factor may, in principle, be evolved using the Friedmann equations in \(~\)\eqref{eq_Hu}. We note, however, that contrary to canonical scenarios, the r.h.s. of the Friedmann equations depend on time derivatives of \(a(\eta)\), given the expression of \(\bar\rho_\phi(\eta), \bar p_\phi(\eta)\).

Alternatively, the scale factor evolution can be obtained from the trace of the energy-momentum tensor of the non-minimally coupled field, \(T_\phi\). Given the traced Einstein equations \(R = -\frac{1}{m_p^2}g^{\mu\nu}\left(T^\phi_{\mu\nu}+T^{\rm m}_{\mu\nu}\right)\) = \(-\frac{1}{m_p^2}\left(T_\phi+T_{\rm m}\right)\), it follows that the background curvature satisfies (Ref. ⁴)

\[ m_p^2\bar R = (1-6\xi) \left[ \left\langle \partial^\mu\phi\,\partial_\mu\phi \right\rangle +\xi\bar R\left\langle\phi^2\right\rangle \right] -6\xi\left\langle\phi V_{,\phi}\right\rangle +4\left\langle V\right\rangle -\left\langle T_{\rm m}\right\rangle\,, \label{eq_EFEtrBack} \]

where \(\langle\cdots\rangle\) denotes again volume-averaging of the corresponding local spatial inhomogeneities, over length scales much larger than the inverse gradient-scales of the problem. Solving the above expression for \(\bar R\), then gives

\[ \begin{eqnarray} \bar R = \frac{F(\phi)}{m_p^2} \left[ (1-6\xi) \left\langle \partial^\mu\phi\,\partial_\mu\phi \right\rangle +4\left\langle V\right\rangle -6\xi\left\langle\phi V_{,\phi}\right\rangle -\left\langle T_{\rm m}\right\rangle \right]\,. \label{eq_eomR}\\ ~~~~{\rm where}~~~ F(\phi) \equiv \frac{1}{ 1+(6\xi-1)\xi\left\langle\phi^2\right\rangle/m_p^2 }\,.~~~~~~~~~~~~~~~ \label{eq_Fphi} \end{eqnarray} \]

Using Eq.\(~\)\eqref{eq_cosmic_R}, we then obtain the following differential equation for the scale factor,

\[ \frac{a''}{a} +(1-\alpha)\left(\frac{a'}{a}\right)^2 = \frac{a^{2\alpha}F(\phi)}{6m_p^2} \left[ (1-6\xi) \left\langle \partial^\mu\phi\,\partial_\mu\phi \right\rangle +4\left\langle V\right\rangle -6\xi\left\langle\phi V_{,\phi}\right\rangle -\left\langle T_{\rm m}\right\rangle \right]\,. \label{eq_piadot} \]

This equation can be evolved simultaneously with the equations of motion of the non-minimally coupled scalar and the remaining matter fields. In \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) we actually follow this method.

Axion-Gauge interactions

Axion-like particles (ALPs) enjoy a shift-symmetry \(\phi=\phi+C\) that allows them to couple derivatively to gauge fields through Chern–Simons terms, as \(\phi F\tilde{F}\) (Abelian) or \(\phi G\tilde{G}\) (non-Abelian). Since \(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) \(\tt v2.0\), the code considers an axion–\(U(1)\) sector described by the action

\[ S = \int d^4x\,\sqrt{-g} \left[ \frac{1}{2}m_p^2R -\frac{1}{2}\partial_\mu\phi\,\partial^\mu\phi -V(\phi) -\frac{1}{4}F_{\mu\nu}F^{\mu\nu} +\frac{1}{4}\frac{\phi}{\Lambda} F_{\mu\nu}\tilde F^{\mu\nu} \right]\,, \label{eq_AxionAction} \]

where \(\phi\) is a real pseudo-scalar field, the axion-like particle (ALP), \(V(\phi)\) is its interaction potential (which might break the shift symmetry explicitly in some scenarios), and the field-strength tensor of the Abelian gauge field \(A_\mu\) is defined as in standard canonical theories, \(F_{\mu\nu}\equiv\partial_\mu A_\nu-\partial_\nu A_\mu\). To characterize the ALP-gauge coupling strength, we define a dimensionless parameter \(\alpha_\Lambda\equiv m_p/\Lambda\) associated to the scale \(\Lambda\). The dual field-strength tensor is defined as \(\tilde F_{\mu\nu}\equiv\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}\), where \(\epsilon_{\mu\nu\rho\sigma}\) is the four-dimensional Levi-Civita pseudotensor in curved spacetime, normalized as \(\epsilon_{0123}=1/\sqrt{-g}\). Varying the action with respect to \(\phi\) and \(A_i\) in a FLRW background, and considering the temporal gauge, \(A_0=0\), leads to the equations of motion in \(\alpha\)-time as

\[ \begin{eqnarray} \phi'' +(3-\alpha)\mathcal{H}\phi' -a^{2(\alpha-1)}\vec{\nabla}^{\,2}\phi +a^{2\alpha}V_{,\phi} &=& \frac{\alpha_\Lambda}{m_p} a^{\alpha-3}\vec E\cdot\vec B\,, \label{eq_axion_eom} \\ E_i' +(1-\alpha)\mathcal{H}E_i +a^{2(\alpha-1)}\epsilon_{ijk}\partial_jB_k &=& -\frac{\alpha_\Lambda}{m_p}a^{\alpha-1} \left( \phi'B_i-\epsilon_{ijk}\partial_j\phi\,E_k \right)\,. \label{eq_axion_gauge_eom} \end{eqnarray} \]

where primes denote derivatives with respect to \(\alpha\)-time, \(\eta\), and \(E_i \equiv F_{0i} = A_i'\). The system is subject to the Gauss constraint, which, in the absence of external charged currents, reads as

\[ \partial_iE_i = -\frac{\alpha_\Lambda}{m_p} a^{\alpha-1}\partial_i\phi\,B_i\,. \label{eq_axion_gauss} \]

The local energy density and pressure of the combined axion–gauge sector are

\[ \begin{eqnarray} \rho &=& K_\phi+G_\phi+V(\phi)+K_{U(1)}+G_{U(1)}\,, \label{eq_axion_rho} \\ p &=& K_\phi-\frac{1}{3}G_\phi-V(\phi) +\frac{1}{3} \left[ K_{U(1)}+G_{U(1)} \right]\,. \label{eq_axion_p} \end{eqnarray} \]

The individual kinetic and gradient contributions of the scalar and vector fields coincide with those defined in Eq.\(~\)\eqref{eq_energy-contributions}. If the axion–\(U(1)\) gauge sector dominates the energy budget of the Universe, the expansion is governed by the Friedmann equations, Eqs.\(~\)\eqref{eq_FriedmannHub} and \eqref{eq_FriedmannD2a}, retaining only the contributions from this sector.

Gravitational Waves

\(\mathcal{C}\mathtt{osmo}\mathcal{L}\mathtt{attice}\) can compute the generation of gravitational waves (GWs) sourced by singlet-scalar, Abelian gauge-scalar, and Abelian axion-gauge sectors. Technically, GWs are transverse and traceless tensor perturbations \(h_{ij}\) of a spacetime metric \(g_{ij}\). Over the FLRW background, the perturbed line element in \(\alpha\)-time reads

\[ \begin{align} \label{eq_GWmetric} ds^2 = -a^{2\alpha}(\eta)d\eta^2 +a^2(\eta)\left(\delta_{ij}+h_{ij}\right)dx^idx^j\,, \end{align} \]

with \(\partial_i h_{ij}=0\) (transverse) and \(h_{ii}=0\) (traceless) Linearizing the Einstein equations in \(h_{ij}\), leads to the equation of motion of the GWs as

\[ \begin{align} \label{eq_GWEOMcontinuum} h_{ij}'' +(3-\alpha){a'\over a}h_{ij}' -a^{-2(1-\alpha)}\nabla^2h_{ij} = {2\over m_p^2a^{2(1-\alpha)}}\Pi_{ij}^{\rm TT} \,. \end{align} \]

GWs are sourced by the transverse-traceless (TT) part of the anisotropic stress tensor, \(\Pi_{ij}^{\rm TT}\). For a generic background fluid, this anisotropic stress tensor takes the form \(\Pi_{ij} \equiv T_{ij}-\bar p\,g_{ij}\). In practice, however, it is more convenient to define an effective anisotropic stress tensor containing only those contributions to \(\Pi_{ij}\) that have a non-zero TT projection. For example, for a generic model consisting of canonically normalized scalars and Abelian gauge fields, this takes the form (Ref. ⁵)

\[ \begin{align} \label{eq_GWEffectiveAnisotropicStress} \Pi_{ij}^\mathrm{eff}=\sum_a \nabla_i\phi_a \nabla_j\phi_a+2\sum_b \text{Re}\left[ \left(D_i^A\varphi_b\right)^*D_j^A\varphi_b \right]-a^{-2\alpha}E_iE_j-a^{-2}B_i B_j \,, \end{align} \]

from which the source of GWs is obtained after TT projection. This operation is non-local in real space, but corresponds to an algebraic relation in Fourier space,

\[ \begin{align} \label{eq_GWTTsourceProjection} \Pi_{ij}^{\rm TT}({\bf k},\eta) = \Lambda_{ij,lm}(\hat{\bf k})\Pi_{lm}^{\rm eff}({\bf k},\eta)\,, \end{align} \]

where the Fourier-space TT projector takes the form,

\[ \begin{align} \label{eq_GWTTprojector} \Lambda_{ij,lm}(\hat{\bf k}) \equiv P_{il}(\hat{\bf k})P_{jm}(\hat{\bf k}) -{1\over2}P_{ij}(\hat{\bf k})P_{lm}(\hat{\bf k})\,,\quad\quad \text{with}\quad\quad P_{ij}(\hat{\bf k}) \equiv \delta_{ij}-\hat k_i\hat k_j\,, \end{align} \]

and \(\hat k_i\equiv {k_i\over k}\). This guarantees the projected tensor is both transverse, \(k_i\Pi_{ij}^{\rm TT}=0\), and traceless, \(\Pi_{ii}^{\rm TT}=0\).

Finally, it is worth mentioning GW observables. The most relevant quantity related to GWs is the energy density of the GW background,

\[ \begin{align} \label{eq_GWrhoContinuum} \rho_{\rm GW}(\eta) &= {m_p^2\over 4a^{2\alpha}V} \int_V d^3{\bf x}\, h'_{ij}({\bf x},\eta)h'_{ij}({\bf x},\eta) \simeq {m_p^2\over 4a^{2\alpha}V} \int_V {d^3{\bf k}\over(2\pi)^3} h'_{ij}({\bf k},\eta)h_{ij}^{\prime *}({\bf k},\eta) \nonumber \\ &\equiv \int \left({d\rho_{\rm GW}\over d\log k}d\log k\right) \,, \end{align} \]

from which the spectral density can be defined as

\[ \begin{align} \label{eq_GWrhoSpectrumContinuum} \left({d\rho_{\rm GW}\over d\log k}\right)(k,\eta) = {m_p^2 k^3\over 8\pi^2a^{2\alpha}V} \int {d\Omega_k\over4\pi}\, h'_{ij}(\hat{\bf k},k,\eta)h_{ij}^{\prime *}(\hat{\bf k},k,\eta)\,. \end{align} \]

For stochastic sources the volume average can be replaced by an ensemble average \(\langle...\rangle\) over the independent realizations of the tensor fluctuations,

\[ \begin{align} \rho_{\rm GW}(\eta) &= \dfrac{m_p^2}{4a^{2\alpha}} \left\langle h'_{ij}({\bf x},\eta) h_{ij}^{\prime *}({\bf x},\eta)\right\rangle = \dfrac{m_p^2}{4a^{2\alpha}} \int \dfrac{\text{d}^3{\bf k}}{(2\pi)^3} \dfrac{\text{d}^3{\bf k'}}{(2\pi)^3} e^{-i {\bf x}\cdot({\bf k} - {\bf k'})} \left\langle h'_{ij}({\bf k},\eta) h_{ij}^{\prime *}({\bf k'},\eta)\right\rangle \nonumber \\ &\equiv \dfrac{m_p^2}{8\pi^2a^{2\alpha}} \int\dfrac{\text{d}k}{k} k^3 P_{h'}(k,\eta)\,, \label{eq_stochasticPS} \end{align} \]

where we have introduced the power spectrum of the time derivative of \(h_{ij}\),

\[ \begin{align} \left\langle h'_{ij}({\bf k},\eta)h_{ij}^{\prime *}({\bf k'},\eta) \right\rangle = (2\pi)^3 P_{h'}(k,\eta)\delta^{(3)}({\bf k} - {\bf k'})\,. \label{eq_stochasticPS_2} \end{align} \]

In addition to the energy density, one usually also defines the fractional GW energy density power spectrum as

\[ \begin{align} \label{eq_GWOmegaContinuum} \Omega_{\rm GW}(k,\eta) \equiv {1\over\rho_\text{c}}{d\rho_{\rm GW}\over d\log k} = {k^3\over 24\pi^2\mathcal H^2}P_{h'}(k,\eta) = {k^3\over 24\pi^2a^{2\alpha}H^2}P_{h'}(k,\eta)\,. \end{align} \]

where \(\rho_\text{c}=3m_p^2 H^2\) is the critical energy density.

D. G. Figueroa, A. Florio, F. Torrenti, and W. Valkenburg. The art of simulating the early universe – part i: integration techniques and canonical cases. JCAP, 04:035, 2021. arXiv:2006.15122, doi:10.1088/1475-7516/2021/04/035. ↩↩↩↩↩↩
M. Pippig. Pfft: an extension of fftw to massively parallel architectures. SIAM J. Sci. Comput., 35:C213, 2013. ↩
J. Baeza-Ballesteros, D. G. Figueroa, A. Florio, J. Lizarraga, N. Loayza, K. Marschall, T. Opferkuch, B. A. Stefanek, F. Torrentí, and A. Urio. The art of simulating the early universe. part ii. non-canonical cases & gravitational waves. JCAP, 06:087, 2026. arXiv:2512.15627, doi:10.1088/1475-7516/2026/06/087. ↩↩
Daniel G. Figueroa, Adrien Florio, Toby Opferkuch, and Ben A. Stefanek. Lattice simulations of non-minimally coupled scalar fields in the jordan frame. SciPost Phys., 15(3):077, 2023. arXiv:2112.08388, doi:10.21468/SciPostPhys.15.3.077. ↩↩
Daniel G. Figueroa. Aspects of Reheating. PhD thesis, 2010. ↩