Searching for local features in primordial power spectrum using genetic algorithms

We present a novel methodology for exploring local features directly in the primordial power spectrum using a genetic algorithm (GA) pipeline coupled with a Boltzmann solver and Cosmic Microwave Background data (CMB). After testing the robustness of our pipeline using mock data, we apply it to the latest CMB data, including Planck 2018 and CamSpec PR4. Our model-independent approach provides an analytical reconstruction of the power spectra that best fits the data, with the unsupervised machine learning algorithm exploring a functional space built off simple ``grammar'' functions. We find significant improvements upon the simple power-law behaviour, by $\Delta \chi^2 \lesssim -21$, consistently with more traditional model-based approaches. These best-fits always address both the low$\ell$ anomaly in the TT spectrum and the residual high$\ell$ oscillations in the TT, TE and EE spectra. The proposed pipeline provides an adaptable tool for exploring features in the primordial power spectrum in a model-independent way, providing valuable hints to theorists for constructing viable inflationary models that are consistent with the current and upcoming CMB surveys.


INTRODUCTION
The most cogent case for inflation comes from its ability to make sense of the homogeneity at very large angular scales in the Cosmic Microwave Background.These otherwise puzzling measurements are readily explained by an accelerated expansion that took place very early in the Universe's evolution.Additionally, cosmic inflation provides a natural mechanism for quantum fluctuations, which give rise to density fluctuations that eventually grow into the large-scale structure we see today.
The simplest inflationary paradigm delivering a viable cosmology consists of a scalar field slowly rolling down an almost flat potential, resulting in adiabatic scalar perturbations with a nearly scaleinvariant power spectrum.The conventional way to characterise this scalar power spectrum is via a power law of the form: where  s ,  s and  * = 0.05 Mpc −1 are respectively the amplitude, spectral index, and pivot scale of the primordial spectrum.Similarly, one may show that, in the simplest scenario, tensor fluctuations are produced from vacuum fluctuations and engender a primordial tensor power spectrum that is also well described by a power law   () =  t (/ * )  t .Planck's 2018 data analysis has accurately measured the scalar power law spectrum to be slightly ★ E-mail: shafieloo@kasi.re.kr red-tilted and placed constraints on the maximum amount of tensor modes at CMB scales.Assuming the power law, such measurements reveal that the inflationary models must correctly predict ln(10 10  s ) = 3.044±0.014, s = 0.9649±0.0042at 68% confidence level, and  0.002 ≃  t / s < 0.10 at 95% confidence level (Planck Collaboration 2020a), with even tighter constraints on primordial tensor modes when combined with CMB data from the BICEP-Keck collaboration:  0.05 < 0.035 (Ade et al. 2022).
The motivation for exploring features beyond this simple power law in the primordial scalar power spectrum is manifold.Firstly, they provide valuable insights into the dynamics of inflation itself.Indeed, several compelling models of inflation comprise transient deviations from the slow-roll dynamics.These features in the inflationary trajectory may be due to localised bumps, dips, saddles points, etc., in the scalar potential or its derivatives or to a sudden turn in multidimensional field space, resulting in some features and oscillations in the scalar power spectrum (see, e.g.Starobinsky (1992); Chung et al. (2000); Adams et al. (2001); Gong (2005); Joy et al. (2009); Achúcarro et al. (2011); Achucarro et al. (2014); Mizuno et al. (2014); Braglia et al. (2020); Fumagalli et al. (2021); Dalianis et al. (2021)).This family of features is called sharp features, for the deviation from smoother slow-roll dynamics is only transient.Another family of features, so-called resonant features, may be produced by tiny oscillations overimprinted on the slow-roll potential, or by the oscillations of an excited massive scalar field in a multifield context (Chen 2012;Battefeld et al. 2013;Gao et al. 2013;Chen & Namjoo 2014;Meerburg et al. 2014;Palma 2015).Contrary to sharp features, resonant features result in oscillations in log , which may be local or global.Other kinds of features may be encountered, leading to peaks and troughs in the scalar power spectrum, broken power laws or plateaus, see (Slosar et al. 2019) for a recent review on primordial features.Thus, identifying and characterising features in the primordial power spectrum is then a key step towards gaining a deeper understanding of the inflationary mechanism, its microphysics and particle content.
Our analysis is further motivated by the fact that, despite the success of the simple power law power spectrum, there are several issues with the standard ΛCDM model, which include the  L anomaly (Planck Collaboration 2020b), the Hubble (Riess et al. 2019) and  8 tensions (see Di Valentino et al. (2021) for reviews).The phenomenological lensing consistency parameter,  L , was introduced to test the consistency of the amount of lensing-induced smoothing in the angular power spectrum data with the theoretical prediction of  L = 1 (Calabrese et al. 2008).Analyses of Planck data have a preference for  L > 1 at the 3 level when considering only power spectrum data, which decreases to 2 when including lensing reconstructions (Planck Collaboration 2020c).The presence of  L > 1 "anomaly" suggests additional smoothing of spectra compared to the predictions of the standard cosmological model.The "Hubble tension" is a 3.4 inconsistency between the measurement of the current expansion rate of the Universe in our close-by environment and the one extrapolated from precision CMB observations at the time of photon decoupling.The presence of features in the primordial power spectrum could potentially offer insights into the underlying causes of these tensions by inducing shifts in other cosmological parameters and reconciling the existing tensions within the standard ΛCDM paradigm (see , e.g.Hazra et al. (2019Hazra et al. ( , 2022)); Ballardini & Finelli (2022); Antony et al. ( 2023)).
We study the form of the primordial power spectrum using the CMB temperature and polarisation data from Planck which offers the most sensitive avenue to search for these signatures.At the CMB scales, the modes exhibit sufficient linearity, allowing us to utilise forward modelling to isolate the primordial information.However, the original signal is significantly suppressed by Silk damping resulting from photon diffusion, particularly beyond  ≳ O (10 −1 ) Mpc −1 , and tends to be overshadowed by foregrounds (Chluba et al. 2015).Besides the CMB, large-scale structure (LSS) surveys offer another opportunity to probe deviations from a featureless primordial power spectrum.Despite the recent advancements in handling nonlinearities for global (not localised in k) oscillatory features (Beutler et al. 2019;Mergulhão et al. 2023), the modelling uncertainties in the nonlinear regime constrain our ability to comprehensively test arbitrary forms of features using large-scale structure data (Hu et al. 2014;Hunt & Sarkar 2015;Benetti & Alcaniz 2016;Zeng et al. 2019).
In order to search for features in the primordial spectrum, two fundamental approaches may be considered: parametric and nonparametric ones.
Parametric approaches involve extending the power law primordial spectrum with a specific functional template and fitting the parameters of this extended model to observational data.This includes a model with varying amplitude and tilt of density perturbations as a power law of the scale of the perturbation (Cline & Hoi 2006;Planck Collaboration 2020c).Another subclass are features where deviations from a smooth power spectrum are confined to a particular range of scales.These "local features" include broken power-law spectrum, sharp peaks and dips, and resonant oscillations in the spectrum at a specific scale (Chluba et al. 2015;Slosar et al. 2019).Parametric approaches have the advantage of being interpretable and computa-tionally efficient to test with the data.However, the downside is that they may not be flexible enough to capture more complex features in the primordial spectrum.
Genetic algorithms (GAs) are meta-heuristic optimisation techniques that are inspired by the process of natural evolution.Since the introduction of GA in the mid-1970s (Holland 1975), they have been widely adopted in various fields, such as engineering, economics, accelerator physics and cosmology, to search for the optimal solution to a particular problem.GAs have been coupled with Wilson-Devinney (W-D) code to optimise the light curve of eclipsing binaries (Metcalfe 1999) and with radiative transfer codes to fit dust spectra from AGB stars (Baier et al. 2010).In the context of cosmology, GAs were used for reconstructing the expansion history of the Universe in a model-independent fashion (Bogdanos & Nesseris 2009;Nesseris & Garcia-Bellido 2012;Nesseris & García-Bellido 2013;Gangopadhyay et al. 2023), to test the validity of ΛCDM in combination with Om-statistics (Nesseris & Shafieloo 2010), growth rates (Arjona et al. 2022), constraining deviations from general relativity (Alestas et al. 2022) and finding approximate analytical functions for transfer functions (Orjuela-Quintana et al. 2022) and sound horizon at drag epoch (Aizpuru et al. 2021).
In this work, we present a novel approach to directly study features in the primordial power spectrum using a genetic algorithm.Our method is designed to explore this high-dimensional feature space by introducing combinations of local and global features in the primordial power spectrum, utilising a grammar that captures these essential aspects.GAs have been recently applied to the search for primordial features in the scalar potential of a single scalar field driving inflation (Kamerkar et al. 2022;Abel et al. 2023), providing a complementary approach to the one we shall pursue here.To validate our approach, we test it on mock data before applying it to the Planck 2018 likelihood and the latest release from CamSpec PR4.Our results demonstrate the effectiveness of our proposed methodology in identifying essential features in the primordial power spectrum, which can provide valuable insights into our understanding of the early Universe.
The paper is organised as follows: In Section 2, we provide a concise overview of the genetic algorithm and pipeline utilised to search for features that we validate with mocks.Section 3 describes the datasets used for the analysis and presents the results.We finally conclude the paper with a summary and discussion of possible future directions in Section 4.

METHODOLOGY
To compute the angular power spectrum of the CMB, one needs to solve the Einstein-Boltzmann equations, which describe the evolution of photon-baryon perturbations in the early Universe.This hierarchy of ordinary differential equations can be solved numerically using cosmological Boltzmann solvers, which take the primordial power spectrum as an input.The solution of the Boltzmann equation yields transfer functions that describe the evolution of the baryon-photon plasma and the growth of the perturbations.The transfer functions depend on various cosmological parameters, such as the density of dark matter, baryons, and dark energy, as well as the Hubble constant and the amplitude of primordial perturbations.We can calculate the CMB temperature and polarisation anisotropies in terms of the primordial power spectrum using the relation: where () is the primordial scalar power spectrum,   ℓ () is the transfer function where  = ,  correspond to the temperature and E-mode polarisation, respectively, and ℓ is the multipole moment.While the constraints on primordial  modes from experiments such as BICEP-Keck (Ade et al. 2022) can provide valuable information on the tensor counterpart of the primordial scalar power spectrum, the fact that Planck 18 data used in this work does not include  modes and only provides a weak constraint on primordial tensors,  0.002 < 0.10, motivates us to leave reconstruction of tensor modes for future works.

Genetic Algorithms
In our approach, we represent the desired power spectrum as chromosomes, akin to the composite structure found in genetics, while individual features can be seen as analogous to 'genes.'A fitness function analogous to a loss function in machine learning acts as a mapping between chromosomes and the properties of the desired solution.It serves as a means to assess the proximity of a candidate solution to the available data.In our case, the fitness function is defined as the sum of  2 values derived from different likelihoods.Through the utilisation of genetic algorithms, we aim to optimise the selection of features and iteratively refine the solution space to identify the most promising configurations that best agree with the given data.
The basic methodology of a GA can be summarised as follows: • Initialisation: A population of individuals is randomly generated, each representing a potential solution to the optimisation problem.
• Fitness Evaluation: The fitness of each individual in the population is evaluated using a fitness function, which measures how well the individual solves the optimisation problem.
• Selection: The fittest individuals are selected from the population to produce the next generation of individuals.The likelihood of reproduction is based on the fitness scores of the individuals.
Table 1.The priors for the Eq. ( 4) that determine the amplitude, shape, width and position of the feature.The initial population for GA analysis is created by randomly sampling this parameter space.These ranges are explicitly enforced by means of constrained crossover and mutation operators.

Parameter
Prior Range A [-2,4] B [0, len(grammar)] C [log(0.8), log(10)] D [log(0.01), log( 5)] • Crossover: The selected individuals are combined to form a new generation of offspring.This is done by taking bits of genetic information from each parent and combining them to form a new individual.
• Mutation: A small percentage of the offspring are randomly mutated to introduce new genetic diversity into the population.
• Termination: The algorithm stops when a specified number of generations is reached, or the individual fulfils pre-defined fitness criteria.

Pipeline
We first generate four sets of four random numbers (A, B, C, D), each from a pre-defined prior range as presented in Table 1.Each set represents a distinct feature, comprising the amplitude, shape from grammar, scale, and position of the feature.These features are linearly combined to generate an individual using Eq (4).Utilizing the Sympy library (Meurer et al. 2017), we transform each individual into an analytical expression, which is then incorporated into the primordial module of CLASS (Lesgourgues 2011) to calculate the corresponding angular power spectrum ( ℓ s).Subsequently, the computed  ℓ s are passed to 'clik' (Planck Collaboration 2020a) to compute the log-likelihood value.Due to the computational complexity of our approach, we keep the background cosmological and nuisance parameters fixed to the BOBQYA (Powell 2009) best-fit values from the vanilla ΛCDM model.
The deviation to power law spectra is modelled by where  = log   * and  () is the linear combination of features in logarithmic space, in our case, we limit its length to four terms, see Eq. ( 4).The selection of combining four features linearly in our methodology is based on preliminary trial runs.If we include fewer features, the GA's effectiveness would be diminished, see (Kamerkar et al. 2022).Conversely, incorporating an excessive number of features could hinder optimisation.These trials aimed to strike a suitable compromise, providing sufficient flexibility in the P(k) space while minimising the risk of poor convergence.The presence of the modulus in Eq. ( 3) is motivated by our will always to have P(k) positive and avoid unphysical solutions with negative power.Note that the GA looks for features that are functions of the log() instead of  to better parametrise the relative power changes across many scales.
So, the overall function reconstructed by the GA is given by where and   () are functions in the grammar as described in Table 2.These two classes, local and global, are motivated by inflationary scenarios where the deviation from the slow-roll dynamics is nearly instantaneous, or on the contrary, present during the whole evolution.For example, if a slow-roll potential is superimposed by a tiny but persistent, periodic modulation,  slow−roll () →  slow−roll () [1 + sin(/Λ)] with  ≪ 1 and Λ ≪  Pl , the power-law primordial spectrum is modified to feature a global oscillatory component, with oscillations linear in log  space across all scales (Chen et al. 2008), therefore also justifying the use of the variable .
Another possibility is that the oscillatory component is localised around a particular region of the inflationary potential or that the resonant feature comes from the temporary excitation of a heavier degree of freedom, such as background oscillations of a second scalar field in multifield inflation, see, e.g., (Chen et al. 2015) for explicit two-field realisations.Such a localised, so-called resonant feature results in oscillations in log  space again, but this time with a finite extension, explaining our local grammar.We comment that yet another possibility is one of the so-called sharp features, where the transient deviation from a slow roll is very short but very violent, and where no particular oscillation is present in the background, see (Slosar et al. 2019;Martin et al. 2023) for a recent review on primordial features and their relations to inflationary physics.Sharp features result in localised oscillations linear in  space and are therefore not encompassed by our present study; we leave for future work the search for sharp features in the CMB data with GA.The fitness of each individual is then evaluated based on its negative log-likelihood (− log L) value.To form a mating pool, we employ a binary tournament selection process where the fittest individual receives the maximum number of copies.From this pool, we randomly select two individuals for each pairing, generating a total of  pop ×  c offsprings via a crossover, with the crossover individuals being mutated with a probability of  m .We set crossover and mutation probability to 0.8 and 0.4, respectively.The remaining individuals ( pop × (1 −  c )) for the next generation are selected from the fittest individuals of the previous generations, ensuring their continued contribution to the evolving population.
For the crossover operator, we adopt a combination of one-point crossover and a bounded version of Simulated Binary crossover (SBC) crossover (Deb et al. 1995) for our case through which values of   ,   and   are generated by where () and u is a random number in the range ∈ [0, 1] If a parameter exceeds the upper or lower bounds of the prior, it is clipped to be within the permissible range.The   is a hyperparameter which can be tuned to improve the performance, but we fix  c to 2. The motivation for using SBC comes from its ability to consistently generate unique offspring from the same set of parents.
For mutation, we employ the Non Uniform Mutation (NUM) operator (Michalewicz 1996).Specifically, we randomly select a feature, say   (), from the individual and tweak one of its parameters.We allow the mutation operator to alter the amplitude (  ), width (  ) or position (  ) of the feature but leave the shape unaltered.Specifically, we draw two random numbers (, ) in the range ∈ [0, 1].The first random number () determines the sign, while the second one () determines the magnitude.
Here,   represents the parameter's original value in consideration, while  ′  denotes the mutated value.The term range() refers to the parameter's allowed range of prior values.The mutation index, denoted as   , is fixed at 2 for all further analyses.
We evolved our individuals for 800 generations since there was no significant improvement in the  2 during our initial runs.We reiterate the fact that GAs are stochastic in nature and do not guarantee finding the global optimum.The quality of the solution is contingent upon factors such as grammar, hyper-parameters, as well as the choice of crossover and mutation operators utilised.

Validation with Mocks
We borrow the mock dataset from the authors of Sohn et al. (2022) to validate our pipeline and optimise our hyperparameters.The mock data are generated by adding features to the mean of BOBYQA best-fit power law spectrum to CamSpec PR3 unbinned TT, TE andEE data in (30,2500), (30, 2000) and (30,2000) ℓ range respectively in addition to lowT (commander) and lowE (simall).These features were generated using five parameters which are not explicitly included in our grammar.The 100 realisations were generated assuming cosmic variance limited  2 distribution for   ℓ in ℓ < 30 range.For all other data, we draw  ℓ s from a multivariate normal distribution with the CamSpec covariance matrix.
These tests serve two purposes: 1. Compare the results of the GA to the known "true" values of the mock dataset to see how well the algorithm found the optimal solution.2. Test the robustness of the GA by introducing variations or noise into the mock dataset and seeing how well the algorithm is able to handle these variations.
Due to the computational complexity inherent in our pipeline, we were constrained to conduct tests on a limited scale; we assessed our methodology's performance using five realisations of each of the four mock features.In these tests, a single chain with an initial population of 100 individuals was run for 500 generations.The outcome of one of these runs is illustrated in Fig. 1.It is important to note that the mock features employed in these tests were generated using fiveparameter feature models, which were not explicitly included in our grammar.Nevertheless, our pipeline exhibited a broad capacity to detect these features, albeit with slight variations in amplitude and the range of wave numbers corresponding to the true underlying features.Notably, in three instances, the GA identified the presence of a solitary feature in the mock and utilised the null function from our grammar to replace the other terms in the feature generation equation, given by Eq. ( 4).Furthermore, we conducted an additional set of tests by running four supplementary chains with different seeds to investigate the impact of the initial population while fixing the mock dataset.GA consistently reproduced the shape and position of the features across all cases, demonstrating the robustness of our approach, as seen in Fig. 2.
Overfitting is a common pitfall where the model becomes too complex and excessively fits the noise in the data rather than the underlying signal.To mitigate this issue in the case of GA, we take several steps to mitigate the risk of overfitting.Firstly, we limit both the length and maximum number of generations of the GA to prevent extra tailoring of features.Additionally, we employ Bayesian Model selection after post-processing (see Appendix) to assess whether the data warrant additional degrees of freedom proposed by GA.Secondly, one can use  2 statistics using noise simulations of the data to quantify overfitting in terms of the probability to exceed (PTE) (Planck Collaboration 2020a) but this is computationally expensive with GA and out of scope for this study.
Our results using simulations indicate that we can recover various forms of local features in the primordial spectrum.However, to assign significance to such recovered features (addressing the issue of noise-fitting) and, in particular, dealing with the real data, we need to perform some comprehensive analysis that is beyond the scope of this work.We should clarify here that in this work, we are not performing any model selection or model comparison, and we employ our approach to provide valuable insights into credible and plausible forms of the primordial spectrum that are consistent with CMB observations.

Data
In this work, we utilise the latest CMB power spectra measurements from the Planck satellite, specifically the temperature and polarisation data sets provided by the Planck Collaboration (Planck Collab-  gen = 500.Individual coloured lines correspond to four independent runs with different random seeds for the initial conditions.The grey line shows the power law prior; the dashed black line shows the true feature used to create mock data.We find best-fit power spectra from different chains and qualitatively find true features irrespective of the initial population.oration 2020c).We make use of two likelihoods, namely Plik and CamSpec, to analyse the data.
(i) plik: We use the unbinned version of the plikTT-TEEE+lowl+lowE+lensing data set.The high-ℓ likelihood, referred to as Plik, provides power spectra in the temperature range 30 ≤ ℓ ≤ 2509 and polarisation range 30 ≤ ℓ ≤ 1997 by utilising cross-spectra between 100-, 143-, and 217-GHz maps Planck Collaboration (2020c).On the other hand, the low-ℓ likelihoods (Commander and SimAll) cover the multipole range 2 ≤ ℓ ≤ 30 for the TT and EE spectra.Additionally, the CMB lensing likelihood provides measurements of the lensing trispectrum, estimating the power of the lensing potential    ℓ over 8 ≤  ≤ 40.These likelihoods are accessed via the  library, which provides the log-likelihood values used in our analysis.We refer to this data combination as "Planck" throughout the paper.The background cosmological parameters, nuisance parameters, and the tensor-to-scalar ratio  are fixed to specific values obtained from the BOBYQA best-fit of plikTTTEEE+ lowl + lowE + lensing, with  set to 0.
(ii) CamSpec: In CamSpec analysis, we utilise the unbinned and co-added  ℓ s obtained from the CamSpec 12.6 version of the likelihood.CamSpec employs a correction procedure for each TE and EE cross-frequency spectrum, applying a fixed dust and temperature-topolarisation leakage template before co-adding the spectra to form the EE and TE components of its data vector, see Efstathiou & Gratton (2021).The CamSpec likelihood incorporates Planck's 100-, 143-, and 217-GHz maps, which were processed using the NPIPE pipeline.To further constrain the features in the low-frequency region (ℓ ∈ [2,29]), we include lowTT, lowTE and lowEE data from the Planck Legacy Archive (PLA) with a diagonal covariance matrix defined by uncertainties due to cosmic variance.The background cosmological and nuisance parameters are fixed to the values obtained from the BOBYQA best-fit of CamSpecTTTEEE+ lowl + lowE with  = 0, ensuring consistency across our analysis.

Planck Runs
We apply our method to Planck 2018 data using the plik likelihood as specified in the previous section.Fig. 3 shows our results for eight chains with an initial population of 100 after 800 generations.In the top-left Fig 3a, we show the fitness of the best individual (combined  2 from high-TTTEEE lowTT, lowEE and lensing likelihoods) as a function of generations.The horizontal black dotted and brown dotted lines are the best-fit of the power law model, and the Starobinsky model to Planck Data is provided as a reference.The Fig. 3b shows the best-fit features in () space, where the best-fit individual from each chain is plotted in solid bold lines, and other individuals in the final generation are plotted in grey lines.The best-fit power law is shown as a black dotted line.
To further investigate the spectra obtained from the GA, we plot the residual  ℓ with respect to the best-fit power law in Fig 3c; points with error bars are binned Planck data, whereas grey points represent unbinned Planck data.In the TT spectra, we see features produce considerable shifts in range (20 < ℓ < 500), whereas in TE and EE spectra, we see the GA fitting peaks and troughs in the range (300 < ℓ < 1200).In addition to the residuals, in Table 3, we show the breakdown improvement in  2 between different likelihoods and further breakdown high-TTTEEE into high-TT, high-TE and high-EE assuming non-diagonal components do not contribute significantly.All chains improve greatly upon lowTT and high-TTTEEE data but at a small cost of low-EE data.While all chains address the low-ℓ anomaly, we also note that some chains perform much better than other ones with respect to high-ℓ data.This motivates us to apply the GA pipeline to high-ℓ data only, a procedure which we present in greater detail in App. A.

CamSpec Runs
Similar to the previous section, we perform our analysis on CamSpec PR4 data by taking co-added    ℓ s and   ℓ s along with   ℓ obtained following the method of Planck Collaboration et al. (2016) for the maximum likelihood solution.The nuisance and background parameters, in addition to   and   , were fixed to the values taken from the CamSpecTTTEEE+lowl+simallE best-fit parameters.In Fig. 4a, we show the evolution of total  2 for eight chains through generations.We observe all eight chains outperform the power law shown in the black dotted line.Next, Fig. 4b shows best-fit individuals from each chain in the () space.The GA finds considerable deviations from the best-fit power-law spectrum in the ℓ < 1000 region.
In Fig. 4c, we plot the  ℓ residuals for TT, TE and EE spectra, which makes it apparent where a significant share of improvement in  2 is coming from lowT and EE data.We also see our pipeline   trying to fit the peak in EE residual space around ℓ ≈ 400.In Table 4, we quantify the improvement in  2 from the best-fit individuals.All eight chains have more than 25  2 improvements, with the best one Δ 2 ≈ −35 compared to the fiducial power law model.The analytically best-fit spectrum obtained from the GA with the CamSpec run is given by:  Our analysis reveals that in both the CamSpec and Planck runs, a significant portion of the improvement is attributed to fitting the lower power in the spectra around ℓ ≈ 20 along with the small periodic deviations in the EE spectra at higher multipoles in CamSpec runs, which further improve the fit to the data as compared to the Planck runs.However, we do not observe the features obtained by the Planck run in the TT spectra within the multipole range of ℓ = 750 − 1000 in CamSpec runs.The global feature found with the Planck run is also consistent with the first application of GA to inflation Kamerkar et al. (2022), although the best-fits found in this previous work did not address the low-ℓ anomaly at all.Furthermore, the oscillations in between / * = 10 −2 and 10 −1 of the best chain (red line in Figure 4c) are of similar frequencies and phase seen in the regularised MRL results (Sohn et al. 2022) on the CamSpec data.This concurrence highlights the robustness and reliability of our results and strengthens the validity of the observed features in the primordial power spectrum.

SUMMARY AND DISCUSSIONS
Our study introduces a novel GA methodology to effectively explore and identify features within the primordial power spectrum.By incorporating a high degree of freedom that can be interpreted analytically, our approach enables efficient and targeted searches for these features.Investigating features is essential as they can provide valuable insights into the underlying physics of the early Universe, potentially alleviate cosmological tensions, and help us validate or rule out inflationary models.
Through tests on mock data with known true features, we evaluated the performance of the GA pipeline.Despite not explicitly incorporating the specific five-parameter featured models used to generate the mock data in our grammar, the GA demonstrated its ability to broadly detect these features, although with some amplitude and -range shifts.Further refinement is needed in the pipeline, particularly in improving the efficiency of computing CMB angular power spectra (C ℓ ) and its likelihood, to enable more accurate and comprehensive feature searches by incorporating variations in background and nuisance parameters.
Our study showcases the effectiveness of genetic algorithms in directly identifying features in the primordial power spectra using Cosmic Microwave Background observations.By applying our pipeline to Planck and CamSpec PR4 data, we identify potential features using a simple motivated grammar made of both global and local features.These features may be due, e.g., to transient deviations from the slow-roll dynamics during inflation or persistent small features atop a smooth primordial scalar potential.Another advantage of GA lies in their ability to provide best-fit analytical solutions that consist of combinations of individual features.This allows us to evaluate the effects of these features individually, helping us to identify the source of improvement in fitness to data more precisely than other non-analytic methods.
By disentangling the effects of these individual features, we can construct a model with fewer parameters.For instance, in one of our high-TTTEEE runs, which are further explained in the appendix, we discovered a feature that significantly improved fitness to data, reducing the  2 value by Δ 2 ≈ −12.This finding can help us design a featured model with just three parameters (A, C, D), whose significance can be assessed by employing an MCMC analysis and calculating the Bayesian evidence.This quantification will provide additional insights into the reliability and robustness of the features found by GA.We should emphasise here that in this work, we are not performing any model selection analysis, and we use GA mainly to suggest some plausible and viable form of the primordial spectrum that can fit CMB observations well.In our future work, we plan to extend our research by incorporating more complex grammatical structures and evaluating their effectiveness on various CMB surveys, including SPT (Benson et al. 2014), ACT (Aiola et al. 2020), and their combinations with Planck data.Moreover, we expect the straightforward adaptability of our methodology to forthcoming survey data, such as the Simons Observatory (SO) (Ade et al. 2019) and CMB-S4 (Abazajian et al. 2016).
The regime of the applicability of the techniques developed in this work goes well beyond testing primordial features with current and future CMB probes.Upcoming galaxy surveys such as DESI (Levi et al. 2013), EUCLID (Laureĳs et al. 2011), andLSST (Ivezić et al. 2019), in addition to 21-cm tomographic observations, exhibit raw statistical power comparable to CMB measurements and would provide valuable complementary constraints (Chen et al. 2016a,b;Palma et al. 2018;Chandra & Pal 2022;Chandra 2022).From a theoretical standpoint, it is crucial to emphasise the potential of our analysis methodology to explore deviations from scale-invariance, particularly in the tensor power spectrum.Our approach facilitates a model-independent investigation of such deviations through the primordial B-modes, with consideration for their potential detection in the forthcoming generation of CMB experiments, such as SO, CMB-S4, and LiteBIRD.One could also naturally expand the scope of our analysis to include features atop a GW signal at different scales, such as the one recently detected by pulsar timing array experiments (Agazie et al. 2023;Antoniadis et al. 2023;Reardon et al. 2023;Xu et al. 2023).Our approach to features would be equally effective also in testing primordial physics with future gravitational wave detectors at small scales.

Figure 1 .
Figure1.Samples from the final generation of GA explorations around power law with  pop = 100 evolved till  gen = 500 for four different feature mocks.The thick grey line shows the power law prior; the dashed black line shows the true feature used to create mock data, and the solid colour line shows the best-fit qualitatively detects the true feature in the mock.Additionally, silver lines depict the population of the final generation.

Figure 2 .
Figure 2. Best-fit power spectra for a featured mock with  pop = 100 and Evolution of the  2 as a function of the generations.
from the best-fit Planck Data.

Figure 3 .
Figure3.Results from the GA search around a power-law power spectrum with Planck data after 800 generations with a population of 100.The individual coloured lines correspond to eight independent runs with different random seeds for the initial conditions.In the top left, we show the  2 evolution of the best-fit individual in the population as a function of the number of generations.The best-fit  2 from the power law and Starobinsky model are shown with dotted black and brown lines, respectively.We plot the best-fit power spectra from different chains in the top-right.The grey lines depict the rest of the population at the final generation, visually representing the distribution of individuals at the end of the evolutionary process.Finally in bottom subplot, we illustrate the residuals in   ℓ ,    ℓ and   ℓ .The binned and unbinned data are shown in red and grey , respectively.

Figure 4 .
Figure 4. Results from the GA search using CamSpec data after 800 generations with a population of 100.Each solid-coloured line corresponds to one of the eight independent runs with different random seeds for the initial conditions.The top left plot displays the evolution of the  2 value for the best-fit individual in the population as the number of generations increases.The dotted black line represents the best-fit  2 from the power law model.In the top-right plot, we showcase the best-fit power spectra obtained from different chains and grey lines illustrating the population distribution in the final generation.In the bottom plot, we illustrate the residuals in   ℓ ,    ℓ , and   ℓ are illustrated.

Table 2 .
Overview of the grammar used in GA analysis, including global and local classes, along with the corresponding functions associated with each class.The global class encompasses trigonometric functions and polynomials, while the local class comprises various combinations of trigonometric, polynomial, and higher-order trigonometric functions.

Table 3 .
Summary of the best Δ 2 ( =  2 GA - 2 prior ) values for different Planck 2018 likelihoods across eight chains.The best chain has a Δ 2 value of -21.356, indicating a significant improvement in the fit high-TTTEEE and lowl although at a small cost of lowE

Table 4 .
Summary of the best Δ 2 ( =  2 GA - 2 prior ) values for co-added CamSpec data across eight chains.All eight chains have a Δ 2 value < -25, with best-fit improving  2 by 35.296, signifying a notable improvement compared to the power law model.Agencia Estatal de Investigación) through the Grant IFT Centro de Excelencia Severo Ochoa No CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033.AS would like to acknowledge the support by National Research Foundation of Korea 2021M3F7A1082056, and the support of the Korea Institute for Advanced Study (KIAS) grant funded by the government of Korea.A.S. would like to thank IFT for the hospitality during the early stages of this work.