A double take on early and interacting dark energy from JWST

The very first light captured by the James Webb Space Telescope (JWST) revealed a population of galaxies at very high redshifts more massive than expected in the canonical ΛCDM model of structure formation. Barring, among others, a systematic origin of the issue, in this paper, we test alternative cosmological perturbation histories. We argue that models with a larger matter component Ω m and/or a larger scalar spectral index ns can substantially improve the fit to JWST measurements. In this regard, phenomenological extensions related to the dark energy sector of the theory are appealing alternatives, with Early Dark Energy emerging as an excellent candidate to explain (at least in part) the unexpected JWST preference for larger stellar mass densities. Conversely, Interacting Dark Energy models, despite producing higher values of matter clustering parameters such as σ 8, are generally disfavored by JWST measurements. This is due to the energy-momentum flow from the dark matter to the dark energy sector, implying a smaller matter energy density. Upcoming observations may either strengthen the evidence or falsify some of these appealing phenomenological alternatives to the simplest ΛCDM picture.


Introduction
The minimal ΛCDM model of cosmology, described by only six fundamental parameters, has successfully explained a large number of cosmological observations at different scales.Nevertheless, in recent years, several anomalies have emerged, challenging the model across all cosmological epochs.
Currently, the most significant problem -known as the Hubble tension [1][2][3] -is represented by a mismatch at the level of ∼ 5σ between the value of the present-day expansion rate of the Universe (H 0 ) inferred from Planck-2018 CMB observations (H 0 = 67.4± 0.5 km s −1 Mpc −1 [4]) and the value of the same parameter measured directly from local distance ladder measurements using Type-Ia Supernovae calibrated with Cepheid variable stars (H 0 = 73 ± 1 km s −1 Mpc −1 [5]).
Another conundrum, albeit less significant than the Hubble tension, is the so-called S 8 tension.This parameter is closely related to the clustering of matter in the Universe, and its value can be inferred both from the measurements of CMB anisotropies as those from Planck, and -more directly -from the measurements of galaxy lensing made by experiments such as the Dark Energy Survey (DES) [6,7], and the Kilo-Degree Survey (KiDS) [8,9].While Planck data favors larger values of S 8 , KiDS and DES seem to prefer lower ones, leading to a well-documented discrepancy ranging between 2 and 3 standard deviations [10]. 1  Yet another emerging anomaly -which is the focal point of this study -pertains to the observations recently released by the James Webb Space Telescope (JWST) of galaxies at very high redshifts (much higher than those commonly explored by large-scale structure galaxy surveys).Interestingly, the preliminary JWST results point towards a population of surprisingly massive galaxy candidates (see e.g.[12][13][14][15][16][17][18][19]) with stellar masses of the order of M ≥ 10 10.5 M ⊙ .In a recent study [20], it was pointed out that the JWST data indicates a higher cumulative stellar mass density in the redshift range 7 < z < 11 than predicted by the ΛCDM model, questioning one more time the canonical cosmological picture [21][22][23].
Despite the fact that recent comparisons of photometric and spectroscopic redshifts in overlapping samples of galaxies solidify the evidence for a high space density of bright galaxies at z ≳ 8 compared to theoretical model predictions [24], the JWST anomalous observations could still hint a lack of accuracy in extracting the intrinsic galaxy properties [25][26][27][28][29][30][31].Other possibilities recently explored include the hypothesis of unusual dense regions due to the currently limited JWST observations (which only cover an area of approximately 38 square arcminutes), as well as the possible presence of unknown systematics in CMB Planck polarization data [32,33].
Although it is certainly premature to draw any definitive conclusions from these preliminary observations, if neither of the aforementioned possibilities can account for the discrepancies between the JWST results and the theoretical predictions of a baseline ΛCDM cosmology, it may be necessary to consider modifications to the model itself  or to the galaxy formation process [31,[56][57][58][59].In this work, we take a step forward in this direction by testing alternative models where galaxy evolution could be notably different than in ΛCDM.
In particular, we consider extensions related to the dark energy sector of the theory as possible phenomenological alternatives to explain the JWST preliminary findings.We test Early Dark Energy (EDE) and Interacting Dark Energy (IDE) cosmologies as both these extended scenarios, featuring modifications in the growth of structure, might predict a different evolution of perturbations, potentially resulting in the formation of more massive galaxies.We demonstrate that while EDE emerges as an excellent candidate to explain (at least partially) the unexpected JWST preference for larger stellar mass densities, IDE is generally disfavored by JWST measurements, despite yielding higher values of matter clustering parameters σ 8 and S 8 .
The structure of the paper is as follows: in section 2 we further elaborate on the reasons why we choose Early Dark Energy and Interacting Dark Energy scenarios and also we define the physics and the basic equations governing the cosmological evolution in these models.In section 3 we describe the main observables considered in this work, namely, the cumulative stellar mass density, the observational data and the likelihoods exploited in the numerical analyses.section 4 contains our main results.We conclude in section 5.

Extended Dark Energy scenarios
Despite the undeniable uncertainties surrounding the preliminary findings from JWST, one might wonder whether these new emerging anomalies could be somehow linked to other wellknown longstanding problems in cosmology, such as the Hubble tension.This raises the question of whether they could both originate from a common issue related to our current theoretical understanding of the Universe and, on a broader scale, what kind of beyond-ΛCDM phenomenology (if any) can increase the present-day expansion rate of the Universe while also leading to a higher cumulative stellar mass density at earlier times.
At first glance, this question can even appear misplaced, as these observations are often believed to imply an older universe compared to the ΛCDM predictions.Since the age of the Universe is roughly ∝ 1/H 0 , this would suggest that increasing the Hubble constant could worsen the discrepancy with the observations released by JWST.However, structure formation is influenced by the evolution of primordial density perturbations and the underlying cosmology.Models addressing the Hubble tension often propose modifications at both the background level and in perturbation dynamics.This could allow for comparable or greater structure formation in a younger Universe (various N-body simulations of extension to ΛCDM show such variations, e.g.[60][61][62][63][64][65]). Furthermore, in beyond-ΛCDM models different correlations among cosmological parameters can shift their fitting values.These effects may significantly impact parameters related to structure formation, such as the matter density Ω m , and the other matter clustering parameters σ 8 and S 8 = σ 8 • (Ω m /0.3) 1/2 .These correlations are crucial as they could affect the amplitude and shape of the matter power spectrum.
An exercise certainly useful for understanding which kind of phenomenology could hit two targets with one arrow -increasing H 0 and aligning more closely with JWST -is breaking down the problem into smaller parts.In particular, we focus on the baseline ΛCDM model fixing all parameters to the best-fit values provided by Planck and altering each one individually within a 4-standard-deviation range.Through this analysis, focusing on the JWST likelihood (χ 2 JWST ), we identify physical adjustments for better consistency with observations.From Figure 1 we can derive a quite significant amount of information: • First and foremost, we observe that the parameter on which χ 2 JWST is most sensitive is the matter density parameter Ω m .In particular, a larger fraction of matter in the Universe will considerably improve the quality of the fit to JWST observations by facilitating structure formation. 2 • Secondly, we can clearly note that increasing the amplitude of the primordial perturbations A s or considering a larger tilt n s results in a significant reduction in χ 2 JWST .This is in line with previous findings documented in Ref. [32], where it was argued that relaxing the Planck constraints on polarization, which in turns allows τ to reach considerably higher values [33], can substantially improve the agreement between JWST and CMB data3 .Similarly, larger n s can substantially increase the power in the matter power spectrum on small scales, also facilitating more structure to form.
• Finally, concerning H 0 , we observe that a straightforward increase in the value of this parameter while simultaneously keeping σ 8 constant worsens the fit to JWST data.This aligns with the argument that, when fixing structure formation parameters, a younger Universe reduces the number of structures that can form.However, the impact of H 0 on χ 2 JWST is relatively small and it proves how the Hubble constant plays only a partial role in a more complex interplay of various parameters.Summarizing these results, from a phenomenological standpoint, an effective model to increase the value of H 0 and improve the agreement with preliminary JWST data should predict a higher spectral index, along with a greater quantity of matter in the Universe and possibly higher values of σ 8 and S 8 .On one side, this phenomenology is common in proposals aimed at resolving the Hubble tension by introducing new physical components that act before recombination.For this reason, we explore extensions related to the early Universe and, as a case study, analyze EDE cosmology.On the other hand, larger values of matter clustering parameters σ 8 and S 8 can also be achieved within late-time solutions of the Hubble tension that attempt to modify physics after recombination, influencing the value of H 0 derived from the angular distance from the CMB.Therefore, in the spirit of not leaving anything untried, we also test IDE cosmologies where both the growth of perturbations and the matter clustering are significantly different than in ΛCDM.Below, we summarize the basic aspects of the theoretical models we will be considering.

Early Dark Energy
Early Dark Energy models are a natural hypothesis of dark energy, see e.g., Refs.[68][69][70][71][72]. Deviating from the traditional cosmological constant framework, EDE models account for a non-negligible contribution from dark energy in the early Universe.In addition, these EDE models can be based on a generic dark energy fluids which are inhomogeneous.Their density and pressure vary over time, leading to a non-static equation of state.The phenomenological analyses of these inhomogeneous dark energy models usually require additional dark energy clustering parameters, the dark energy effective sound speed and the dark energy anisotropic stress.The effective sound speed determines the clustering properties of dark energy and consequently it affects the growth of matter density fluctuations.Therefore, in principle, its presence could be revealed in large scale structure observations.The growth of perturbations can also be affected by the anisotropic stress contributions which lead to a damping in the velocity perturbations.
Recently, EDE models have garnered significant attention, particularly due to their potential role in addressing some of the aforementioned cosmological tensions [75,91,100].Our analysis will concentrate on the EDE implementation detailed in [84].This model proposes that, in the early Universe, a light scalar field deviates from its potential minimum and, constrained by Hubble friction, is functionally similar to a cosmological constant.As soon as, at some particular redshift z ⋆ , the Hubble parameter reduces to be less than the mass of the field, the scalar field rolls down its potential and begins to oscillate about the minimum.To avoid spoiling late-time cosmology, the vacuum energy must redshift away quicker than matter (i.e., faster than a −3 ), and the field should behave as a subdominant component.A typical set of parameters used in this model is: the fractional contribution to the total energy density of the Universe, f EDE (z) ≡ ρ EDE (z)/ρ tot (z) evaluated at the critical redshift z c at which it reaches the maximum value, and θ i , which is the parameter that usually describes the initial field displacement.This particular behavior implies a larger amount of energy-density in the early Universe (just prior to recombination), a reduction of the sound horizon and, consequently, a larger value of the Hubble constant inferred by CMB observations.This is the reason why EDE models have been proposed as a possible solution to the Hubble constant tension.

Interacting Dark Energy
Interacting Dark Energy models describe a phenomenological scenario where the dark fluids of the Universe interact with each other by allowing a transfer of energy and/or momentum between them, see e.g., Refs. .Instead, the other components of the Universe (such as radiation and baryons) remain unaffected.The background evolution for dark energy and dark matter is modified, as the continuity equations for the single component present an interaction function Q whose sign governs the energy-momentum flow.A negative value of the interaction rate, Q < 0, implies a transfer of energy and/or momentum from pressureless dark matter to dark energy, while the opposite, refers to an energy-momentum flow from the dark energy sector to the dark matter one.In order to solve the background evolution, one would need a specific interaction function Q. Depending on such a function, it can be solved either analytically or numerically, together with the equation for the Hubble rate evolution.In what follows we shall use the well-known interaction rate [101,102,[156][157][158]: where ξ is a dimensionless coupling parameter.The equation governing the evolution of the density perturbations for the dark sector can be found in [101,102,157].IDE models may suffer from instabilities in the perturbation evolution [101,156].Our analysis adheres to the criterion of Ref. [102] in terms of the so-called doom factor which is required to be negative in order to avoid instabilities.Consequently this stability condition for our case is translated into a stable parameter space in which (1 + w) and ξ must have opposite sign [102].Therefore, in the phantom regime in which (1 + w) is a negative quantity the dimensionless coupling ξ must be positive.On the other hand, in the quintessence region ξ must be negative.For earlier studies, see [101,116,156,157,[159][160][161][162][163][164][165][166].
We conclude this section with a final remark: even if the interaction scenario considered here is a pure phenomenological model, some studies have shown that using a multi-scalar field action, the coupling function (2.1) can be derived [167].Therefore, the interaction model of the form given by Eq. (2.1) also benefits from a solid theoretical motivation.

Theory
To explore the extended dark energy scenarios in relation to the JWST observations, we strictly follow the methodology detailed in Ref. [32].We compute the predicted Cumulative Stellar Mass Density (CSMD hereafter), which is given by where we have is the cosmic baryon fraction and ϵ is the star formation efficiency of converting baryons into stars.For our analysis, we opt for a conservative approach and fix ϵ = 0.2 following Ref.[168].However, as suggested by Ref. [168], in principle star formation efficiency can be a function of the halo mass and further adjustments to star formation physics might be needed for more precise computations [169].Nonetheless, for the mass scale we are working with, ϵ can only vary smoothly as a function of mass, following a power law.Therefore, given the short mass range we are using in our analyses, this does not change in a significant way our main conclusions.Notice also that we consider the cosmic baryon fraction instead of computing the baryon evolution in different halos [64,[170][171][172][173][174][175] 4 .All these methodological choices and simplifications are widely used in the literature and allow us to present conservative and credible results that, without sacrificing generality, can be directly compared with similar works following the same approach.Furthermore, dn(M,z) dM is the comoving number density of the collapsed objects between a mass range from M to M + ∆M at a certain redshift, given by where ρ m is the comoving matter density of the Universe, F (ν) the Sheth-Tormen (ST) Halo mass function [176,177] and σ(R, z) is the variance of linear density field smoothed at a scale R, where we assume a Top-Hat window function.The ST formalism is theoretically motivated in terms of the collapse of halos [178,179] and has been exhaustively tested by N-body simulations for different dark energy models taking different priors for Ω m and Ω Λ [180].In this study, we assume that the halo mass function of the ST type applies to both theoretical frameworks examined.Nevertheless, analyses of future JWST data may need to consider a more sophisticated fit.

Observational data
Regarding the JWST observations, given their preliminary nature and the wide research interest they have generated, a multitude of datasets obtained following different methodologies have been released.Consequently, a choice of which dataset to use becomes necessary.In this study, we consider four independent datasets summarized in Table 1 and listed below: Notice that many additional measurements have been released, see, e.g., Ref. [26,181].In future, these measurements could potentially serve as independent tests of our preliminary findings.
• Confirmed CSMD measurements taken from Ref. [20].We find this dataset particularly well-suited for our analysis as it provides explicit constraints on the CSMD that are directly linked to cosmological structure formation.For this dataset, the errors reported in Table 1 are assumed to follow a Log Normal distribution.We refer to this dataset as JWST-CEERS .
• Five observational CSMD coming from the photometric data of JWST coverage of the UKIDSS Ultra Deep Survey (UDS) and Hubble Ultra Deep Field (HUDF) [182].The errors for HUDF & UDS displayed in Table 1 are conservatively taken as the maximum value, while the stellar masses are set to be equal to M = 10 8 M ⊙ as suggested in Table 6 of Ref. [182].We refer to this dataset as JWST-HUDF &UDS .
• Two values of the observational CSMD coming from optical data in the GLASS-ERS 1324 program [12].The stellar mass for GLASS datapoints is set to the average between the mass interval reported in Table 1 where we refer to this dataset as JWST-GLASS .
• JWST FRESCO NIRCam/grism survey [26].This dataset spans an area of 124 arcmin 2 , covering a survey volume of approximately 1.2 × 10 6 Mpc 3 within the redshift range z ∈ [5,9].We refer to it as JWST-FRESCO.Taking the JWST-FRESCO data at face value (and barring any potential systematic errors), we consider 3 obscured galaxies located within densely dusty regions, with redshifts in the range 5 ≲ z ≲ 6.
Referring to Fig. 3 of Ref. [26], we can see that these galaxies show exceptionally extreme properties such as dark matter halo masses of log (M halo /M ⊙ ) = 12.88 +0.11 −0.13 , 12.68 +0. 23 −0.17 , and 12.54 +0.17 −0.18 .Notice also that the quantity measured by JWST-FRESCO is the cumulative comoving number density of dark matter halos, not the CSMD.For any given model of cosmology, the cumulative comoving number density of dark matter halos can be computed as: where the integrand in Eq. (3.3) is defined by Eq. (3.2).Regarding the error for the FRESCO dataset, we use the approach of Poissonian approximation for small numbers of observed events (the interested reader can refer to our Appendix A for more details).
We highlight that in Eq. (3.1), represents the modeldependent comoving volume of the Universe between redshifts z 1 and z 2 (with R(z) being the comoving radius at z).In the values of ρ obs given in Table 1, the comoving volume has been derived assuming the best-fit ΛCDM model (Planck TTTEEE+lowE+lensing CMB measurements with h = 0.6732, Ω m = 0.3158, n s = 0.96605, σ 8 = 0.8120, see [20]), making the data points model-dependent, too.As a result, these values need to be appropriately rescaled before interfacing them with the predicted values for dark energy scenarios.That is, they must be corrected by a factor equal to the ratio of the two comoving volumes (V Λ /V DE ).Similarly, the conversion between luminosity fluxes and distances for the stellar masses has been inferred, again, under the underlying assumption of the best-fit Planck ΛCDM model.A similar rescaling has to be applied, using the squared ratios of the luminosity distances.

Numerical analyses
For our analyses, we perform a Monte Carlo Markov Chain (MCMC) using the publicly available package Cobaya [183] and generate theoretical predictions exploiting with a modified    version of the software CLASS [184,185] to address the IDE scenario, while, for EDE, the publicly available software CLASS EDE, 5 which solves the evolution of the background and of the perturbations in the presence of a scalar field by means of the Klein-Gordon equation.We investigate the posterior distributions of our parameter space through the MCMC sampler developed for CosmoMC [186,187] and tailored for parameter spaces with a speed hierarchy which also implements the "fast dragging" procedure [188].
The likelihood used for the MCMC analysis are: • CMB temperature and polarization power spectra from the legacy Planck release [4,189] with plik TTTEEE+low-ℓ+lowE.
• Lensing Planck 2018 likelihood [190], reconstructed from measurements of the power spectrum of the lensing potential.
In the following discussion, we will refer to the combinations of these two datasets simply as CMB.The convergence of the chains obtained with this procedure is assessed using the Gelman-Rubin criterion [191] setting a convergence threshold at R − 1 ≲ 0.02.
Once the chains have converged, 6 for each sampled model we calculate the CMSD (or the comoving number density of dark matter halos for the JWST-FRESCO mesurements) corresponding to different parameter combinations explored in the MCMC analysis.In particular, we set the cosmology by varying the cosmological parameters {Ω b , Ω c , θ s , τ }, the inflationary parameters {n s , A s } and the extra parameters from EDE {f EDE , log 10 z c , θ i } and IDE {w, ξ}.Afterwards, CLASS is employed to compute the linear matter power spectrum P (k, z) and subsequently the variance and, taking its derivative, we estimate the Halo mass function, see Eq. (3.2).Eventually, through double integration over mass and redshift, Eq. (3.1), we arrive at ρ ⋆ (M ).For each point in the MCMC chains obtained within a given cosmological model, we calculate the χ 2 JWST for various datasets listed in Table 1.Each dataset corresponds to a distinct χ 2 JWST .Specifically, for any given dataset, χ 2 JWST is determined straightforwardly by where, xi is the observed value, σ 2 i is its error, x th i the theoretically predicted quantity and i runs over the number of data points of the specific JWST dataset.We then obtain the updated constraints on the cosmological parameters by re-weighting the MCMC chains, i.e. performing an importance sample, using the package getdist.

Results
We now discuss the results obtained for the various extended models of Dark Energy considered in this study.To promote better organization of our findings, we divide this section into two minor subsections: in subsection 4.1 we focus on the results obtained for EDE, while in subsection 4.2 we present the findings related to IDE.

Results for Early Dark Energy
We start by examining EDE.For this model, we summarize the best-fit values related to the most relevant parameters in tables 2 to 5, distinguishing the results obtained by considering the four JWST datasets independently.
First and foremost, we examine the best-fit values obtained by exclusively considering CMB measurements from the Planck satellite (indicated as CMB in the table).In this case, the values reported in tables 2 to 5, simply represent the combination of cosmological parameters for which χ 2 CMB (= 2772) acquires its minimum value among those obtained within the MCMC analysis.It is worth noting that we retrieve results widely documented in the literature.In particular, the Planck data, while not showing any substantial preference for a non-vanishing fraction of EDE, produce a best-fit value of f EDE = 0.06.This leads to a present-day expansion rate of the Universe H 0 = 69.45km/s/Mpc, which is generally higher than the best-fit value obtained for this parameter within the standard cosmological model.Another point that is worth emphasizing is that for the inflationary spectral index we get a best-fit value n s = 0.981, confirming once more the trend of EDE models in predicting a spectrum of primordial perturbations closer to the scale-invariant case than what is observed in ΛCDM and predicted by the most typical inflationary potentials [88,[192][193][194][195][196][197][198]     As a second step, following the methodology outlined in the previous section, for each combination of parameters in the MCMC chains (i.e., for each collected model), we calculate the χ 2 against the four different JWST datasets listed in Table 1.
We start discussing the results obtained re-weighting the chains in light of χ 2 JWST−CEERS resulting from the JWST-CEERS dataset.This dataset has been recently analyzed in many similar studies and allows us for direct comparison with the existing findings in the literature [21,35].In this case, we summarize the results in Table 2, distinguishing between the low (z low ) and high (z high ) redshift bins (also see Table 1).Similar to the CMB analysis, we present the specific combination of parameters that minimizes χ 2 JWST−CEERS .Furthermore, in Figure 2 we provide a triangular plot showing the distribution of sampled models and the correlations among different parameters, together with a color-map representing the value of JWST so that the color pattern in the figure represents the direction towards which one needs to move in the parameter space to improve the fit to JWST data.χ 2 JWST−CEERS .For the sake of comparison, in the same figure, we also depict the predictions of ΛCDM.A few intriguing conclusions can be derived from both Table 2 and Figure 2. Firstly, there are no significant differences between the results obtained for the high and low redshift bins.Secondly, as for the best-fit values of cosmological parameters, we now find a pronounced preference for a non-vanishing fraction of EDE, f EDE = 0.151.We also get higher σ 8 = 0.85 and observe the same trend towards higher values of inflationary spectral index n s = 0.997, now essentially consistent with a Harrison-Zel'dovich spectrum.As pointed out in section 2, this is exactly the kind of phenomenology one needs to increase the agreement with JWST data.Therefore -not surprisingly -the minimum value of significantly better than what we get in ΛCDM (where χ 2 JWST−CEERS ∼ 17 [32]).This suggests that EDE stands as a valid phenomenological alternative to explain (at least partially) the preliminary measurements released by JWST-CEERS.Furthermore, regarding H 0 , the JWST-CEERS best-fit value reads H 0 = 72.60 km/s/Mpc.Therefore, not only within the context of EDE we can improve the agreement between the theoretical predictions of the model and the JWST-CEERS data, but to achieve this, we move through the parameter space in the same direction needed to solve the Hubble tension, as well.This is also clearly confirmed by the color pattern in Figure 2, underscoring that it is indeed possible to address both issues within the same theoretical framework.
Finally, always in Table 2, we present the results inferred by summing up the χ 2 values of CMB and JWST-CEERS.We observe that the combination of parameters that minimizes the total χ 2 CMB+JWST is the same as that minimizing the fit to only the Planck data χ 2 CMB .At first glance, this implies that the cost of improving the fit to the JWST-CEERS data is an overall deterioration in the fit to the CMB.On the other hand, such deterioration is entirely expected, given the strong preference of Planck data for a ΛCDM cosmology and the general disagreement between JWST-CEERS data and the latter.As extensively documented in the literature and confirmed by our analysis, individual Planck data do not provide clear evidence in favor of an EDE cosmology.In any case, the best-fit parameters suggest an inclination towards a model where the fraction of EDE remains modest and well below f EDE ≲ 0.1.In contrast, reconciling JWST-CEERS with an EDE cosmology would require an EDE fraction f EDE ≳ 0.1.Forcing such an EDE fraction into the model would source significant effects in the CMB spectra that can only be partially compensated by the observed shift in the fitting values of other cosmological parameters.Just to provide an illustrative example, the increase in the expansion rate of the Universe before recombination due to a substantial EDE component leads to a significant reduction in the value of the sound horizon at the combination, forcing the value of H 0 in the direction of SH0ES, which is certainly not the direction favored by CMB data.Additionally, since EDE does not alter the physics of postrecombination, a higher H 0 implies a lower angular diameter distance from the CMB, D A .In turn, this leads to a shift in the wavenumber associated with the damping tail k D , as these two parameters are related by the relationship ℓ D ∼ k D D A , where the multipole ℓ D is also fixed by CMB measurements.In an attempt to maintain a good fit to the damping scale, the value of n s is shifted towards a scale-invariant primordial spectrum (i.e., n s → 1) which certainly improves the agreement between JWST-CEERS and EDE but is again highly disfavored by Planck (by over 8σ in ΛCDM).Overall, all these effects and shifts in the fitting values of cosmological parameters seem to favor JWST-CEERS observations.However, although they partially compensate for each other, they still remain somewhat disfavored based solely on CMB data, leading to a deterioration in the fit.
When analyzing the other JWST datasets listed in Table 1, all the conclusions we have drawn so far remain mostly true.For instance, by comparing the results obtained for JWST-CEERS in Figure 2 with those obtained for JWST-FRESCO in Figure 3, at first glance, we can spot the very same color patterns, indicating that a non-vanishing fraction of EDE could, in principle, help to reduce χ 2 FRESCO while also yielding higher values of H 0 .However, paying closer attention to the colorbar scale, we observe that as we approach the limit f EDE → 0 moving towards the ΛCDM cosmology, we get χ 2 JWST−FRESCO ∼ 50 for 3 data points.This value can be reduced all the way down to min(χ 2 JWST−FRESCO ) ∼ 21.42 when f EDE ∼ 0.15 and H 0 ∼ 73 km/s/Mpc, as seen in Table 5.On one hand, this lends weight to the idea that EDE could potentially pave the way to partially explaining more massive galaxies and higher values of H 0 .On the other hand, it is important to note a nearly threefold increase in χ 2 JWST−FRESCO compared to the results obtained for JWST-CEERS.Taking the large χ 2 at face value, we must draw the conclusion that the JWST-FRESCO dataset remains in strong disagreement with the theoretical predictions of the standard cosmological model and -to a lesser extent -with EDE as well.
Looking at the χ 2 of JWST-GLASS and JWST-HUD&UDS, a similar conclusion can be drawn.As a result, we conclude that while EDE represents a phenomenological possibility to partially address the JWST data, it falls short of being exhaustive in fully addressing issues, leaving the quest for a more comprehensive solution wide open.
Having that said, it is worth keeping in mind some caveats surrounding the joint analyses.
For instance, the total χ 2 CMB+JWST is obtained by considering the sum of χ 2 for each sampled model in the MCMC chains afterward and not through a joint analysis of the two experiments from the outset.Additionally, only CMB data are taken into account in the MCMC analysis, which we know do not favor high values of f EDE and H 0 .Considering other datasets, such as the measurements of the expansion rate provided by the SH0ES collaboration, could lead to significantly different results in terms of the χ 2 analysis, as typically pointed out by the EDE community (see, e.g., the discussion on page 25 of Ref. [75]).Hence, a full joint likelihood analysis of all these datasets (which is beyond the aim of this work) is needed before deriving any definitive conclusions on this matter.

Results for Interacting Dark Energy
We now move to the study of IDE.In this case, we consider three different models: the usual IDE cosmology with a fixed dark energy equation of state w ≃ −1, and wIDE models where the equation of state parameter w is free to vary, although limited either in quintessence (w > −1) or phantom (w < −1) regime.For the sake of simplicity, in this subsection, we present only the results obtained from JWST-CEERS.Several well-motivated reasons underpin this decision.Firstly, no significant differences emerged in the results for EDE when analyzing the four different JWST datasets listed in Table 1.Overall, all these observations converge on anomalous galaxies that are more massive than predicted by the standard cosmological model.Therefore, no significant disparities are anticipated when analyzing the same four datasets across the various IDE models proposed in this section.Yet another motivation involves noting that addressing these JWST anomalous observations requires a somewhat clear beyond-ΛCDM phenomenology that none of the three IDE models proposed here can offer.To streamline the analysis and emphasize this point, we focus on JWST-CEERS, which will provide the phenomenological guidelines applicable directly to all datasets not explicitly mentioned, without exception.
Table 6 displays the results for the IDE model with a fixed dark energy equation of state.Similar to EDE, we consider three different combinations of data: CMB, JWST-CEERS, and CMB+JWST-CEERS.In the table, we always show the combination of parameters that minimizes the χ 2 for these three datasets.When focusing solely on the Planck CMB data, we note that the best-fit value for the parameter encapsulating new physics -i.e., the coupling ξ -reads ξ = −0.28.This suggests a quite significant transfer of energy-momentum from the dark matter sector to the dark energy sector.As widely documented in the literature, such a transfer of energy-momentum leads to a higher present-day expansion rate of the Universe, whose best-fit value reads H 0 = 70.98 km/s/Mpc (significantly higher than in the standard cosmological model).Furthermore, while we do not observe significant differences in the value of the spectral index, we notice a tendency toward higher values of σ 8 = 1.09.This makes the model potentially interesting for JWST.Nevertheless, the results we obtain from JWST-CEERS data seem to indicate precisely the opposite.In contrast to the CMB fit (which prefers ξ < 0), when considering only the JWST-CEERS likelihood, the coupling parameter ξ tends towards ξ → 0. Consequently, we lose any ability to increase the presentday expansion rate, getting a best fit value H 0 = 67.43km/s/Mpc.Figure 4 further reinforces our conclusions: it is not n s but Ω m that now assumes a critical role.When moving in the direction ξ < 0, the matter density undergoes a significant decrease due to the energy transfer from dark matter to dark energy.This leads to a substantial increase in the value of σ 8 to compensate for the reduced Ω m .However, as illustrated by the color pattern in Figure 4, in order to minimize χ  preventing such energy transfer, essentially moving towards ξ → 0. This behaviour is further supported by comparing the best-fit values of Ω m (σ 8 ) for CMB and JWST: while in the former case Ω m = 0.214 (σ 8 = 1.09), for JWST-CEERS, we get back to more typical values Ω m = 0.301 (σ 8 = 0.878).Consequently, this model fails to provide a satisfactory fit to the JWST-CEERS observations.The best-fit values of cosmological parameters for the wIDE model with w > −1 (i.e., confined to the quintessence regime) are displayed in Table 7. Qualitatively, these results mirror those previously obtained for w = −1.There are no differences between the high and low redshift bins, and there is no preference for ξ ̸ = 0 from JWST-CEERS data.The behaviors of the parameters depicted in Figure 5 clearly indicate that introducing a coupling while leaving w free to vary does not improve the fit to JWST-CEERS observations.Once more, the reason behind this phenomenon is the decrease in matter density resulting from the energy flow within the interacting model.
The situation becomes somewhat more intricate when we turn to the case where the dark energy equation of state is confined to the phantom regime w < −1.In this case, the best-fit values of parameters are summarized in Table 8 for the usual combinations of datasets.When considering only the best-fit values from the CMB, we find the welldocumented Planck preference for a phantom equation of state [208,209], with the best-fit value reading w = −2.04.Given the well-known degeneracy between the effects produced by a phantom w and increasing the present-day expansion rate of the Universe, interacting phantom models can provide a much larger value of H 0 [147]   Universe, the question of whether (and to what extent) a phantom wIDE model could effectively contribute to explaining the anomalies observed in JWST remains a topic of debate.
Looking at the brighter side, when we consider the fit to JWST-CEERS data, we observed that both for the high (χ 2 JWST−CEERS = 16.3) and low (χ 2 JWST−CEERS = 10.8)redshift bins, χ 2 JWST−CEERS slightly improves compared to the standard cosmological model.Furthermore, the JWST-CEERS data seem to point toward a non-zero coupling ξ ∼ 0.2 − 0.3, which is substantially higher than that preferred by the Planck data (whose best-fit value is ξ ≃ 0.05).The reason underlying this preference is that to ensure the stability of the perturbations, ξ is now required to be positive.This results in a shift of the energy flow from the dark energy sector to the dark matter one.Consequently, increasing the value of the coupling means injecting more power into the matter sector, facilitating structure formation, and improving the JWST-CEERS fit.However, looking at Figure 6 -where, as usual, we plot the correlations among different parameters together with a color-map representing the value of χ 2 JWST−CEERS -it becomes really difficult to identify a color pattern representing the direction in which we need to move in the parameter space to improve the fit to JWST-CEERS data.Despite this, with a good degree of imagination, we can speculate that by moving towards larger values of the coupling ξ → 0.4 (corresponding to values of Ω m → 0.4), the value of χ 2 JWST−CEERS seems to undergo a general improvement for the highlighted reasons, and its value becomes as good as the ΛCDM one.Interestingly, looking again at Figure 6, we note that the value of the Hubble rate corresponding to such a coupling is H 0 ∼ 70 km/s/Mpc, i.e., close to the result provided by SH0ES collaboration.
Therefore, in the context of IDE cosmology, the only potential scenario in which we can simultaneously slightly improve the agreement between the model predictions and JWST data while obtaining values of H 0 in line with local distance ladder measurements seems to involve considering a phantom component of the equation of state of dark energy.said, the ability of this model to address these two issues remains somewhat limited, above all when compared to the competing EDE solutions discussed in the previous section.

Discussion and Conclusions
Cosmological data have entered an era of high precision.Advanced telescopes with higher sensitivity are unveiling measurements at previously uncharted cosmological scales and epochs.However, this is a double-edged sword: while high-precision parameter extraction became achievable, a range of unexplained anomalies emerged, some of them growing in significance rather than diluting as the error bars decreased.The most recent anomaly involves the JWST observations of a population of very highmass galaxies at previously unexplored redshifts of z ∼ 7 − 10.This suggests a higher cumulative stellar mass density in the redshift range 7 < z < 11 than predicted by the standard ΛCDM model of cosmology.Assuming no systematic issue behind these preliminary findings, one might wonder whether these emerging anomalies could be somehow connected to other long-standing cosmological puzzles, such as the Hubble tension.All these issues might hint at a shared limitation in our current comprehension of the Universe, motivating the need to consider alternative theoretical scenarios.
In this paper, taking the four independent JWST datasets at face value, we first explore the correlation between the JWST likelihood and the fundamental six-ΛCDM parameters.We argue that, from a phenomenological standpoint, models with a larger matter component Ω m , a higher amplitude of primordial inflationary fluctuations A s and a bigger scalar spectral index n s are able to predict larger cumulative stellar mass densities, thus providing a better fit to JWST data.As a next step, we notice that part of this phenomenology is very common in models featuring new physics in the dark energy sector (both at early and late times) that have been recently proposed as possible solutions to the Hubble tension.In such models the behavior of the dark energy fluid can also influence the growth of structure, potentially leading to the formation of more massive galaxies able to account for the preliminary JWST measurements.Inspired by this idea, we explore whether in the context of Early Dark Energy or Interacting Dark Energy, the JWST findings could be explained, or at least, if within these scenarios, the fit to JWST observations is improved with respect to the one in the minimal ΛCDM framework.
We find that EDE, which leads to an increased Ω m and n s concurrent with a rising EDE fraction f EDE , constitutes an excellent candidate.Not only we can improve the agreement between the theoretical predictions of the model and the JWST data (i.e., the minimum value of χ 2 JWST obtained for both the high and low redshift bins is significantly better than what we get in ΛCDM), but to achieve this, we move through the parameter space in the same direction needed to solve the Hubble tension.This underscores that it is indeed possible to address both issues within the same theoretical framework.
Conversely, wIDE models featuring a dark energy equation of state w ≥ −1 are generally disfavored from JWST, despite yielding higher values of matter clustering parameter σ 8 .This is due to the energy flow from the dark matter sector to the dark energy one, implying a smaller Ω m .On the other hand, when the equation of state is confined to the phantom regime w < −1, the situation becomes somewhat more intricate.Whether, and to what extent, a phantom wIDE model could effectively contribute to explaining the anomalies observed in JWST findings remains an open question.The energy-momentum dynamics and parameter degeneracies can lead to a significant increase in the matter component, which in turn can slightly improve the agreement between the model predictions and the JWST data while also yielding values of H 0 in line with local distance ladder measurements.However, the ability of this model to address these two issues simultaneously remains limited compared to EDE.
We conclude by pointing out that other studies have recently investigated dark energy and dark matter alternatives to reconcile the JWST anomalies.Just to mention a few appealing possibilities, in Ref. [35] it was argued that adjusting the dark energy sector to accommodate a dynamic equation of state can be compatible with existing observations.However, this would require considering somewhat exotic dynamical dark energy models that necessitate specific configurations of dark energy parameters.Always in relation to the dark energy sector of the cosmological model, in Refs.[34,47], it was argued that an evolving DE component with positive energy density on top of a negative cosmological constant can be consistent with JWST observations.Other different yet interesting possibilities explored to approach the structure formation dilemma involve considering extensions related to particle physics, such as Fuzzy Dark Matter models consisting of ultra-light axions [37] or Warm Dark Matter [39,40].Axion models have the potential to mitigate the formation of smaller structures, aligning JWST's stellar mass density data with the CMB optical depth near z ≈ 8. Instead, Warm Dark Matter models introduce a structural formation delay, predicting fewer low-mass systems at higher redshifts.However, it is worth noting that up until now, observations do not conflict with predictions from either cold dark matter or Warm Dark Matter theories for particles heavier than 2 keV [39].Finally, another category of solutions competitive with those examined in this study frequently proposes new physics in the gravitational sector.For instance, if long-range attractive forces stronger than gravity are realized in nature, the formation of cosmic halos could begin during the radiation-dominated era, providing a possible explanation for the excess of heavy galaxies observed by JWST.
Future JWST data will either reinforce or diminish the need for exploring physics beyond where, as usual, S ≃ 1.645 and β(S), γ(S) and δ(S) are given by Eqs.( 9), ( 10) and ( 12) of Ref. [211], respectively.We stress that this statistical methodology, while accurate, necessarily introduces an additional layer of approximation.For this reason, we remain conservative proving the uncertainties on log n(> M halo ) at 95% CL (corresponding to fixing S ≃ 1.645).

Figure 1 .
Figure 1.Changes in χ 2 JWST when varying a single parameter while keeping the others fixed to the Planck ΛCDM best-fit values (bold points in the figure): Ω b h 2 = 0.02238, Ω c h 2 = 0.1201, H 0 = 67.32km/s/Mpc, n s = 0.9659, A s × 10 9 = 2.1, τ = 0.0543, and Ω m h 2 = 0.143.In the top panel's third plot, when H 0 is free to vary, σ 8 is kept fixed by rescaling A s accordingly.

Figure 2 .
Figure 2. Triangular plot showing the distribution of points and the correlations among the most relevant parameters of EDE.The color map refers to the value of χ 2JWST so that the color pattern in the figure represents the direction towards which one needs to move in the parameter space to improve the fit to JWST data.

Figure 3 .
Figure 3. Triangular plot showing the distribution of points and the correlations among the most relevant parameters of EDE.The color map refers to the value of χ 2FRESCO so that the color pattern in the figure represents the direction towards which one needs to move in the parameter space to improve the fit to JWST FRESCO data.

Figure 4 .
Figure 4. Triangular plot showing the distribution of points and the correlations among the most relevant parameters of IDE, when w is fixed to w ≃ −1.The color map refers to the value of χ 2 JWST so that the color pattern in the figure represents the direction towards which one needs to move in the parameter space to improve the fit to JWST data.

Figure 5 .
Figure 5. Triangular plot showing the distribution of points and the correlations among the most relevant parameters of wIDE, when w is free to vary in the quintessence region w > −1.The color map refers to the value of χ 2JWST so that the color pattern in the figure represents the direction towards which one needs to move in the parameter space to improve the fit to JWST data.

Figure 6 .
Figure 6.Triangular plot showing the distribution of points and the correlations among the most relevant parameters of wIDE, when w is free to vary in the phantom region w < −1.The color map refers to the value of χ 2 JWST so that the color pattern in the figure represents the direction towards which one needs to move in the parameter space to improve the fit to JWST data.

Table 1 .
Observational points for JWST for the four different dataset.The values here must be rescaled by the corresponding comoving volume and luminosity distance for the Planck bestfit ΛCDM model.

Table 3 .
Results for EDE.We provide the best-fit values of cosmological parameters, namely the combination that minimizes the χ 2 of the fit to the CMB data alone (χ 2

Table 4 .
Results for EDE.We provide the best-fit values of cosmological parameters, namely the combination that minimizes the χ 2 of the fit to the CMB data alone (χ 2 CMB ), JWST-GLASS data alone (χ 2 JWST−GLASS ), and CMB+JWST-GLASS data (χ 2 CMB+JWST ).
2JWST−CEERS , it becomes necessary to revert to the ΛCDM framework by

Table 6 .
Results for IDE, with the dark energy equation of state fixed to w ≃ −1.We provide the best-fit values of cosmological parameters, namely the combination that minimizes the χ 2 of the fit and, in fact, we obtain a best-fit value H 0 = 103.8km/s/Mpc.This essentially indicates that, without including datasets at lower redshifts, breaking this line of degeneracy proves challenging and values of H 0 in line with those measured by SH0ES collaboration can always be reintroduced by considering a sufficiently phantom dark energy component.In addition, due to a combination of correlations among different parameters such as H 0 , w, and Ω m , we observe that for the latter parameter, the best-fit value reads Ω m = 0.139 and is compensated by a substantial increase in σ 8 = 1.026.In light of these effects on parameters governing the matter clustering in the

Table 7 .
Results for wIDE, with w > −1 free to vary in the quintessence regime.We provide the best-fit values of cosmological parameters, namely the combination that minimizes the χ 2 of the fit That being

Table 8 .
Results for wIDE with w < −1 free to vary in the phantom regime.We provide the best-fit values of cosmological parameters, namely the combination that minimizes the χ 2 of the fit