Tracking the Multifield Dynamics with Cosmological Data: A Monte Carlo approach

. We introduce a numerical method specifically designed for investigating generic multifield models of inflation where a number of scalar fields ϕ K are minimally coupled to gravity and live in a field space with a non-trivial metric G IJ ( ϕ K ) . Our algorithm consists of three main parts. Firstly, we solve the field equations through the entire inflationary period, deriving predictions for observable quantities such as the spectrum of scalar perturbations, primordial gravitational waves, and isocurvature modes. We also incorporate the transfer matrix formalism to track the behavior of adiabatic and isocurvature modes on super-horizon scales and the transfer of entropy to scalar modes after the horizon crossing. Secondly, we interface our algorithm with Boltzmann integrator codes to compute the subsequent full cosmology, including the cosmic microwave background anisotropies and polarization angular power spectra. Finally, we develop a novel sampling algorithm able to efficiently explore a large volume of the parameter space and identify a sub-region where theoretical predictions agree with observations. In this way, sampling over the initial conditions of the fields and the free parameters of the models, we enable Monte Carlo analysis of multifield scenarios. We test all the features of our approach by analyzing a specific model and deriving constraints on its free parameters. Our methodology provides a robust framework for studying multifield inflation, opening new avenues for future research in the field.


Introduction
Inflation is a period of rapid expansion in the early Universe initially proposed to account for various observational phenomena, including spatial flatness, the horizon and entropy problems, and the apparent lack of topological defects [1][2][3].However, it became quickly clear that inflation also offers an elegant mechanism to explain the physical origins of the first fluctuations in the Universe, which eventually gave rise to the observed structures such as galaxies and clusters of galaxies [4][5][6][7].According to the theory of inflation, the initial seeds are generated by quantum fluctuations in the inflaton field (or fields, if inflation is driven by multiple fields) that are stretched to superhorizon scales during inflation, eventually sourcing density fluctuations in other matter species.
Due to its ability to account for the origin of the structures observed in the presentday Universe, inflation is widely accepted as the leading theory of the very early Universe.However, despite its remarkable success, inflation is not free from limitations.From an observational perspective, the absence of a definitive detection of B-mode polarization [8,9] and the emerging discrepancies among different Cosmic Microwave Background (CMB) experiments [10][11][12][13][14][15][16][17][18] pose new challenges in determining precise predictions for the inflationary models/mechanisms that best explain observational data [19], and various studies suggest that modifications to the inflationary sector of the cosmological model might play a relevant role in addressing (part of) the tension between the value of the expansion rate of the Universe H 0 [20][21][22][23][24] measured through direct local distance ladder measurements and the value inferred from CMB observations, see e.g., [25][26][27][28][29][30][31][32][33][34].On the other hand, from a theoretical standpoint, despite a plethora of proposed models and mechanisms [35], the nature of the inflation field (or fields) remains still unknown and embedding inflation in a more fundamental theory remains an open problem [36][37][38].
The simplest dynamical models of inflation involve a single scalar field minimally coupled to gravity whose evolution should be governed by a potential enough flat to induce a phase of slow-roll evolution.However, several non-standard realizations of inflation have been proposed both in the context of extensions to the Standard Model of particle physics and in modified gravity theories and tested against a wide range of available data, including CMB, Big Bang Nucleosynthesis, and Gravitational Wave measurement. 1 Although a large portion of inflation models or theories is predominantly shaped by a single scalar field, low-energy effective field theories inspired by theories of particle physics beyond the Standard Model or quantum gravity, often incorporate multiple scalar degrees of freedom and suggest that inflation could be driven by multiple fields [108], potentially featuring non-minimal couplings [109][110][111][112][113][114][115][116][117][118][119][120][121][122][123][124].When inflation is driven by multiple scalar fields, after the inflationary period, they have to decay into various standard model particles such as dark matter, baryons, neutrinos, and other species.Furthermore, in multifield models, both adiabatic perturbations and isocurvature modes play crucial roles during inflation [125].These modes can persist immediately after the end of inflation and have a significant impact on the evolution of perturbations during the radiation-dominated epoch, giving rise to a rich phenomenology that can be tested and constrained with cosmological and astrophysical data.Therefore, it is important to extract clues about the period of inflation from cosmological observations to determine at least some of the properties of the field(s) driving inflation in the very early Universe [121,[124][125][126][127][128][129][130][131][132][133][134][135][136].In this regard, it is worth noting that both the amount of isocurvature and adiabatic modes can be accurately constrained by recent measurements of the anisotropies and polarization in the cosmic microwave background radiation [165][166][167][168][169][170][171], offering a valuable opportunity for experimental validation of multifield inflation.However, obtaining precise predictions from generic multifield models/theories is not always easy because observational quantities often depend on various factors.For instance, it is widely known that different initial conditions for the fields lead to different trajectories in field space which can produce slightly different results for observables such as the amplitude of the scalar and tensor perturbations and the spectral index for scalar perturbations, or isocurvature modes [172].This makes a comparison between theory and observations more challenging than in single-field inflation, and most tools employed in cosmological data analyses, such as typical Boltzmann integrator codes and samplers, are either unaware of the physics of inflation or assume single-field potentials. 2 As a result, constraining the multifield landscape in light of current observational data represents an ongoing challenge in the field.
In this paper, we take a first step to tackle down this difficulty and introduce a numerical method to precisely calculate the predictions resulting from generic multifield models of inflation where fields are minimally coupled to gravity and the field space metric is allowed to be non-trivial.To this end, our method is made of three key components.First, we numerically solve the complete field equations throughout the entire inflationary period.Once we have ensured that the fields undergo a phase of slow-roll evolution, we numerically solve the background dynamics by adopting a first-order slow-roll approximation and derive predictions for observable quantities such as the spectrum of scalar perturbations, primordial gravitational waves, and isocurvature modes.We also track the behavior of adiabatic and isocurvature modes on super-horizon scales and the transfer of entropy to scalar modes after crossing the horizon using transfer matrix formalism.Secondly, we interface our algorithm with well-known Boltzmann integrator codes and use the inflationary predictions as initial conditions to compute the subsequent full cosmology, including the CMB anisotropies and polarization angular power spectra.Finally, we narrow down the viable parameter space of the model and derive constraints on its free parameters by introducing a novel sampling algorithm designed to efficiently explore a large parameter volume and identify regions where predictions agree with observations.
The paper is organized as follows.In section 2, we give an overview of the method, starting with a description of the theoretical parameterization (subsection 2.1), discussing the numerical integration scheme adopted for tracking the multifield dynamics (subsection 2.2), and introducing the sampling algorithm and its interfacing with Boltzmann integrator codes such as CAMB [174,175] or CLASS [176] (subsection 2.3).In section 3, we apply this method to an example inflationary model, obtaining constraints on the parameters of the theory.Our conclusions are presented in section 4.

Probing the Multifield Landscape
In this section, we provide a comprehensive overview of the methodology developed for probing the multifield landscape of inflation.In subsection 2.1 we describe the theoretical parameterization adopted for the background equations and the formalism used to track the super-horizon evolution of perturbations.In subsection 2.2 we introduce our integration algorithm, explaining the various analyses performed to reconstruct the multifield dynamics throughout the full inflationary epoch, and the methodology used for calculating observables such as the spectra of primordial scalar and tensor perturbations, and the entropy transfer after horizon crossing.Finally, in subsection 2.3, we explain how our algorithm can be interfaced with standard Boltzmann integrator codes to calculate the subsequent full cosmology of the model and describe the sampling algorithm we designed to explore the parameter-space of generic multifield models and compare theoretical predictions with observations.

Parametrizing the Multifield Dynamics
We start by describing generic multifield models where a number of scalar fields ϕ K are minimally coupled to gravity and live in a field space with a non-trivial metric G IJ (ϕ K ).This geometry is reflected in the kinetic part of the Lagrangian Note that in this this paper we always work in units 8πG = 1.Considering a spatially flat Friedmann-Lemaître-Robertson-Walker metric, the inflationary dynamics can be described by the generalized Klein-Gordon equations obtained from the variation of the action with respect to the fields ϕ K : where g is the determinant of the metric tensor, V ≡ V (ϕ K ), and the notation , K indicates the derivative with respect to the fields.On the other hand, the evolution of the scale factor a(t) is governed by the Friedmann equations: To facilitate the interpretation of the evolution of linear cosmological perturbations, we adopt the formalism of Refs.[177,178] and perform a rotation in the field space.To do so, we define an orthonormal basis in the field space {e I n } (with n = 1, 2) as (where | φK | ≡ √ 2L kin denotes the length of the velocity vector φK containing the fields as components), and the decomposition of their perturbations [179] as In this way, we can introduce the comoving curvature perturbation [177,180] which encodes the adiabatic perturbations, as those along the background trajectory in the basis e K 1 .It is also worth noting that the whole dynamics of the background fields is encapsulated only in e K 1 and ėK 1 as E K ≡ E 2 e K 2 vanishes by construction.In this scenario, the rate of Equation 2.6 clearly depends on the entropy perturbations resulting from the presence of the orthogonal field to the homogeneous trajectory: where L ,2 ≡ e K 2 L ,K is the projection of L ,K (i.e., of the derivative of total Lagrangian L with respect to the fields) on e K 2 (i.e., on the direction of the entropic projection of the field acceleration).Clearly, with the presence of the entropy perturbations, the change of ζ could be significant in addition to the general geometry of the field space.
As previously said, we aim to precisely reconstruct the whole inflationary dynamics and compute the cosmological observables.Therefore, we start by requesting the slow-roll conditions In this regime, we can calculate the primordial scalar spectrum predicted by inflation, evaluated at Hubble radius exit, as where c s is the sound speed of the adiabatic perturbation, and is the usual slow-roll parameter.Since during the slow-roll phase the expansion rate H is almost constant, we expect the spectrum of primordial perturbations to be almost scaleinvariant.Therefore one can expand in terms of the small scale-dependence as (2.13) The very same procedure can be repeated for the power spectrum of primordial tensor modes, eventually achieving another useful quantity to constraint general models of inflation, namely the so-called tensor-to-scalar-ratio which quantifies the fraction of primordial gravitational waves produced by the super-adiabatic amplifications of zero-point quantum fluctuations.
Turning our attention to super-Horizon scales, as seen from Equation 2.7, once k 2 /a 2 → 0, we eventually get with A and B model-dependent functions of time, whose expressions will be explicated in subsection 3.1 for a specific model, and a dimensionless gauge-invariant quantity.In order to account for the correlation between curvature and isocurvature modes, as well as to estimate the transfer of entropy from the latter to the former during the time from soon after the Hubble exit to the end of inflation, we make use of the transfer matrix formalism [181] ζ where the transfer functions relate the power spectrum at the end of inflation with the power spectrum at the Hubble exit: with Θ being defined as the transfer angle.

Integration Scheme
After specifying the initial conditions for the fields ϕ K , their velocities, and the metric G IJ (ϕ K ) of the field space, we can integrate the equations of motion, Equation 2.2.The integration process is carried out for a maximum number of e-folds, which is set to N max = 10000.
During the integration, we dynamically calculate the slow-roll parameter ϵ by Equation 2.10 and continue until the condition ϵ = 1 is satisfied.If this condition is not met within N max efolds, the model is rejected as inflation does not end.On the other hand, if the condition is satisfied during integration, the point ϵ = 1 in the parameter trajectory is considered as a possible end of inflation.
To confirm that it represents the actual end of inflation, we check whether the fields are still active enough to begin a second stage of expansion.Specifically, we test whether the normalized field values with respect to their initial conditions do not exceed a threshold value. 3f the field values satisfy the condition, the multifield dynamics can be considered effectively complete, and the point ϵ = 1 reached during integration is regarded as the actual end of inflation.Instead, if the field values do not meet this condition, the integration continues until they fall below the specified threshold.During this stage, we monitor the value of ϵ to test whether inflation restarts (i.e., if we get back to ϵ < 1).If inflation does not restart, then the original point ϵ = 1 is set as the end of inflation.If inflation does restart, the new end of inflation is determined by the joint conditions of ϵ = 1 and the field values, which must satisfy the selected threshold 4 .
After identifying the end of inflation, we proceed to calculate the total interval of efold ∆N between the start of integration and the end of inflation.We make sure that this interval is greater than a threshold value, which we set at ∆N ≥ 100.This threshold is crucial to ensure that inflation lasts long enough to account for the observed homogeneity and isotropy of the Universe, and to establish the appropriate initial conditions for the subsequent Hot Big Bang Theory evolution.When the total number of e-folds between the start and end of integration is less than this threshold value, we perform a diagnostic test aimed at determining whether the smaller number of e-folds is due to the initial conditions being chosen too close to the end of inflation or if the model, for that specific combination of parameters, is unable to sustain slow-roll dynamics for an adequate duration.To conduct this test, starting from the same original initial point of integration with the very same initial conditions, we integrate backwards in time from the remaining missing e-folds.During this additional integration, we check the inflationary fields and parameters for any problematic behavior, such as exponentially divergent trajectories, and verify that the model can indeed smoothly support the desired number of e-folds of expansion.If the model fails to meet these requirements, it is deemed unsuitable for explaining the observations and rejected.Note that this test is primarily included to maintain a conservative approach and save a limited number of configurations where models do not respect the threshold value due to the slight shift caused by the random choice of initial conditions while not showing pathological behavior of the backwards-in-time trajectory.This also reduces the impact on the results from the initial conditions themselves, as argued in section 3. 5After ensuring that the model supports a satisfactory number of e-folds of expansion, and having carefully reconstructed the field dynamics during the entire inflationary phase, we can accurately calculate the entire evolutionary history.This includes how the slow-roll parameters and observables evolve as a function of N .By doing so, we can obtain the value of all the slow-roll parameters at horizon crossing, which we choose to be N ⋆ = 55 e-folds before the end of inflation.Additionally, we take into account the evolution on super-Horizon scales and the transfer of entropy between isocurvature and scalar perturbations.In this way, we can accurately calculate the spectrum of primordial scalar modes (and all the relative observables) at the end of inflation, by means of the transfer matrix formalism detailed in subsection 2.1.

Sampling Method
Once the integration process successfully ends, we can access all the observable predictions of the model, including the amplitude of the scalar perturbation spectrum (A s ), its spectral index (n s ), the scalar running (α s ), the scalar running of running (β s ), and the amplitude of the tensor perturbations (r).To calculate the subsequent cosmology, we interface our algorithm with standard Boltzmann integrator codes.Specifically, we use the "Code for Anisotropies in the Microwave Background" CAMB [174,175]. 6We input the observable predictions of the multifield inflation as initial conditions for CAMB, along with any other relevant cosmological parameters, such as the other standard ΛCDM parameters: Ω b h 2 , Ω c h 2 , θ MC , and τ .This allows us to relate the predictions of the multifield model to the usual observable quantities such as the angular power spectra of cosmic microwave background anisotropies and polarization and the matter power spectrum, in both standard and non-standard cosmological backgrounds.
The next crucial step is to explore the parameter space of generic multifield models by sampling over the initial conditions and the free parameters using Monte Carlo techniques.This enables us to compare the theoretical predictions with observational data and derive observational constraints.To accomplish this goal, we have developed a novel sampling algorithm able to explore a sufficiently large volume of the parameter space and identify a sub-region where the model's predictions agree with observations.Our sampling algorithm works as follows.We input a large number of Monte Carlo steps (∼ 10 6 steps).At each step, we randomly select the values of the initial conditions and model parameters within some prior ranges for all sampled parameters.Subsequently, we integrate the model for these initial values, performing all the consistency tests detailed in the previous subsection and keeping only models in which the end of inflation is clearly identified and inflation lasts for a sufficiently long number of e-folds able to explain homogeneity and isotropy.If the model satisfies these initial prerequisites, we calculate the observable predictions for all the inflationary parameters, such as the predicted values of A s , n s , α s , β s , and r.We then test whether these values fall in reasonable ranges, retaining only models that simultaneously satisfies the following conditions: where we have conservatively chosen the previous ranges around the values measured by the most recent cosmic microwave background experiments. 7If the model falls within these ranges, we calculate its full cosmology by means of CAMB.In this case, we save as output all the relevant information, including the initial conditions, the values parameters, and all observables.If instead the model fails to satisfy any of these conditions, it is rejected.
By following this procedure, we generate a chain of models that are equivalent to those produced by typical Markov Chains Monte Carlo (MCMC) techniques.During the sampling process, each saved model is assigned a likelihood based on how well it agrees with the most recent observations of the Cosmic Microwave Background (CMB).Specifically, our reference datasets include: • The Planck 2018 temperature and polarization (TT TE EE) data, which also includes low multipole data (ℓ < 30) [193][194][195].
• The Planck 2018 lensing data, derived from measurements of the power spectrum of the lensing potential [196].
• The latest CMB B-modes power spectrum likelihood cleaned from the foreground contamination as released by Bicep/Keck Array X Collaboration [9].
To extract a likelihood for each model, starting from these observations we develop an analytical likelihood based on a multi-dimensional normal distribution: where µ and Σ represent the mean values and covariance matrix of parameters obtained by a joint analysis of the aforementioned experiments.We validate our methodology by verifying that both our analytical likelihood and sampler produce consistent results compared to those obtained by using the original experiment likelihoods and publicly available samplers.We refer to Appendix B for further details.Consequently, we can obtain informative posterior distributions for the most relevant parameters to be inferred from observations.

A Case Study
In this section, we provide a working example of the potentiality of our method.As a case study, we analyze a specific model detailed in subsection 3.1.In subsection 3.2 we demonstrate how our algorithm enables us to track the entire fields dynamics and how we can access all the corresponding observables.Finally, in subsection 3.3, we explore the parameter space of the model by sampling over its free parameters and the initial conditions of the fields, deriving observational constraints.

Model
We focus on the simple case where two scalar fields ϕ K = (ψ, χ) are coupled through the field metric and the action reads as where I, J ∈ {1, 2} and g is the determinant of the metric g µν .These models have been recently studied in Ref. [124] where it was argued that a non-trivial field metric could induce significant changes in the effective mass of the entropy perturbations, opening up a rich phenomenology that we can test for further validations.For the aim of this section, it is worth noting that the fields dynamics is dictated by Equation 2.2 which eventually simplifies to [124] By performing a rotation in field space, the evolution of the vector along the homogeneous trajectory (the so-called adiabatic field σ) is given by where σ = √ 2L kin .Instead, the entropy part of the equations of motion is given by the rate of change of the angle between the initial field basis (ϕ, χ) and the adiabatic/entropy one (σ, s): To consistently discuss the perturbations of a given cosmological configuration, we take into account the perturbations of the metric Φ = Φ(t, x) in the longitudinal gauge [197] and the corresponding fluctuations of the sources such as ψ = ψ(t) + δψ k (t, x) and χ = χ(t) + δχ k (t, x).By considering Equation 3.7 and the perturbed Einstein equations, we write the comoving curvature perturbation Equation 2.6 in terms of the metric fluctuations as [124] ζ and its evolution as where S = H σ δs is the so-called isocurvature perturbation, a gauge-invariant quantity.Note that from Equation 2.16 it follows that E 2 = δs.Clearly, the quantity labelled as adiabatic perturbation Equation 3.9 on super-Hubble scales (i.e.k 2 /a 2 → 0) is solely fed by the entropy perturbation δs, namely by the orthogonal field to the background trajectory in field space.Indeed, even if the perturbations in the entropy field evolve independently from the perturbation in the adiabatic field, the large-scale entropy perturbations do impact the evolution of the adiabatic one when the value of the potential curvature is non-zero, i.e. η ,σs ̸ = 0. Furthermore, their coupling (encoded in the term V ,s ) does not vanish even when a flat field metric is considered.In other words, θ ≡ 0. For this model, the time-dependent dimensionless functions A and B in the transfer function, Equation 2.17 are [124] where η ,M N = V ,M N /3H 2 describes the curvature of the potential in terms of the entropy and adiabatic fields.These functions encode the coupling between adiabatic and entropy modes and enter in the expressions for the spectral index and its runnings at the end of inflation that, in this general two-field model, can be derived from Equation 2.12, Equation 2.13, and Equation 2.18 reading [198] n s ≃ n * − 2 sin Θ(A * cos Θ + B * sin Θ), From the expressions above we see that constraints on the spectral parameters can be translated into constraints on the geometry of the field metric, namely the curved trajectory in field space between Hubble radius exit and the end of inflation.

Predictions
To be concrete and maintain control over the results, we test all the features of our method by analyzing a simple case where the coupling between the two fields is given by with a self-coupling potential of the form where m ψ and m χ are the mass terms of the fields and g the coupling constant.Once the field metric and the self-coupling potential are specified, the formalism developed in the previous section can be applied and the multifield dynamics numerically solved by means of the integration scheme discussed in subsection 2.2.In particular, the integration of the equations of motion, and therefore the trajectory of the fields, will depend on the model's free parameters (i.e., m ψ , m χ , g and c) and the initial conditions of the fields (i.e., ψ ini , ψini , χ ini and χini ).In this subsection, we analyze separately the contribution of these parameters to both cosmological observables and the inflationary dynamics, in order to have a comprehensive understanding of their effect and to better interpret the results of the full Monte Carlo analysis performed in the subsequent subsection.
We start by examining the robustness of the outcomes of the model when subjected to variations in the initial conditions of the fields.Specifically, we evaluate how the trajectories in the field space change by randomly varying the starting points ψ ini and χ ini , while keeping the model's free parameters fixed and setting ψini = χini ≃ 0. The results are depicted in Figure 1.Considering different initial conditions for the fields can result in a period of inflation that is more or less long, while essentially keeping the model's predictions for cosmological observables unchanged.Therefore, we do not expect physical observables to drastically depend on the initial conditions of the fields.
On the other hand, the variation of the free parameters holds greater significance.For instance, the effects of c which are encoded in the curved field manifold play a substantial role in what concerns the interplay between isocurvature and curvature modes between the horizon crossing and the end of inflation.In particular, when c is gradually negatively reduced, the χ field velocity remains relatively constant and it becomes stuck for a significant amount of time until it reintegrates into the dynamics when the ψ field reaches the minimum of the potential, see also Figure 2.This effect is translated into the amount of isocurvature mode feeding the curvature one.On the other hand, if c is progressively incremented will lead to a single field scenario as the χ field drops immediately. 8In other words, the combination of curvature and isocurvature perturbations may lack correlation.However, if the path in field space exhibits curvature from the Hubble exit until the end of inflation, it can introduce an indefinite level of correlation between them impacting the ultimate predictions and potentially causing their growth or reduction, see Figure 3. Shifting our focus to the effects of the masses, from Figure 3 we can appreciate how changes in m χ and m ψ lead to increase or decrease of power in the spectrum of temperatures anisotropies.This can be explained by making explicit Equation 2.8 as and taking into account a self-coupling potential of the form Equation 3.13, one can easily get by means of dH dχ = Ḣ χ .Finally, we achieve the following relation Therefore, on the slow-roll trajectory defined by Equation 3.15 there is a correlation between the mass m χ and the coupling b ,χ .Such a correlation does not exist for m ψ , for which we obtain stating that an increase in m ψ leads to a rapid change in the Hubble rate ( σ2 ≃ − Ḣ) which translates into a shorter period of inflation and hence a reduction of the amplitude of the power spectrum as seen in Figure 3.The same figure shows that increasing either c and m χ results in a larger amplitude of the initial power spectrum, as expected from Equation 3.16.

Monte Carlo analysis and parameter constraints
We use the sampling algorithm detailed in subsection 2.3 to explore the parameter space of our multifield model.The sampling involves 4 free parameters (m ψ , m χ , g, and c) and the initial conditions of the fields (ψ ini , χ ini ), which determine the field trajectory during inflation, as pointed out in subsection 3.2.To ensure that we are able to explore a sufficiently large volume of the parameter space, we randomly vary the model parameters and initial conditions in the large uniform priors listed in parameters and initial conditions as large as the number of steps in the Monte Carlo analysis.
For each step, we use the integration algorithm described in subsection 2.2 and compute the full evolution of the fields as well as any observable quantities of the model, including all the primordial scalar and tensor spectrum parameters such as A s , n s , α s , β s and r.We test the agreement of each combination of parameter and initial conditions against data by using our likelihood Equation 2.19 built on the Planck-2018 [193][194][195][196] and BK18 [9] observations of the cosmic microwave background temperature and polarization anisotropies.
Following this methodology, we perform a full Monte Carlo analysis, repeating the process for ≳ 5 × 10 6 steps and collecting about 2 × 10 4 sampled models, each of one with its own likelihood.In this way, we derive informative posterior distributions for all the parameters involved in the sampling, as well as for any derived quantities of interest, including those associated with entropy transfer on super-horizon scales.The results are summarized in Table 1, while Figure 4 depicts the distribution of the sampled models in the 4D parameter space, represented by a box with dimensions corresponding to the size of the prior volume.
Figure 5 shows instead the 68% and 95% CL contour plots for all the quantities of interest in the model.
The first important test we can perform is to study the constraints on the inflationary observables, such as the amplitude of the primordial scalar spectrum A s , the spectral index n s , its runnings α s and β s , and the amplitude of primordial tensor modes r.Notice that when dealing with multifield inflation, fixing a model does not necessarily establish consistency relations between these parameters as it usually happens in single-field inflation.Consequently, the specific values of these observables depend on the interplay between the free parameters of the model and the initial conditions of the fields.Additionally, it is important to account for corrections that arise due to the super-horizon evolution of adiabatic and isocurvature perturbations, as discussed in previous sections.Regarding the amplitude of the scalar spectrum, we obtain A s = ( 2.109 ± 0.033 ) • 10 −9 at 68% CL, in perfect agreement with the model-independent analysis performed with the full Planck and BK18 likelihoods (see also Appendix B).Similarly, for the spectral index we get n s = 0.9621 +0.0053 −0.0047 at 68% CL, while the for the amplitude of primordial gravitational waves we obtain an upped bound r < 0.04 at 95% CL.All these results are consistent with the most recent joint analyses of Planck and BK18 data, as well.Regarding the higher-order runnings of the primordial scalar spectrum, we obtain at 68 % CL α s = −0.74+0.37  −0.32 × 10 −3 and β s = −0.103+0.088 −0.004 × 10 −3 , respectively.Therefore, assuming this particular multifield model, smaller negative values are favored for both the running and the running of the running, although both of them remain consistent with null values at 95% CL.
One significant aspect of our approach is the ability to obtain observational constraints on the model parameters.For instance, within this particular model, we are able to place a 95% CL upper bound on the mass values of the fields that (in Planck units) read m ψ < 2.30 • 10 −6 and m χ < 1.01 • 10 −5 , respectively.Similarly, for the coupling parameter g we obtain g < 9.72 • 10 −7 at 95% CL, while for the parameter c controlling the curvature of the filed space we obtain c < −0.0211 always at 95% CL.It is worth noting that these upper bounds assume values significantly far from the limits of the priors adopted for these parameters.In this regard, we have verified that the priors chosen for parameter exploration are sufficiently large to provide uninformative ranges, without introducing any unwanted bias in the parameter constraints.In order to study the correlation between the different parameters, we can refer to Figure 3 for the effects on the angular spectra, Figure 4 for the correlation in the 4-D parameter space of the model, and Figure 5 for the contour plot with all sampled and derived parameters.From Figure 3, one would expect a correlation between the effects of the parameters g and c on the spectrum of CMB temperature anisotropies as considering more negative values of c leads to a power amplification, while increasing g reduces the amplitude of the spectrum.Therefore, we expect that more negative values of c can be allowed only for larger values of g, and this is clearly confirmed by both Figure 4 and Figure 5. Notice also that this model strongly prefers highly negative values of c.The reason behind this preference is that when c becomes positive, the curvature of the field space G IJ (ϕ K ) is exponentially suppressed, and one of the two fields essentially becomes a spectator field.As a result, the model reduces to a single-field model with a quadratic potential, which is well-known to predict a significantly large amount of primordial gravitational waves, higher than the current observational limit.Therefore, the parameter space with positive values of c is severely constrained by the B-mode polarization measurements by BK18, in combination with the Planck temperature and polarization measurements.
Finally, our algorithm allows us to derive constraints on any relevant physical quantities in the model, including parameters and functions that govern the transfer of entropy between adiabatic and isocurvature perturbations.For example, we derive a 95% CL upper bound on the angle Θ < −0.686 appearing in Equation 3.11 that weights the corrections acquired by the inflationary parameters between the Hubble exit and the end of inflation.Similarly, we can constrain the time-dependent functions A(t) and B(t) involved in the transfer matrix formalism used to account for the correlation between curvature and isocurvature modes, as well as to estimate the transfer of entropy from the latter to the former during the time from soon after the Hubble exit to the end of inflation.We evaluate these functions at the Hubble exit, finding at 95%CL the following limits A ⋆ > −1.71 and B ⋆ > −0.341.Our results confirm that in multi-field models with non-flat metric spaces, the interplay between isocurvature and adiabatic modes plays a crucial role in the final observable predictions, as already argued in Ref. [124].

Conclusion
Embedding inflation within a more fundamental framework still remains an open problem, and numerous models and mechanisms have been proposed.The simplest dynamical models of inflation involve a single scalar field minimally coupled to gravity whose evolution is governed by a potential that should be enough flat to induce a phase of slow-roll evolution.However, in certain low-energy effective field theories inspired by string theory, the scalar field sector often comprises multiple scalar fields with non-standard kinetic terms or dynamical couplings.When inflation is driven by multiple scalar fields with non-standard kinetic terms, the interplay between adiabatic perturbations and isocurvature modes becomes significant, influencing the observable predictions and giving rise to a rich phenomenology that can be tested using current cosmological and astrophysical data.In principle, precise measurements of the cosmic microwave background radiation can be used, providing stringent constraints on the abundance of both adiabatic and isocurvature modes and offering a valuable opportunity for experimental validation of these models/theories.However, in practice, obtaining precise predictions from multifield inflation is not always easy as observational quantities often depend on various factors such as the initial conditions of the fields and the specific model assumed.
Different trajectories in field space can lead to different results, changing the amplitude of scalar and tensor primordial perturbations.For this reason, most tools employed in cosmological data analyses, such as typical Boltzmann integrator codes and samplers, are either unaware of the physics of inflation or assume single-field potentials.As a result, constraining the multifield landscape with current CMB data represents an ongoing challenge.
In this work, we take a first step to bridge this gap and introduce an algorithm specifically designed to investigate generic multifield models of inflation where a number of scalar fields ϕ K are minimally coupled to gravity and live in a field space with a non-trivial metric G IJ (ϕ K ).In section 2, we describe both our theoretical parameterization and our numerical method.In particular, after specifying the initial conditions for the fields ϕ K , their velocities, the metric G IJ (ϕ K ) of the field space, and the self-coupling potential, our algorithm is able to reconstruct the dynamics of the fields throughout the entire inflationary period precisely determining the end of inflation and performing several consistency tests to ensure the stability of the model.This comprehensive analysis includes more intricate scenarios where distinct field dynamics govern various stages of expansion at different times, such as double inflation or punctuated inflation.An illustrative example can be found in Appendix A. By numerically solving the full field dynamics , we can calculate precise predictions for observable quantities in the slow-roll limit, such as the spectrum of scalar perturbations, primordial gravitational waves, and isocurvature modes.We can also track the super-horizon dynamics of adiabatic and isocurvature modes, determining the transfer of entropy to scalar modes after horizon crossing using the transfer matrix formalism described in subsection 2.1.Once the integration process successfully concludes, we can access all the observable predictions of the model and set the initial conditions to compute the subsequent cosmology.This is done by interfacing our algorithm with standard Boltzmann integrator codes such as CAMB or CLASS that allow us to directly translate the model's predictions in terms of the cosmic microwave background anisotropies and polarization angular power spectra.
Based on this algorithm, we also introduce a novel sampler which is specifically designed to explore a sufficiently large volume of the parameter space of a generic multifield models and identify a sub-region where the model's predictions are in agreement with observations.This allows us to efficiently sample over the initial conditions of the fields and the free parameters of the model, enabling Monte Carlo analysis to compare theoretical predictions with observations.In this work, we make use of the most recent CMB data provided by the Planck collaboration, as well as the latest B-modes power spectrum likelihood released by the Bicep/Keck Array X Collaboration and extract a likelihood for each sampled combinations of parameters.To do so, we develop an analytical likelihood based on these observations that has been extensively tested and proven to reproduce the results obtained from the real likelihoods of different experiments for inflationary parameters, see Appendix B.
In section 3, we provide a detailed illustration of our approach by analyzing a specific case study model where two scalar fields are coupled through the field metric by Equation 3.12, with a self-coupling potential given by Equation 3.13.Focusing on this model, in subsection 3.2 we use our integration scheme to investigate how the observable predictions change with its parameters and initial conditions.We refer to Figure 1 and Figure 2 for the impact of the field trajectories and the coupling function, respectively.Instead Figure 3 illustrates the impact of the different parameters on the CMB angular power spectra of temperature anisotropies.Finally, in subsection 3.3, we employ our sampler to perform a comprehensive Monte Carlo analysis, deriving observational constraints on the free parameters.The results of our analysis are summarized in Table 1, while Figure 4 and Figure 5 show the distribution of sampled models in the parameter space and the correlation among the different parameters, respectively.We are able to derive compelling and precise constraints on both the model's parameters and the inflationary observables such as the primordial power spectra of scalar and tensor perturbations.In addition, we can place constraints on the interplay between curvature and isocurvature modes by accurately accounting for the entropy transfer from isocurvature to curvature perturbations on superhorizon scales.For sake of completeness, in Appendix A, we also discuss the limit where the model reduces to the case of a double quadratic potential with a canonical kinetic term, discussing our ability to recover results already documented in the literature.
Our work provides a robust framework for exploring multifield inflation and opens up exciting opportunities for future research focused on the rich phenomenology expected in both standard and non-standard multifield models of inflation and gravity.

A Double Quadratic Potential & Double Inflation
In this appendix, we make use of our code to briefly examine a simplified two-field model compared to the one analyzed in section 3. Specifically, we investigate the case of a doublefield quadratic potential with a canonical kinetic term which falls within the parameterization adopted in section 3 once we fix c = 0 and g = 0.
Notice that this model has been already discussed in the literature and a comprehensive systematic review of its properties is beyond the goal of our work.That being said, this appendix serves a dual purpose: on one side, we aim to use this simplified case as a safety check to demonstrate our ability to recover many of the results already documented in the literature.On the other hand, we take this opportunity to provide a working example of some features of our algorithm, such as its ability to correctly identify the end of inflation and effectively handle scenarios of double inflation where two stages of expansions are driven by distinct fields at distinct times 9 .
Regarding the topic of interest for our discussion, an analysis similar to the one we carry on here is given in Ref. [203] where, using a parameterization for the inflationary dynamics very close to the one we have developed in this paper, the authors point out the characteristics of this two-field model, deriving predictions for observables such as the amplitude of primordial spectra.As highlighted in Section III.C of Ref. [203], a distinctive feature of this potential is the prediction of a tensor-to-scalar-ratio r ≳ 0.13, regardless of the initial conditions for the fields and the values of the mass ratio.Given that a similar value for the tensor amplitude does not appear to be in agreement with the most recent measurements of CMB B-mode polarization released by the Keck Array collaboration, we do not perform a data analysis (which would be inconclusive given the model's inability to reconcile with the latest observations).However, we perform a model sampling as a cross-validation of the predictions for r.As seen in Figure 6, collecting approximately 10 4 models where the amplitude of the models where the amplitude of the scalar spectrum A s and its tilt n s fall within ranges consistent with our observational data.scalar spectrum A s and its tilt n s fall within ranges consistent with our observational data, we find that the the tensor amplitude r remains consistently above the 95% limit resulting from the joint Planck+BK18 analysis.Furthermore, the value of r predicted by this potential remains largely unchanged both with respect to the mass ratio (across multiple orders of magnitude) and with respect to the initial conditions (which are randomly chosen from model to model) as well as in excellent agreement with Figure 8 of Ref [203].
A closer analysis of the figure reveals the presence of a minor dispersion of patterns that deviate from the degeneration line between the mass ratio and the value of r.A detailed analysis of these models has revealed that this small dispersion is due to double inflation.While theoretically investigated, in Ref. [203] these models were not given thorough consideration and the inflationary dynamics were computed by treating double inflation as if it consists of only one inflationary phase.However, as explained in subsection 2.2, our algorithm has been designed to accurately identify all inflationary phases and comprehensively reconstruct their dynamics.As a result, we can provide precise predictions for scenarios involving double inflation, as well.
To put it more quantitatively and clarify our treatment of double inflation, we focus on one specific model among those depicted in Figure 6.Namely, we consider the one which shows the smallest tensor amplitude (r ∼ 0.05, yet outside the Planck-BK18 range) and provide full details of the dynamics of the field and the background evolution in Figure 7.As evident from the top panel of the figure, in this model, the parameter ϵ reaches the critical value ϵ = 1 for the first time after about 55 e-folds (grey dashed line).During this initial phase, the expansion of the Universe is driven by the field ψ, whose evolution is depicted in the second panel of the same figure.As one can observe, when ϵ = 1, the dynamics of ψ is mostly completed, while the second field χ (whose evolution is shown in the third panel) remains mostly inactive during this phase.Once the inflationary stage guided by ψ terminates, our algorithm monitors the dynamics of the fields and, since χ is still very active, the integration of the equation of motion continues until the second field also completes its dynamics.This allows us to identify a second inflationary phase, this time guided by χ.In both of these Cobaya [204] in combination with the full Planck 2018 [193][194][195][196] and BK18 [9] likelihoods (referred to as 'Real likelihoods'), and the results for the same model derived using our sampling algorithm in combination with our analytical likelihood (Equation 2.19) (referred to as 'This work').
stages, the background dynamics, represented by the evolution of the Hubble parameter H (bottom panel), is accurately computed together with all observational predictions, including the tensor amplitude.

B Sampling and Likelihood Validation
In this appendix, we provide a step by step explanation of the methodology used to build our likelihood based on the joint analysis of B-Modes polarization data from BK18 and the Planck-2018 measurements of temperature and polarization anisotropies.Most importantly, we prove that our method/likelihood is able to reproduce the same results obtained by the most commonly Markov Chain Monte Carlo (MCMC) analyses performed in the literature.
To do this, through all this appendix, we do not consider any explicit inflationary models, and remain completely agnostic about the physics of inflation, as customary in the literature.
1) As a first step, we consider an extension to the standard cosmological model which includes three additional parameters: the tensor amplitude r, the running of the spectral index α s , and the running of the running β s .We refer to this model as ΛCDM+α s +β s +r and perform a full Monte Carlo Markov Chain (MCMC) analysis using the publicly available sampler Cobaya [204], and Boltzmann integrator code CAMB [174,175].Cobaya explores the posterior distributions of the parameter space using the MCMC sampler developed for CosmoMC [205], which has been specifically adapted for parameter spaces with a hierarchy of speeds by implementing the "fast dragging" procedure introduced in Ref. [206].The baseline datasets involved in our MCMC analysis consist of the Planck 2018 temperature and polarization (TT TE EE) likelihood, which includes low multipole data (ℓ < 30) [193][194][195], as well as the lensing likelihood obtained from measurements of the power spectrum of the lensing potential [196].Additionally, we include the most recent CMB B-modes power spectrum likelihood cleaned from foreground contamination, as released by the Bicep/Keck Array X Collaboration [9].The resulting constraints are presented in Table 2, and we also take this opportunity to update the current bounds on this cosmological extension with the latest data.[193][194][195][196] and BK18 [9] likelihoods (grey), and our sampling algorithm in combination with our analytical likelihood (red).
2) As a second step, using the results from our MCMC analyses, we construct our analytical likelihood based on Equation 2.19.In particular, we adopt the generalized covariance matrix Σ and the mean values µ of parameters obtained for the ΛCDM+α s +β s +r model.
3) Finally, we test that our likelihood is able to reproduce the same constraints as real likelihoods from experiments.This is a crucial step needed to validate the results obtained when our likelihood is applied to the study of physical models of multifield inflation.To prove the equivalence of our method and the MCMC analysis, we perform a consistency check for the same ΛCDM+α s +β s +r cosmological model.Specifically, we perform a new run varying the inflationary parameters {A s , n s , α s , β s } using our sampling algorithm.In this case, the sampling is performed directly in the parameter space of inflationary observables and using the same prior ranges used for the MCMC analyses.Since sampling on the inflationary parameters requires significantly fewer computational resources, avoiding also some non-physical behaviors that may arise when exploring non-standard inflationary models or solving the inflationary dynamics, we are able to accumulate ∼ 10 5 models, each of them evaluated using our analytical likelihood Equation 2.19.In Figure 8, we directly compare the results obtained using our approach (labeled as 'this work') to the results obtained using the real likelihoods from experiments and the sampler Cobaya (labeled as 'Real likelihoods').The two methods

Figure 1 :
Figure 1: Field trajectories (as well as their projections in different 2D-planes in grey) as functions of the e-folds of expansions between the beginning of the integration and the end of inflation.The integration process begins at N = 0 with randomly selected initial conditions, represented by a black star-like dot in the figure.The end of inflation (marked by another black star-like dot) is determined using the method explained in subsection 2.2.The color-bar shows the value of the Hubble parameter along the field trajectories.For all trajectories, the model's free parameters are fixed to: m ψ = 1.58 × 10 −6 , m χ = 3.86 × 10 −6 , c = −0.06,and g = 2 × 10 −8 .

Figure 2 :
Figure 2: Effects of the field metric parameter c on the field evolution in the ψ vs χ plane.The field values have been normalized to their initial conditions, so the starting point of the trajectories is represented by the coordinates (1, 1), while the coordinates (0, 0) ideally represent the end of inflation.The red (blue) curves correspond to negative (positive) values of c as indicated in the legends, while the green curve corresponds to a flat filed metric G IJ = diag{1, 1}, (i.e., c = 0).The other model's free parameters are fixed to: m ψ = 1.58 × 10 −6 , m χ = 3.86 × 10 −6 , and g = 2 × 10 −8 .

Figure 3 :
Figure 3: Effects on the CMB angular power spectrum resulting from variations in the model's parameters as indicated in the different panels/legends of the figure.The baseline model (black line) corresponds to the following parameter combination: m ψ = 1.58×10 −6 , m χ = 3.86×10 −6 , c = −0.06,and g = 2 × 10 −8 .

Figure 4 :
Figure 4: How the models distribute in the 4-dimensional parameter space.The box has the size of the prior volume

2 Figure 5 :
Figure 5: Marginalized 2D and 1D posteriors distributions for all the model's parameters and the quantities of interests in this study.

Figure 6 :
Figure 6: Predictions for the tensor amplitude r in terms of the mass ratio log 10 (m χ /m ψ ) for ∼ 10 4

Figure 7 :
Figure 7: Evolution of the parameter ϵ (top panel), the scalar fields ψ and χ (second and third panel, respectively), and the Hubble parameter H (bottom panel) in a model of double inflation characterized by two stages of expansion.
≡ P ζ (k * ) is the amplitude of the scalar spectrum at the pivot scale k ⋆ , which we fix to k ⋆ = 0.05 Mpc −1 throughout this work.In the simplest parameterization any residual scale dependence of the spectrum is parameterized solely in terms of the scalar spectral index n s defined as:

Table 1 .
This enables us to obtain a number of combinations of

Table 2 :
Results for the ΛCDM+α s +β s +r model obtained using the publicly available sampler